From 676386f2dc5f69e98f1b5a401bd8a310b50f96f4 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Thu, 9 May 2013 12:55:58 +0000 Subject: [PATCH] Extend MatrixTools::create_mass_matrix and MatrixTools::create_laplace_matrix. git-svn-id: https://svn.dealii.org/trunk@29486 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 + .../include/deal.II/numerics/matrix_tools.h | 295 +++++---- .../deal.II/numerics/vector_tools.templates.h | 565 ++++++------------ deal.II/source/numerics/matrix_tools.cc | 92 ++- deal.II/source/numerics/matrix_tools.inst.in | 72 ++- .../create_laplace_matrix_constraints_01.cc | 139 +++++ .../cmp/generic | 7 + .../create_laplace_matrix_constraints_02.cc | 143 +++++ .../cmp/generic | 7 + .../create_laplace_matrix_constraints_02b.cc | 136 +++++ .../cmp/generic | 4 + .../create_laplace_matrix_constraints_03.cc | 140 +++++ .../cmp/generic | 7 + .../create_laplace_matrix_constraints_04.cc | 144 +++++ .../cmp/generic | 7 + .../create_mass_matrix_constraints_01.cc | 139 +++++ .../cmp/generic | 7 + .../create_mass_matrix_constraints_02.cc | 143 +++++ .../cmp/generic | 7 + .../create_mass_matrix_constraints_02b.cc | 136 +++++ .../cmp/generic | 4 + .../create_mass_matrix_constraints_03.cc | 140 +++++ .../cmp/generic | 7 + .../create_mass_matrix_constraints_04.cc | 144 +++++ .../cmp/generic | 7 + .../create_laplace_matrix_constraints_02.cc | 157 +++++ .../cmp/generic | 7 + .../create_laplace_matrix_constraints_02b.cc | 149 +++++ .../cmp/generic | 4 + tests/hp/create_mass_matrix_constraints_02.cc | 157 +++++ .../cmp/generic | 7 + .../hp/create_mass_matrix_constraints_02b.cc | 149 +++++ .../cmp/generic | 4 + 33 files changed, 2551 insertions(+), 584 deletions(-) create mode 100644 tests/deal.II/create_laplace_matrix_constraints_01.cc create mode 100644 tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic create mode 100644 tests/deal.II/create_laplace_matrix_constraints_02.cc create mode 100644 tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic create mode 100644 tests/deal.II/create_laplace_matrix_constraints_02b.cc create mode 100644 tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic create mode 100644 tests/deal.II/create_laplace_matrix_constraints_03.cc create mode 100644 tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic create mode 100644 tests/deal.II/create_laplace_matrix_constraints_04.cc create mode 100644 tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic create mode 100644 tests/deal.II/create_mass_matrix_constraints_01.cc create mode 100644 tests/deal.II/create_mass_matrix_constraints_01/cmp/generic create mode 100644 tests/deal.II/create_mass_matrix_constraints_02.cc create mode 100644 tests/deal.II/create_mass_matrix_constraints_02/cmp/generic create mode 100644 tests/deal.II/create_mass_matrix_constraints_02b.cc create mode 100644 tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic create mode 100644 tests/deal.II/create_mass_matrix_constraints_03.cc create mode 100644 tests/deal.II/create_mass_matrix_constraints_03/cmp/generic create mode 100644 tests/deal.II/create_mass_matrix_constraints_04.cc create mode 100644 tests/deal.II/create_mass_matrix_constraints_04/cmp/generic create mode 100644 tests/hp/create_laplace_matrix_constraints_02.cc create mode 100644 tests/hp/create_laplace_matrix_constraints_02/cmp/generic create mode 100644 tests/hp/create_laplace_matrix_constraints_02b.cc create mode 100644 tests/hp/create_laplace_matrix_constraints_02b/cmp/generic create mode 100644 tests/hp/create_mass_matrix_constraints_02.cc create mode 100644 tests/hp/create_mass_matrix_constraints_02/cmp/generic create mode 100644 tests/hp/create_mass_matrix_constraints_02b.cc create mode 100644 tests/hp/create_mass_matrix_constraints_02b/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ace3025a09..26c06d08d1 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -117,6 +117,15 @@ this function.
    +
  1. Improved: The functions MatrixTools::create_mass_matrix and +MatrixTools::create_laplace_matrix take now an optional ConstraintMatrix +argument that allows to directly apply the constraints. This also helps +VectorTools::project. Note that not providing constraints remains the default +and recommended way to ensure consistency when several matrices are added. +
    +(Martin Kronbichler, 2013/05/08) +
    2. +
  2. New: The classes TrilinosWrappers::SparseMatrix and TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with deal.II's own vector classes Vector and diff --git a/deal.II/include/deal.II/numerics/matrix_tools.h b/deal.II/include/deal.II/numerics/matrix_tools.h index 38be9336e2..abb91c8194 100644 --- a/deal.II/include/deal.II/numerics/matrix_tools.h +++ b/deal.II/include/deal.II/numerics/matrix_tools.h @@ -26,6 +26,8 @@ DEAL_II_NAMESPACE_OPEN // forward declarations template class Quadrature; +class ConstraintMatrix; + template class Vector; template class FullMatrix; template class SparseMatrix; @@ -105,17 +107,27 @@ namespace TrilinosWrappers * the matrix wasn't empty before. Therefore you may want to clear the matrix * before assemblage. * - * All created matrices are `raw': they are not condensed, + * By default, all created matrices are `raw': they are not condensed, * i.e. hanging nodes are not eliminated. The reason is that you may * want to add several matrices and could then condense afterwards * only once, instead of for every matrix. To actually do computations * with these matrices, you have to condense the matrix using the * ConstraintMatrix::condense function; you also have to * condense the right hand side accordingly and distribute the - * solution afterwards. + * solution afterwards. Alternatively, you can give an optional argument + * ConstraintMatrix that writes cell matrix (and vector) entries with + * distribute_local_to_global into the global matrix and vector. This way, + * adding several matrices from different sources is more complicated and + * you should make sure that you do not mix different ways of applying + * constraints. Particular caution is necessary when the given + * constraint matrix contains inhomogeneous constraints: In that case, the + * matrix assembled this way must be the only matrix (or you need to + * assemble the same right hand side for every matrix + * you generate and add together). * * If you want to use boundary conditions with the matrices generated - * by the functions of this class, you have to use a function like + * by the functions of this class in addition to the ones in a possible + * constraint matrix, you have to use a function like * ProblemBase<>::apply_dirichlet_bc to matrix and right hand * side. * @@ -210,154 +222,141 @@ namespace TrilinosWrappers namespace MatrixCreator { /** - * Assemble the mass matrix. If no - * coefficient is given, it is assumed - * to be unity. + * Assemble the mass matrix. If no coefficient is given, it is assumed to be + * unity. + * + * If the library is configured to use multithreading, this function works + * in parallel. * - * If the library is configured - * to use multithreading, this - * function works in parallel. + * The optional argument @p constraints allows to apply constraints on the + * resulting matrix directly. Be careful when combining several matrices and + * using inhomogeneous constraints. * - * See the general doc of this class - * for more information. + * See the general doc of this class for more information. */ template void create_mass_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Calls the create_mass_matrix() - * function, see above, with + * Calls the create_mass_matrix() function, see above, with * mapping=MappingQ1@(). */ template void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Assemble the mass matrix and a - * right hand side vector. If no - * coefficient is given, it is - * assumed to be unity. + * Assemble the mass matrix and a right hand side vector. If no coefficient + * is given, it is assumed to be unity. + * + * If the library is configured to use multithreading, this function works + * in parallel. * - * If the library is configured - * to use multithreading, this - * function works in parallel. + * The optional argument @p constraints allows to apply constraints on the + * resulting matrix directly. Be careful when combining several matrices and + * using inhomogeneous constraints. * - * See the general doc of this - * class for more information. + * See the general doc of this class for more information. */ template - void create_mass_matrix (const Mapping &mapping, - const DoFHandler &dof, + void create_mass_matrix (const Mapping &mapping, + const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Calls the create_mass_matrix() - * function, see above, with + * Calls the create_mass_matrix() function, see above, with * mapping=MappingQ1@(). */ template - void create_mass_matrix (const DoFHandler &dof, + void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); + /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template void create_mass_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template void create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template - void create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, - const hp::QCollection &q, + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template - void create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Assemble the mass matrix and a - * right hand side vector along - * the boundary. + * Assemble the mass matrix and a right hand side vector along the boundary. * - * The matrix is assumed to - * already be initialized with a - * suiting sparsity pattern (the - * DoFHandler provides an - * appropriate function). + * The matrix is assumed to already be initialized with a suiting sparsity + * pattern (the DoFHandler provides an appropriate function). * - * If the library is configured - * to use multithreading, this - * function works in parallel. + * If the library is configured to use multithreading, this function works + * in parallel. * - * @arg @p weight: an optional - * weight for the computation of - * the mass matrix. If no weight - * is given, it is set to one. + * @arg @p weight: an optional weight for the computation of the mass + * matrix. If no weight is given, it is set to one. * - * @arg @p component_mapping: if - * the components in @p - * boundary_functions and @p dof - * do not coincide, this vector - * allows them to be remapped. If - * the vector is not empty, it - * has to have one entry for each - * component in @p dof. This - * entry is the component number - * in @p boundary_functions that - * should be used for this - * component in @p dof. By - * default, no remapping is + * @arg @p component_mapping: if the components in @p boundary_functions and + * @p dof do not coincide, this vector allows them to be remapped. If the + * vector is not empty, it has to have one entry for each component in @p + * dof. This entry is the component number in @p boundary_functions that + * should be used for this component in @p dof. By default, no remapping is * applied. * - * @todo This function does not - * work for finite elements with - * cell-dependent shape functions. + * @todo This function does not work for finite elements with cell-dependent + * shape functions. */ template @@ -396,9 +395,7 @@ namespace MatrixCreator /** - * Calls the - * create_boundary_mass_matrix() - * function, see above, with + * Calls the create_boundary_mass_matrix() function, see above, with * mapping=MappingQ1@(). */ template @@ -413,8 +410,7 @@ namespace MatrixCreator std::vector component_mapping = std::vector()); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template @@ -429,8 +425,7 @@ namespace MatrixCreator std::vector component_mapping = std::vector()); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ // // void create_boundary_mass_matrix (const hp::MappingCollection<1,1> &mapping, @@ -443,8 +438,7 @@ namespace MatrixCreator // const Function<1> * const a = 0); /** - * Same function as above, but for hp - * objects. + * Same function as above, but for hp objects. */ template @@ -458,123 +452,122 @@ namespace MatrixCreator std::vector component_mapping = std::vector()); /** - * Assemble the Laplace - * matrix. If no coefficient is - * given, it is assumed to be - * constant one. + * Assemble the Laplace matrix. If no coefficient is given, it is assumed to + * be constant one. + * + * If the library is configured to use multithreading, this function works + * in parallel. * - * If the library is configured - * to use multithreading, this - * function works in parallel. + * The optional argument @p constraints allows to apply constraints on the + * resulting matrix directly. Be careful when combining several matrices and + * using inhomogeneous constraints. * - * See the general doc of this - * class for more information. + * See the general doc of this class for more information. */ template - void create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Calls the - * create_laplace_matrix() - * function, see above, with + * Calls the create_laplace_matrix() function, see above, with * mapping=MappingQ1@(). */ template - void create_laplace_matrix (const DoFHandler &dof, + void create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Assemble the Laplace matrix - * and a right hand side - * vector. If no coefficient is - * given, it is assumed to be - * constant one. + * Assemble the Laplace matrix and a right hand side vector. If no + * coefficient is given, it is assumed to be constant one. + * + * If the library is configured to use multithreading, this function works + * in parallel. * - * If the library is configured - * to use multithreading, this - * function works in parallel. + * The optional argument @p constraints allows to apply constraints on the + * resulting matrix directly. Be careful when combining several matrices and + * using inhomogeneous constraints. * - * See the general doc of this - * class for more information. + * See the general doc of this class for more information. */ template - void create_laplace_matrix (const Mapping &mapping, - const DoFHandler &dof, + void create_laplace_matrix (const Mapping &mapping, + const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Calls the - * create_laplace_matrix() - * function, see above, with + * Calls the create_laplace_matrix() function, see above, with * mapping=MappingQ1@(). */ template - void create_laplace_matrix (const DoFHandler &dof, + void create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Like the functions above, but for hp - * dof handlers, mappings, and quadrature - * collections. + * Like the functions above, but for hp dof handlers, mappings, and + * quadrature collections. */ template - void create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Like the functions above, but for hp - * dof handlers, mappings, and quadrature - * collections. + * Like the functions above, but for hp dof handlers, mappings, and + * quadrature collections. */ template - void create_laplace_matrix (const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Like the functions above, but for hp - * dof handlers, mappings, and quadrature - * collections. + * Like the functions above, but for hp dof handlers, mappings, and + * quadrature collections. */ template - void create_laplace_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** - * Like the functions above, but for hp - * dof handlers, mappings, and quadrature - * collections. + * Like the functions above, but for hp dof handlers, mappings, and + * quadrature collections. */ template - void create_laplace_matrix (const hp::DoFHandler &dof, + void create_laplace_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const a = 0); + const Function *const a = 0, + const ConstraintMatrix &constraints = ConstraintMatrix()); /** * Exception diff --git a/deal.II/include/deal.II/numerics/vector_tools.templates.h b/deal.II/include/deal.II/numerics/vector_tools.templates.h index 129cc5f779..3dec6628b3 100644 --- a/deal.II/include/deal.II/numerics/vector_tools.templates.h +++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h @@ -557,11 +557,11 @@ namespace VectorTools template - void project (const Mapping &mapping, - const DoFHandler &dof, + void project (const Mapping &mapping, + const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VECTOR &vec_result, const bool enforce_zero_boundary, const Quadrature &q_boundary, @@ -607,11 +607,27 @@ namespace VectorTools Vector vec (dof.n_dofs()); SparsityPattern sparsity; + // check if constraints are compatible (see below) + bool constraints_are_compatible = true; + { + for (std::map::iterator it=boundary_values.begin(); + it != boundary_values.end(); ++it) + if (constraints.is_constrained(it->first)) + if (!(constraints.get_constraint_entries(it->first)->size() > 0 + || + (constraints.get_inhomogeneity(it->first) == it->second))) + { + constraints_are_compatible = false; + break; + } + } + // use csp to consume less memory and to // still be fast { CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs()); - DoFTools::make_sparsity_pattern (dof, csp, constraints); + DoFTools::make_sparsity_pattern (dof, csp, constraints, + !constraints_are_compatible); sparsity.copy_from (csp); } @@ -619,52 +635,60 @@ namespace VectorTools SparseMatrix mass_matrix (sparsity); Vector tmp (mass_matrix.n()); - // create mass matrix and rhs at once, - // which is faster. - MatrixCreator::create_mass_matrix (mapping, dof, quadrature, mass_matrix, - function, tmp); - - constraints.condense (mass_matrix); - constraints.condense (tmp); - if (boundary_values.size() != 0) - MatrixTools::apply_boundary_values (boundary_values, - mass_matrix, vec, tmp, - true); - // Allow for a maximum of 5*n - // steps to reduce the residual by - // 10^-12. n steps may not be - // sufficient, since roundoff - // errors may accumulate for badly - // conditioned matrices - ReductionControl control(5*tmp.size(), 0., 1e-12, false, false); + // If the constraint matrix does not conflict with the given boundary + // values (i.e., it either does not contain boundary values or it contains + // the same as boundary_values), we can let it call + // distribute_local_to_global straight away, otherwise we need to first + // interpolate the boundary values and then condense the matrix and vector + if (constraints_are_compatible) + { + const Function* dummy = 0; + MatrixCreator::create_mass_matrix (mapping, dof, quadrature, + mass_matrix, function, tmp, + dummy, constraints); + if (boundary_values.size() > 0) + MatrixTools::apply_boundary_values (boundary_values, + mass_matrix, vec, tmp, + true); + } + else + { + // create mass matrix and rhs at once, which is faster. + MatrixCreator::create_mass_matrix (mapping, dof, quadrature, + mass_matrix, function, tmp); + MatrixTools::apply_boundary_values (boundary_values, + mass_matrix, vec, tmp, + true); + constraints.condense(mass_matrix, tmp); + } + + // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n + // steps may not be sufficient, since roundoff errors may accumulate for + // badly conditioned matrices + ReductionControl control(5*tmp.size(), 0., 1e-12, false, false); GrowingVectorMemory<> memory; - SolverCG<> cg(control,memory); + SolverCG<> cg(control,memory); PreconditionSSOR<> prec; prec.initialize(mass_matrix, 1.2); - // solve + cg.solve (mass_matrix, vec, tmp, prec); - // distribute solution constraints.distribute (vec); - // copy vec into vec_result. we - // can't use ve_result itself - // above, since it may be of - // another type than Vector - // and that wouldn't necessarily go - // together with the matrix and - // other functions + // copy vec into vec_result. we can't use vec_result itself above, since + // it may be of another type than Vector and that wouldn't + // necessarily go together with the matrix and other functions for (unsigned int i=0; i - void project (const DoFHandler &dof, + void project (const DoFHandler &dof, const ConstraintMatrix &constraints, const Quadrature &quadrature, - const Function &function, + const Function &function, VECTOR &vec, const bool enforce_zero_boundary, const Quadrature &q_boundary, @@ -3063,25 +3087,12 @@ namespace VectorTools { 0, 1 } }; - // The projection is - // divided into two steps. - // In the first step we - // project the boundary - // function on the - // horizontal shape - // functions. Then the - // bounary function is - // projected on the - // vertical shape - // functions. We begin - // with the horizontal - // shape functions and - // set up a linear system - // of equations to get - // the values for degrees - // of freedom associated - // with the interior of - // the face. + // The projection is divided into two steps. In the first step we + // project the boundary function on the horizontal shape functions. + // Then the boundary function is projected on the vertical shape + // functions. We begin with the horizontal shape functions and set + // up a linear system of equations to get the values for degrees of + // freedom associated with the interior of the face. for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point) { @@ -4269,19 +4280,16 @@ namespace VectorTools std::vector face_dofs; - // create FE and mapping - // collections for all elements in - // use by this DoFHandler + // create FE and mapping collections for all elements in use by this + // DoFHandler hp::FECollection fe_collection (dof_handler.get_fe()); hp::MappingCollection mapping_collection; for (unsigned int i=0; i face_quadrature_collection; for (unsigned int i=0; i (unit_support_points)); } - // now create the object with which - // we will generate the normal - // vectors + // now create the object with which we will generate the normal vectors hp::FEFaceValues x_fe_face_values (mapping_collection, fe_collection, face_quadrature_collection, update_q_points); - // have a map that stores normal - // vectors for each vector-dof - // tuple we want to - // constrain. since we can get at - // the same vector dof tuple more - // than once (for example if it is - // located at a vertex that we - // visit from all adjacent cells), - // we will want to average later on - // the normal vectors computed on - // different cells as described in - // the documentation of this - // function. however, we can only - // average if the contributions - // came from different cells, - // whereas we want to constrain - // twice or more in case the - // contributions came from - // different faces of the same cell - // (i.e. constrain not just the - // *average normal direction* but - // *all normal directions* we - // find). consequently, we also - // have to store which cell a + // have a map that stores normal vectors for each vector-dof tuple we want + // to constrain. since we can get at the same vector dof tuple more than + // once (for example if it is located at a vertex that we visit from all + // adjacent cells), we will want to average later on the normal vectors + // computed on different cells as described in the documentation of this + // function. however, we can only average if the contributions came from + // different cells, whereas we want to constrain twice or more in case the + // contributions came from different faces of the same cell + // (i.e. constrain not just the *average normal direction* but *all normal + // directions* we find). consequently, we also have to store which cell a // normal vector was computed on typedef std::multimap, @@ -4349,24 +4341,20 @@ namespace VectorTools const FiniteElement &fe = cell->get_fe (); typename DH::face_iterator face = cell->face(face_no); - // get the indices of the - // dofs on this cell... + // get the indices of the dofs on this cell... face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); x_fe_face_values.reinit (cell, face_no); const FEFaceValues &fe_values = x_fe_face_values.get_present_fe_values(); - // then identify which of - // them correspond to the - // selected set of vector - // components + // then identify which of them correspond to the selected set of + // vector components for (unsigned int i=0; i vector_dofs; vector_dofs.dof_indices[0] = face_dofs[i]; @@ -4389,59 +4377,33 @@ namespace VectorTools Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(), ExcInternalError()); - // we need the normal - // vector on this face. we - // know that it is a vector - // of length 1 but at least - // with higher order - // mappings it isn't always - // possible to guarantee - // that each component is - // exact up to zero - // tolerance. in - // particular, as shown in - // the deal.II/no_flux_06 - // test, if we just take - // the normal vector as - // given by the fe_values - // object, we can get - // entries in the normal - // vectors of the unit cube - // that have entries up to + // we need the normal vector on this face. we know that it + // is a vector of length 1 but at least with higher order + // mappings it isn't always possible to guarantee that + // each component is exact up to zero tolerance. in + // particular, as shown in the deal.II/no_flux_06 test, if + // we just take the normal vector as given by the + // fe_values object, we can get entries in the normal + // vectors of the unit cube that have entries up to // several times 1e-14. // - // the problem with this is - // that this later yields - // constraints that are - // circular (e.g., in the - // testcase, we get - // constraints of the form + // the problem with this is that this later yields + // constraints that are circular (e.g., in the testcase, + // we get constraints of the form // // x22 = 2.93099e-14*x21 + 2.93099e-14*x23 // x21 = -2.93099e-14*x22 + 2.93099e-14*x21 // - // in both of these - // constraints, the small - // numbers should be zero - // and the constraints - // should simply be - // x22 = x21 = 0 + // in both of these constraints, the small numbers should + // be zero and the constraints should simply be + // x22 = x21 = 0 // - // to achieve this, we - // utilize that we know - // that the normal vector - // has (or should have) - // length 1 and that we can - // simply set small - // elements to zero - // (without having to check - // that they are small - // *relative to something - // else*). we do this and - // then normalize the - // length of the vector - // back to one, just to be - // on the safe side + // to achieve this, we utilize that we know that the + // normal vector has (or should have) length 1 and that we + // can simply set small elements to zero (without having + // to check that they are small *relative to something + // else*). we do this and then normalize the length of the + // vector back to one, just to be on the safe side Point normal_vector = (cell->face(face_no)->get_boundary() .normal_vector (cell->face(face_no), @@ -4453,9 +4415,8 @@ namespace VectorTools normal_vector[d] = 0; normal_vector /= normal_vector.norm(); - // now enter the - // (dofs,(normal_vector,cell)) - // entry into the map + // now enter the (dofs,(normal_vector,cell)) entry into + // the map dof_to_normals_map .insert (std::make_pair (vector_dofs, std::make_pair (normal_vector, @@ -4463,27 +4424,19 @@ namespace VectorTools } } - // Now do something with the - // collected information. To this - // end, loop through all sets of - // pairs (dofs,normal_vector) and - // identify which entries belong to - // the same set of dofs and then do - // as described in the - // documentation, i.e. either - // average the normal vector or - // don't for this particular set of - // dofs + // Now do something with the collected information. To this end, loop + // through all sets of pairs (dofs,normal_vector) and identify which + // entries belong to the same set of dofs and then do as described in the + // documentation, i.e. either average the normal vector or don't for this + // particular set of dofs typename DoFToNormalsMap::const_iterator p = dof_to_normals_map.begin(); while (p != dof_to_normals_map.end()) { - // first find the range of entries in - // the multimap that corresponds to the - // same vector-dof tuple. as usual, we - // define the range half-open. the - // first entry of course is 'p' + // first find the range of entries in the multimap that corresponds to + // the same vector-dof tuple. as usual, we define the range + // half-open. the first entry of course is 'p' typename DoFToNormalsMap::const_iterator same_dof_range[2] = { p }; for (++p; p != dof_to_normals_map.end(); ++p) @@ -4495,13 +4448,10 @@ namespace VectorTools if (p == dof_to_normals_map.end()) same_dof_range[1] = dof_to_normals_map.end(); - // now compute the reverse mapping: for - // each of the cells that contributed - // to the current set of vector dofs, - // add up the normal vectors. the - // values of the map are pairs of - // normal vectors and number of cells - // that have contributed + // now compute the reverse mapping: for each of the cells that + // contributed to the current set of vector dofs, add up the normal + // vectors. the values of the map are pairs of normal vectors and + // number of cells that have contributed typedef std::map ::active_cell_iterator, @@ -4525,9 +4475,7 @@ namespace VectorTools Assert (old_count > 0, ExcInternalError()); - // in the same entry, - // store again the now - // averaged normal vector + // in the same entry, store again the now averaged normal vector // and the new count cell_to_normals_map[q->second.second] = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1), @@ -4536,8 +4484,7 @@ namespace VectorTools Assert (cell_to_normals_map.size() >= 1, ExcInternalError()); - // count the maximum number of - // contributions from each cell + // count the maximum number of contributions from each cell unsigned int max_n_contributions_per_cell = 1; for (typename CellToNormalsMap::const_iterator x = cell_to_normals_map.begin(); @@ -4546,40 +4493,26 @@ namespace VectorTools = std::max (max_n_contributions_per_cell, x->second.second); - // verify that each cell can have only - // contributed at most dim times, since - // that is the maximum number of faces - // that come together at a single place + // verify that each cell can have only contributed at most dim times, + // since that is the maximum number of faces that come together at a + // single place Assert (max_n_contributions_per_cell <= dim, ExcInternalError()); switch (max_n_contributions_per_cell) { - // first deal with the case that a - // number of cells all have - // registered that they have a - // normal vector defined at the - // location of a given vector dof, - // and that each of them have - // encountered this vector dof - // exactly once while looping over - // all their faces. as stated in - // the documentation, this is the - // case where we want to simply - // average over all normal vectors + // first deal with the case that a number of cells all have + // registered that they have a normal vector defined at the + // location of a given vector dof, and that each of them have + // encountered this vector dof exactly once while looping over all + // their faces. as stated in the documentation, this is the case + // where we want to simply average over all normal vectors case 1: { - // compute the average - // normal vector from all - // the ones that have the - // same set of dofs. we - // could add them up and - // divide them by the - // number of additions, - // or simply normalize - // them right away since - // we want them to have - // unit length anyway + // compute the average normal vector from all the ones that have + // the same set of dofs. we could add them up and divide them by + // the number of additions, or simply normalize them right away + // since we want them to have unit length anyway Tensor<1,dim> normal; for (typename CellToNormalsMap::const_iterator x = cell_to_normals_map.begin(); @@ -4593,8 +4526,7 @@ namespace VectorTools normal[d] = 0; normal /= normal.norm(); - // then construct constraints - // from this: + // then construct constraints from this: const internal::VectorDoFTuple & dof_indices = same_dof_range[0]->first; internal::add_constraint (dof_indices, normal, @@ -4604,56 +4536,27 @@ namespace VectorTools } - // this is the slightly - // more complicated case - // that a single cell has - // contributed with exactly - // DIM normal vectors to - // the same set of vector - // dofs. this is what - // happens in a corner in - // 2d and 3d (but not on an - // edge in 3d, where we - // have only 2, i.e. t; @@ -4666,14 +4569,10 @@ namespace VectorTools ExcMessage("Found a set of normal vectors that are nearly collinear.")); } - // so all components of - // this vector dof are - // constrained. enter - // this into the - // constraint matrix + // so all components of this vector dof are constrained. enter + // this into the constraint matrix // - // ignore dofs already - // constrained + // ignore dofs already constrained for (unsigned int i=0; ifirst.dof_indices[i]) @@ -4689,31 +4588,18 @@ namespace VectorTools } - // this is the case of an - // edge contribution in 3d, - // i.e. the vector is - // constrained in two - // directions but not the - // third. + // this is the case of an edge contribution in 3d, i.e. the vector + // is constrained in two directions but not the third. default: { Assert (dim >= 3, ExcNotImplemented()); Assert (max_n_contributions_per_cell == 2, ExcInternalError()); - // as described in the - // documentation, let us - // first collect what - // each of the cells - // contributed at the - // current point. we use - // a std::list instead of - // a std::set (which - // would be more natural) - // because std::set - // requires that the - // stored elements are - // comparable with - // operator< + // as described in the documentation, let us first collect what + // each of the cells contributed at the current point. we use a + // std::list instead of a std::set (which would be more natural) + // because std::set requires that the stored elements are + // comparable with operator< typedef std::map::active_cell_iterator, std::list > > CellContributions; @@ -4725,24 +4611,13 @@ namespace VectorTools cell_contributions[q->second.second].push_back (q->second.first); Assert (cell_contributions.size() >= 1, ExcInternalError()); - // now for each cell that - // has contributed - // determine the number - // of normal vectors it - // has contributed. we - // currently only - // implement if this is - // dim-1 for all cells - // (if a single cell has - // contributed dim, or if - // all adjacent cells - // have contributed 1 - // normal vector, this is - // already handled above) + // now for each cell that has contributed determine the number of + // normal vectors it has contributed. we currently only implement + // if this is dim-1 for all cells (if a single cell has + // contributed dim, or if all adjacent cells have contributed 1 + // normal vector, this is already handled above) // - // for each contributing - // cell compute the - // tangential vector that + // for each contributing cell compute the tangential vector that // remains unconstrained std::list > tangential_vectors; for (typename CellContributions::const_iterator @@ -4763,81 +4638,24 @@ namespace VectorTools Assert (index == dim-1, ExcInternalError()); } - // calculate the - // tangent as the - // outer product of - // the normal - // vectors. since - // these vectors do - // not need to be - // orthogonal (think, - // for example, the - // case of the - // deal.II/no_flux_07 - // test: a sheared - // cube in 3d, with - // Q2 elements, where - // we have - // constraints from - // the two normal - // vectors of two - // faces of the - // sheared cube that - // are not - // perpendicular to - // each other), we - // have to normalize - // the outer product + // calculate the tangent as the outer product of the normal + // vectors. since these vectors do not need to be orthogonal + // (think, for example, the case of the deal.II/no_flux_07 + // test: a sheared cube in 3d, with Q2 elements, where we have + // constraints from the two normal vectors of two faces of the + // sheared cube that are not perpendicular to each other), we + // have to normalize the outer product Tensor<1,dim> tangent; switch (dim) { case 3: - // take - // cross - // product - // between - // normals[0] - // and - // normals[1]. write - // it in - // the - // current - // form - // (with - // [dim-2]) - // to make - // sure - // that - // compilers - // don't - // warn - // about - // out-of-bounds - // accesses - // -- the - // warnings - // are - // bogus - // since we - // get here - // only for - // dim==3, - // but at - // least - // one - // isn't - // quite - // smart - // enough - // to - // notice - // this and - // warns - // when - // compiling - // the - // function - // in 2d + // take cross product between normals[0] and + // normals[1]. write it in the current form (with [dim-2]) + // to make sure that compilers don't warn about + // out-of-bounds accesses -- the warnings are bogus since + // we get here only for dim==3, but at least one isn't + // quite smart enough to notice this and warns when + // compiling the function in 2d cross_product (tangent, normals[0], normals[dim-2]); break; default: @@ -4852,15 +4670,9 @@ namespace VectorTools tangential_vectors.push_back (tangent); } - // go through the list of - // tangents and make sure - // that they all roughly - // point in the same - // direction as the first - // one (i.e. have an - // angle less than 90 - // degrees); if they - // don't then flip their + // go through the list of tangents and make sure that they all + // roughly point in the same direction as the first one (i.e. have + // an angle less than 90 degrees); if they don't then flip their // sign { const Tensor<1,dim> first_tangent = tangential_vectors.front(); @@ -4872,9 +4684,7 @@ namespace VectorTools *t *= -1; } - // now compute the - // average tangent and - // normalize it + // now compute the average tangent and normalize it Tensor<1,dim> average_tangent; for (typename std::list >::const_iterator t = tangential_vectors.begin(); @@ -4883,11 +4693,8 @@ namespace VectorTools average_tangent += *t; average_tangent /= average_tangent.norm(); - // now all that is left - // is that we add the - // constraints that the - // vector is parallel - // to the tangent + // now all that is left is that we add the constraints that the + // vector is parallel to the tangent const internal::VectorDoFTuple & dof_indices = same_dof_range[0]->first; internal::add_tangentiality_constraints (dof_indices, diff --git a/deal.II/source/numerics/matrix_tools.cc b/deal.II/source/numerics/matrix_tools.cc index 776031f0f5..77fa2af43d 100644 --- a/deal.II/source/numerics/matrix_tools.cc +++ b/deal.II/source/numerics/matrix_tools.cc @@ -209,6 +209,7 @@ namespace MatrixCreator std::vector dof_indices; FullMatrix cell_matrix; dealii::Vector cell_rhs; + const ConstraintMatrix *constraints; }; } @@ -592,10 +593,15 @@ namespace MatrixCreator (data.cell_rhs.size() == dofs_per_cell), ExcInternalError()); - matrix->add(data.dof_indices, data.cell_matrix); if (right_hand_side != 0) - for (unsigned int i=0; idistribute_local_to_global(data.cell_matrix, + data.cell_rhs, + data.dof_indices, + *matrix, *right_hand_side); + else + data.constraints->distribute_local_to_global(data.cell_matrix, + data.dof_indices, + *matrix); } } } @@ -609,7 +615,8 @@ namespace MatrixCreator const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -631,6 +638,7 @@ namespace MatrixCreator 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; WorkStream::run (dof.begin_active(), static_cast::active_cell_iterator>(dof.end()), @@ -648,10 +656,11 @@ namespace MatrixCreator void create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_mass_matrix(StaticMappingQ1::mapping, dof, - q, matrix, coefficient); + q, matrix, coefficient, constraints); } @@ -663,7 +672,8 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -684,6 +694,7 @@ namespace MatrixCreator 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; WorkStream::run (dof.begin_active(), static_cast::active_cell_iterator>(dof.end()), @@ -703,10 +714,12 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_mass_matrix(StaticMappingQ1::mapping, - dof, q, matrix, rhs, rhs_vector, coefficient); + dof, q, matrix, rhs, rhs_vector, coefficient, + constraints); } @@ -716,7 +729,8 @@ namespace MatrixCreator const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -734,6 +748,7 @@ namespace MatrixCreator 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; WorkStream::run (dof.begin_active(), static_cast::active_cell_iterator>(dof.end()), @@ -748,24 +763,27 @@ namespace MatrixCreator template - void create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { - create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); + create_mass_matrix(hp::StaticMappingQ1::mapping_collection, + dof, q, matrix, coefficient, constraints); } template - void create_mass_matrix (const hp::MappingCollection &mapping, - const hp::DoFHandler &dof, + void create_mass_matrix (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -783,6 +801,7 @@ namespace MatrixCreator 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; WorkStream::run (dof.begin_active(), static_cast::active_cell_iterator>(dof.end()), @@ -797,15 +816,16 @@ namespace MatrixCreator template - void create_mass_matrix (const hp::DoFHandler &dof, - const hp::QCollection &q, + void create_mass_matrix (const hp::DoFHandler &dof, + const hp::QCollection &q, SparseMatrix &matrix, - const Function &rhs, + const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_mass_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, - matrix, rhs, rhs_vector, coefficient); + matrix, rhs, rhs_vector, coefficient, constraints); } @@ -1706,7 +1726,8 @@ namespace MatrixCreator const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -1727,6 +1748,7 @@ namespace MatrixCreator 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; void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, SparseMatrix *, @@ -1749,7 +1771,8 @@ namespace MatrixCreator void create_laplace_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_laplace_matrix(StaticMappingQ1::mapping, dof, q, matrix, coefficient); } @@ -1763,7 +1786,8 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -1784,6 +1808,7 @@ namespace MatrixCreator 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; void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, SparseMatrix *, @@ -1808,7 +1833,8 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_laplace_matrix(StaticMappingQ1::mapping, dof, q, matrix, rhs, rhs_vector, coefficient); @@ -1821,7 +1847,8 @@ namespace MatrixCreator const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -1839,6 +1866,7 @@ namespace MatrixCreator 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; void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, SparseMatrix *, @@ -1861,7 +1889,8 @@ namespace MatrixCreator void create_laplace_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, coefficient); } @@ -1875,7 +1904,8 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { Assert (matrix.m() == dof.n_dofs(), ExcDimensionMismatch (matrix.m(), dof.n_dofs())); @@ -1893,6 +1923,7 @@ namespace MatrixCreator 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; void (*copy_local_to_global) (const MatrixCreator::internal::AssemblerData::CopyData &, SparseMatrix *, @@ -1917,7 +1948,8 @@ namespace MatrixCreator SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function *const coefficient) + const Function *const coefficient, + const ConstraintMatrix &constraints) { create_laplace_matrix(hp::StaticMappingQ1::mapping_collection, dof, q, matrix, rhs, rhs_vector, coefficient); diff --git a/deal.II/source/numerics/matrix_tools.inst.in b/deal.II/source/numerics/matrix_tools.inst.in index ad686e7118..07c6301155 100644 --- a/deal.II/source/numerics/matrix_tools.inst.in +++ b/deal.II/source/numerics/matrix_tools.inst.in @@ -21,13 +21,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const Mapping &mapping, @@ -36,7 +38,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, @@ -44,7 +47,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template @@ -53,13 +57,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const Mapping &mapping, @@ -68,7 +74,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const DoFHandler &dof, @@ -76,7 +83,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_boundary_mass_matrix @@ -121,13 +129,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::MappingCollection &mapping, @@ -136,7 +146,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, @@ -144,7 +155,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template @@ -153,13 +165,15 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::MappingCollection &mapping, @@ -168,7 +182,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_mass_matrix (const hp::DoFHandler &dof, @@ -176,7 +191,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); #if deal_II_dimension != 1 @@ -216,14 +232,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const Mapping &mapping, const DoFHandler &dof, const Quadrature &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const Mapping &mapping, @@ -232,7 +250,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const DoFHandler &dof, @@ -240,7 +259,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); // hp versions of create_laplace_matrix template @@ -248,14 +268,16 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS (const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const hp::MappingCollection &mapping, const hp::DoFHandler &dof, const hp::QCollection &q, SparseMatrix &matrix, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const hp::MappingCollection &mapping, @@ -264,7 +286,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); template void MatrixCreator::create_laplace_matrix (const hp::DoFHandler &dof, @@ -272,7 +295,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS SparseMatrix &matrix, const Function &rhs, Vector &rhs_vector, - const Function * const coefficient); + const Function * const coefficient, + const ConstraintMatrix &constraints); diff --git a/tests/deal.II/create_laplace_matrix_constraints_01.cc b/tests/deal.II/create_laplace_matrix_constraints_01.cc new file mode 100644 index 0000000000..3d1fca2ab9 --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_01.cc @@ -0,0 +1,139 @@ +//---------------------------- create_laplace_matrix_constraints_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_01.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FE_Q element(2); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1,1.), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_laplace_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_01/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_laplace_matrix_constraints_02.cc b/tests/deal.II/create_laplace_matrix_constraints_02.cc new file mode 100644 index 0000000000..a0e2e8299f --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_02.cc @@ -0,0 +1,143 @@ +//---------------------------- create_laplace_matrix_constraints_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_02.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_laplace_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_02/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_laplace_matrix_constraints_02b.cc b/tests/deal.II/create_laplace_matrix_constraints_02b.cc new file mode 100644 index 0000000000..8e868db6ad --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_02b.cc @@ -0,0 +1,136 @@ +//---------------------------- create_laplace_matrix_constraints_02b.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_02b.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints (and similar to +// create_laplace_matrix_constraints_02 but without rhs). + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ZeroFunction(2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + MatrixTools:: + create_laplace_matrix (mapping, dof, + quadrature, matrix_ref); + constraints.condense(matrix_ref); + + const Function *const dummy = 0; + MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix, + dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic new file mode 100644 index 0000000000..6a3d892caa --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_02b/cmp/generic @@ -0,0 +1,4 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:3d::Matrix error Frobenius: 0 diff --git a/tests/deal.II/create_laplace_matrix_constraints_03.cc b/tests/deal.II/create_laplace_matrix_constraints_03.cc new file mode 100644 index 0000000000..662185cad8 --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_03.cc @@ -0,0 +1,140 @@ +//---------------------------- create_laplace_matrix_constraints_03.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_03.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FE_Q element(2); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1,1.), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_laplace_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref, + &rhs_function); + constraints.condense(matrix_ref, rhs_ref); + + MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, &rhs_function, + constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_03/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_laplace_matrix_constraints_04.cc b/tests/deal.II/create_laplace_matrix_constraints_04.cc new file mode 100644 index 0000000000..5e0237afe2 --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_04.cc @@ -0,0 +1,144 @@ +//---------------------------- create_laplace_matrix_constraints_04.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_04.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_laplace_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref, + &rhs_function); + constraints.condense(matrix_ref, rhs_ref); + + MatrixTools::create_laplace_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, &rhs_function, + constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic b/tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_laplace_matrix_constraints_04/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_mass_matrix_constraints_01.cc b/tests/deal.II/create_mass_matrix_constraints_01.cc new file mode 100644 index 0000000000..8e5682e50d --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_01.cc @@ -0,0 +1,139 @@ +//---------------------------- create_mass_matrix_constraints_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_01.cc --------------------------- + + +// Same as create_mass_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FE_Q element(2); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1,1.), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_mass_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_mass_matrix_constraints_01/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_01/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_01/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_mass_matrix_constraints_02.cc b/tests/deal.II/create_mass_matrix_constraints_02.cc new file mode 100644 index 0000000000..14971a18d6 --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_02.cc @@ -0,0 +1,143 @@ +//---------------------------- create_mass_matrix_constraints_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_02.cc --------------------------- + + +// Same as create_mass_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_mass_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_mass_matrix_constraints_02/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_02/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_02/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_mass_matrix_constraints_02b.cc b/tests/deal.II/create_mass_matrix_constraints_02b.cc new file mode 100644 index 0000000000..ad7406775a --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_02b.cc @@ -0,0 +1,136 @@ +//---------------------------- create_mass_matrix_constraints_02b.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_02b.cc --------------------------- + + +// Same as create_mass_matrix_0[1234]b but passing an additional constraint +// matrix and comparing results without constraints (and similar to +// create_mass_matrix_constraints_02 but without rhs). + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ZeroFunction(2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + MatrixTools:: + create_mass_matrix (mapping, dof, + quadrature, matrix_ref); + constraints.condense(matrix_ref); + + const Function *const dummy = 0; + MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix, + dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic new file mode 100644 index 0000000000..6a3d892caa --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_02b/cmp/generic @@ -0,0 +1,4 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:3d::Matrix error Frobenius: 0 diff --git a/tests/deal.II/create_mass_matrix_constraints_03.cc b/tests/deal.II/create_mass_matrix_constraints_03.cc new file mode 100644 index 0000000000..0c9d3a0d71 --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_03.cc @@ -0,0 +1,140 @@ +//---------------------------- create_mass_matrix_constraints_03.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_03.cc --------------------------- + + +// Same as create_mass_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FE_Q element(2); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1,1.), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_mass_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref, + &rhs_function); + constraints.condense(matrix_ref, rhs_ref); + + MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, &rhs_function, + constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_mass_matrix_constraints_03/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_03/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_03/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/deal.II/create_mass_matrix_constraints_04.cc b/tests/deal.II/create_mass_matrix_constraints_04.cc new file mode 100644 index 0000000000..385bda3f7d --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_04.cc @@ -0,0 +1,144 @@ +//---------------------------- create_mass_matrix_constraints_04.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_04.cc --------------------------- + + +// Same as create_mass_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + // create a system element composed + // of one Q1 and one Q2 element + FESystem element(FE_Q(1), 1, + FE_Q(2), 1); + DoFHandler dof(tr); + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_mass_matrix (mapping, dof, + quadrature, matrix_ref, + rhs_function, rhs_ref, + &rhs_function); + constraints.condense(matrix_ref, rhs_ref); + + MatrixTools::create_mass_matrix (mapping, dof, quadrature, matrix, + rhs_function, rhs, &rhs_function, + constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/deal.II/create_mass_matrix_constraints_04/cmp/generic b/tests/deal.II/create_mass_matrix_constraints_04/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/deal.II/create_mass_matrix_constraints_04/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/hp/create_laplace_matrix_constraints_02.cc b/tests/hp/create_laplace_matrix_constraints_02.cc new file mode 100644 index 0000000000..5326b7d85e --- /dev/null +++ b/tests/hp/create_laplace_matrix_constraints_02.cc @@ -0,0 +1,157 @@ +//---------------------------- create_laplace_matrix_constraints_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_02.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + hp::FECollection element; + element.push_back (FESystem (FE_Q(1), 1, + FE_Q(2), 1)); + element.push_back (FESystem (FE_Q(2), 1, + FE_Q(3), 1)); + element.push_back (FESystem (FE_Q(3), 1, + FE_Q(4), 1)); + + hp::DoFHandler dof(tr); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(); + cell != dof.end(); ++cell) + cell->set_active_fe_index (rand() % element.size()); + + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_laplace_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_laplace_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/hp/create_laplace_matrix_constraints_02/cmp/generic b/tests/hp/create_laplace_matrix_constraints_02/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/hp/create_laplace_matrix_constraints_02/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/hp/create_laplace_matrix_constraints_02b.cc b/tests/hp/create_laplace_matrix_constraints_02b.cc new file mode 100644 index 0000000000..c348650bb9 --- /dev/null +++ b/tests/hp/create_laplace_matrix_constraints_02b.cc @@ -0,0 +1,149 @@ +//---------------------------- create_laplace_matrix_constraints_02b.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_laplace_matrix_constraints_02b.cc --------------------------- + + +// Same as create_laplace_matrix_0[1234]b but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + hp::FECollection element; + element.push_back (FESystem (FE_Q(1), 1, + FE_Q(2), 1)); + element.push_back (FESystem (FE_Q(2), 1, + FE_Q(3), 1)); + element.push_back (FESystem (FE_Q(3), 1, + FE_Q(4), 1)); + + hp::DoFHandler dof(tr); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(); + cell != dof.end(); ++cell) + cell->set_active_fe_index (rand() % element.size()); + + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + MatrixTools:: + create_laplace_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix_ref); + constraints.condense(matrix_ref); + + const Function *const dummy = 0; + MatrixTools::create_laplace_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix, + dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/hp/create_laplace_matrix_constraints_02b/cmp/generic b/tests/hp/create_laplace_matrix_constraints_02b/cmp/generic new file mode 100644 index 0000000000..6a3d892caa --- /dev/null +++ b/tests/hp/create_laplace_matrix_constraints_02b/cmp/generic @@ -0,0 +1,4 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:3d::Matrix error Frobenius: 0 diff --git a/tests/hp/create_mass_matrix_constraints_02.cc b/tests/hp/create_mass_matrix_constraints_02.cc new file mode 100644 index 0000000000..8b858019bd --- /dev/null +++ b/tests/hp/create_mass_matrix_constraints_02.cc @@ -0,0 +1,157 @@ +//---------------------------- create_mass_matrix_constraints_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_02.cc --------------------------- + + +// Same as create_mass_matrix_0[1234] but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + hp::FECollection element; + element.push_back (FESystem (FE_Q(1), 1, + FE_Q(2), 1)); + element.push_back (FESystem (FE_Q(2), 1, + FE_Q(3), 1)); + element.push_back (FESystem (FE_Q(3), 1, + FE_Q(4), 1)); + + hp::DoFHandler dof(tr); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(); + cell != dof.end(); ++cell) + cell->set_active_fe_index (rand() % element.size()); + + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + Functions::ExpFunction rhs_function; + + Vector rhs (dof.n_dofs()), rhs_ref(dof.n_dofs()); + + MatrixTools:: + create_mass_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix_ref, + rhs_function, rhs_ref); + constraints.condense(matrix_ref, rhs_ref); + + const Function *const dummy = 0; + MatrixTools::create_mass_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix, + rhs_function, rhs, dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/hp/create_mass_matrix_constraints_02/cmp/generic b/tests/hp/create_mass_matrix_constraints_02/cmp/generic new file mode 100644 index 0000000000..f29bc0719c --- /dev/null +++ b/tests/hp/create_mass_matrix_constraints_02/cmp/generic @@ -0,0 +1,7 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:1d::RHS vector error l2: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:2d::RHS vector error l2: 0 +DEAL:3d::Matrix error Frobenius: 0 +DEAL:3d::RHS vector error l2: 0 diff --git a/tests/hp/create_mass_matrix_constraints_02b.cc b/tests/hp/create_mass_matrix_constraints_02b.cc new file mode 100644 index 0000000000..27c8597854 --- /dev/null +++ b/tests/hp/create_mass_matrix_constraints_02b.cc @@ -0,0 +1,149 @@ +//---------------------------- create_mass_matrix_constraints_02b.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- create_mass_matrix_constraints_02b.cc --------------------------- + + +// Same as create_mass_matrix_0[1234]b but passing an additional constraint +// matrix and comparing results without constraints + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + + + +template +void +check () +{ + Triangulation tr; + if (dim==2) + GridGenerator::hyper_ball(tr, Point(), 1); + else + GridGenerator::hyper_cube(tr, -1,1); + tr.refine_global (1); + tr.begin_active()->set_refine_flag (); + tr.execute_coarsening_and_refinement (); + if (dim==1) + tr.refine_global(2); + + hp::FECollection element; + element.push_back (FESystem (FE_Q(1), 1, + FE_Q(2), 1)); + element.push_back (FESystem (FE_Q(2), 1, + FE_Q(3), 1)); + element.push_back (FESystem (FE_Q(3), 1, + FE_Q(4), 1)); + + hp::DoFHandler dof(tr); + + for (typename hp::DoFHandler::active_cell_iterator + cell = dof.begin_active(); + cell != dof.end(); ++cell) + cell->set_active_fe_index (rand() % element.size()); + + dof.distribute_dofs(element); + + // use a more complicated mapping + // of the domain and a quadrature + // formula suited to the elements + // we have here + MappingQ mapping (3); + QGauss quadrature(6); + + ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints (dof, constraints); + VectorTools::interpolate_boundary_values(dof, 0, + ConstantFunction(1.,2), + constraints); + constraints.close (); + + // create sparsity pattern. note + // that different components should + // not couple, so use pattern + SparsityPattern sparsity; + { + Table<2,DoFTools::Coupling> mask (2, 2); + mask(0,0) = mask(1,1) = DoFTools::always; + mask(0,1) = mask(1,0) = DoFTools::none; + CompressedSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern (dof, mask, csp, constraints); + sparsity.copy_from(csp); + } + + SparseMatrix matrix, matrix_ref; + matrix.reinit (sparsity); + matrix_ref.reinit (sparsity); + + MatrixTools:: + create_mass_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix_ref); + constraints.condense(matrix_ref); + + const Function *const dummy = 0; + MatrixTools::create_mass_matrix (hp::MappingCollection(mapping), dof, + hp::QCollection(quadrature), matrix, + dummy, constraints); + + // compute reference: need to cancel constrained entries as these will in + // general get different values + matrix.add(-1., matrix_ref); + for (unsigned int i=0; i (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/hp/create_mass_matrix_constraints_02b/cmp/generic b/tests/hp/create_mass_matrix_constraints_02b/cmp/generic new file mode 100644 index 0000000000..6a3d892caa --- /dev/null +++ b/tests/hp/create_mass_matrix_constraints_02b/cmp/generic @@ -0,0 +1,4 @@ + +DEAL:1d::Matrix error Frobenius: 0 +DEAL:2d::Matrix error Frobenius: 0 +DEAL:3d::Matrix error Frobenius: 0 -- 2.39.5