From: bangerth Date: Mon, 8 Jul 2013 03:37:55 +0000 (+0000) Subject: Revert to a state of the test that really does fail. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7904f5a970672f09073352c325486c4856e6513f;p=dealii-svn.git Revert to a state of the test that really does fail. git-svn-id: https://svn.dealii.org/trunk@29947 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/base/work_stream_03.cc b/tests/base/work_stream_03.cc index 054ae23433..618700df8d 100644 --- a/tests/base/work_stream_03.cc +++ b/tests/base/work_stream_03.cc @@ -20,7 +20,6 @@ // verifying that it indeed works #include "../tests.h" -#include #include #include #include @@ -44,247 +43,134 @@ #include #include #include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include +#include +#include #include #include -#include - -using namespace dealii; - char logname[] = "work_stream_03/output"; - template class F : public Function { public: - F () - : - q(1) + F (const unsigned int q, + const unsigned int n_components) + : + Function(n_components), + q(q) {} virtual double value (const Point &p, - const unsigned int component) const + const unsigned int component) const { - Assert (component == 0, ExcInternalError()); - double val = 0; - for (unsigned int d=0; dn_components == 1), + ExcInternalError()); + double val = 0; + for (unsigned int d=0; d triangulation; -FE_Q<2> fe(3); -DoFHandler<2> dof_handler(triangulation); -ConstraintMatrix constraints; - -namespace -{ - template - struct Scratch - { - Scratch (const ::dealii::hp::FECollection &fe, - const UpdateFlags update_flags, - const ::dealii::hp::QCollection &quadrature) - : - fe_collection (fe), - quadrature_collection (quadrature), - x_fe_values (fe_collection, - quadrature_collection, - update_flags), - rhs_values(quadrature_collection.max_n_quadrature_points()), - update_flags (update_flags) - {} - - Scratch (const Scratch &data) - : - fe_collection (data.fe_collection), - quadrature_collection (data.quadrature_collection), - x_fe_values (fe_collection, - quadrature_collection, - data.update_flags), - rhs_values (data.rhs_values), - update_flags (data.update_flags) - {} - - const ::dealii::hp::FECollection &fe_collection; - const ::dealii::hp::QCollection &quadrature_collection; - - ::dealii::hp::FEValues x_fe_values; - - std::vector rhs_values; - - const UpdateFlags update_flags; - }; - - - struct CopyData - { - std::vector dof_indices; - FullMatrix cell_matrix; - dealii::Vector cell_rhs; - const ConstraintMatrix *constraints; - }; -} - -template - void - mass_assembler(const CellIterator &cell, Scratch &data, - CopyData ©_data) - { - data.x_fe_values.reinit(cell); - const FEValues &fe_values = - data.x_fe_values.get_present_fe_values(); - - const unsigned int dofs_per_cell = fe_values.dofs_per_cell, n_q_points = - fe_values.n_quadrature_points; - const FiniteElement &fe = fe_values.get_fe(); - const unsigned int n_components = fe.n_components(); + virtual void vector_value (const Point &p, + Vector &v) const { - copy_data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell); - copy_data.cell_rhs.reinit(dofs_per_cell); - - copy_data.dof_indices.resize(dofs_per_cell); + for (unsigned int c=0; cget_dof_indices(copy_data.dof_indices); - - data.rhs_values.resize(n_q_points); - F().value_list(fe_values.get_quadrature_points(), - data.rhs_values); - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - const unsigned int component_i = fe.system_to_component_index(i).first; - const double *phi_i = &fe_values.shape_value(i, 0); - - double add_data = 0; - for (unsigned int point = 0; point < n_q_points; ++point) - add_data += phi_i[point] * fe_values.JxW(point) - * data.rhs_values[point]; - copy_data.cell_rhs(i) = add_data; - } - } + private: + const unsigned int q; +}; -template -void copy_local_to_global (const CopyData &data, - MatrixType *matrix, - VectorType *right_hand_side) +template +void do_project (const Triangulation &triangulation, + const FiniteElement &fe, + const unsigned int p, + const unsigned int order_difference) { - const unsigned int dofs_per_cell = data.dof_indices.size(); - data.constraints->distribute_local_to_global(data.cell_matrix, - data.cell_rhs, - data.dof_indices, - *matrix, *right_hand_side); -} - - + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); -template -void create_mass_matrix (const DoFHandler &dof, - const Quadrature &q, - SparseMatrix &matrix, - Vector &rhs_vector) -{ - Assert (matrix.m() == dof.n_dofs(), - ExcDimensionMismatch (matrix.m(), dof.n_dofs())); - Assert (matrix.n() == dof.n_dofs(), - ExcDimensionMismatch (matrix.n(), dof.n_dofs())); + std::cout << "n_dofs=" << dof_handler.n_dofs() << std::endl; - hp::FECollection fe_collection (dof.get_fe()); - hp::QCollection q_collection (q); - Scratch - assembler_data (fe_collection, - update_values | - update_JxW_values | update_quadrature_points, - q_collection); - CopyData copy_data; - copy_data.cell_matrix.reinit (assembler_data.fe_collection.max_dofs_per_cell(), - assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.cell_rhs.reinit (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.dof_indices.resize (assembler_data.fe_collection.max_dofs_per_cell()); - copy_data.constraints = &constraints; + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + constraints.close (); - dealii::WorkStream::run (dof.begin_active(), - static_cast::active_cell_iterator>(dof.end()), - &mass_assembler::active_cell_iterator>, - std_cxx1x::bind(& - copy_local_to_global, Vector >, - std_cxx1x::_1, &matrix, &rhs_vector), - assembler_data, - copy_data, - 2*multithread_info.n_default_threads, - 1); + Vector projection (dof_handler.n_dofs()); + Vector error (triangulation.n_active_cells()); + for (unsigned int q=0; q<=p+2-order_difference; ++q) + { + // project the function + VectorTools::project (dof_handler, + constraints, + QGauss(p+2), + F (q, fe.n_components()), + projection); + // just to make sure it doesn't get + // forgotten: handle hanging node + // constraints + constraints.distribute (projection); + + // then compute the interpolation error + VectorTools::integrate_difference (dof_handler, + projection, + F (q, fe.n_components()), + error, + QGauss(std::max(p,q)+1), + VectorTools::L2_norm); + std::cout << fe.get_name() << ", P_" << q + << ", rel. error=" << error.l2_norm() / projection.l2_norm() + << std::endl; + + if (q<=p-order_difference) + if (error.l2_norm() > 1e-10*projection.l2_norm()) + { + std::cout << "Projection failed with relative error " + << error.l2_norm() / projection.l2_norm() + << std::endl; + Assert (false, ExcInternalError()); + } + } } + +// check the given element of polynomial order p. the last parameter, if +// given, denotes a gap in convergence order; for example, the Nedelec element +// of polynomial degree p has normal components of degree p-1 and therefore +// can only represent polynomials of degree p-1 exactly. the gap is then 1. template -void do_project (const unsigned int p) +void test_no_hanging_nodes (const FiniteElement &fe, + const unsigned int p, + const unsigned int order_difference = 0) { - SparsityPattern sparsity; - { - CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(), dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, - false); - - sparsity.copy_from (csp); - } - - SparseMatrix mass_matrix (sparsity); - Vector tmp (mass_matrix.n()); - - const Function<2>* dummy = 0; - create_mass_matrix(dof_handler, QGauss < dim > (5), mass_matrix, tmp); + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); - double l1 = 0; - for (unsigned int i=0; i void test () { - Threads::TaskGroup<> g; - for (unsigned int p=1; p<48*3; ++p) - g += Threads::new_task (&do_project, 3); - g.join_all (); + test_no_hanging_nodes (FE_Q(1), 1); } @@ -297,13 +183,5 @@ int main () deallog.depth_console(0); deallog.threshold_double(1.e-10); - GridGenerator::hyper_cube (triangulation, 0, 1); - triangulation.refine_global (1); - dof_handler.distribute_dofs (fe); - constraints.close (); - - test<2>(); - - deallog << "OK" << std::endl; + test<3>(); } -