<ol>
+<li> 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.
+<br>
+(Martin Kronbichler, 2013/05/08)
+</li>
+
<li> New: The classes TrilinosWrappers::SparseMatrix and
TrilinosWrappers::BlockSparseMatrix now fully implement vmult and Tvmult with
deal.II's own vector classes Vector<double> and
// forward declarations
template <int dim> class Quadrature;
+class ConstraintMatrix;
+
template<typename number> class Vector;
template<typename number> class FullMatrix;
template<typename number> class SparseMatrix;
* 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 <b>same</b> right hand side for <b>every</b> 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
* <tt>ProblemBase<>::apply_dirichlet_bc</tt> to matrix and right hand
* side.
*
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 <int dim, typename number, int spacedim>
void create_mass_matrix (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, typename number, int spacedim>
void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, typename number, int spacedim>
- void create_mass_matrix (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
+ void create_mass_matrix (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
+ void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *const a = 0,
+ const ConstraintMatrix &constraints = ConstraintMatrix());
+
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, typename number, int spacedim>
void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *const a = 0,
+ const ConstraintMatrix &constraints = ConstraintMatrix());
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, typename number, int spacedim>
void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *const a = 0,
+ const ConstraintMatrix &constraints = ConstraintMatrix());
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
- const hp::DoFHandler<dim,spacedim> &dof,
- const hp::QCollection<dim> &q,
+ void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *const a = 0,
+ const ConstraintMatrix &constraints = ConstraintMatrix());
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
- const hp::QCollection<dim> &q,
+ void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
/**
- * Calls the
- * create_boundary_mass_matrix()
- * function, see above, with
+ * Calls the create_boundary_mass_matrix() function, see above, with
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, int spacedim>
std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, int spacedim>
std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
- * 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,
// const Function<1> * const a = 0);
/**
- * Same function as above, but for hp
- * objects.
+ * Same function as above, but for hp objects.
*/
template <int dim, int spacedim>
std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
/**
- * 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 <int dim, int spacedim>
- void create_laplace_matrix (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, int spacedim>
- void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
- void create_laplace_matrix (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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
* <tt>mapping=MappingQ1@<dim@>()</tt>.
*/
template <int dim, int spacedim>
- void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
- void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
- const hp::DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
- void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
- void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
- const hp::DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *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 <int dim, int spacedim>
- void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+ void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const a = 0);
+ const Function<spacedim> *const a = 0,
+ const ConstraintMatrix &constraints = ConstraintMatrix());
/**
* Exception
template <int dim, class VECTOR, int spacedim>
- void project (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
+ void project (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const Quadrature<dim> &quadrature,
- const Function<spacedim> &function,
+ const Function<spacedim> &function,
VECTOR &vec_result,
const bool enforce_zero_boundary,
const Quadrature<dim-1> &q_boundary,
Vector<double> vec (dof.n_dofs());
SparsityPattern sparsity;
+ // check if constraints are compatible (see below)
+ bool constraints_are_compatible = true;
+ {
+ for (std::map<unsigned int,double>::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);
}
SparseMatrix<double> mass_matrix (sparsity);
Vector<double> 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<spacedim>* 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<double>
- // 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<double> and that wouldn't
+ // necessarily go together with the matrix and other functions
for (unsigned int i=0; i<vec.size(); ++i)
vec_result(i) = vec(i);
}
template <int dim, class VECTOR, int spacedim>
- void project (const DoFHandler<dim,spacedim> &dof,
+ void project (const DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const Quadrature<dim> &quadrature,
- const Function<spacedim> &function,
+ const Function<spacedim> &function,
VECTOR &vec,
const bool enforce_zero_boundary,
const Quadrature<dim-1> &q_boundary,
{ 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)
{
std::vector<unsigned int> 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<dim,spacedim> fe_collection (dof_handler.get_fe());
hp::MappingCollection<dim,spacedim> mapping_collection;
for (unsigned int i=0; i<fe_collection.size(); ++i)
mapping_collection.push_back (mapping);
- // now also create a quadrature
- // collection for the faces of a
- // cell. fill it with a quadrature
- // formula with the support points
- // on faces for each FE
+ // now also create a quadrature collection for the faces of a cell. fill
+ // it with a quadrature formula with the support points on faces for each
+ // FE
hp::QCollection<dim-1> face_quadrature_collection;
for (unsigned int i=0; i<fe_collection.size(); ++i)
{
.push_back (Quadrature<dim-1> (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<dim,spacedim> 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<internal::VectorDoFTuple<dim>,
const FiniteElement<dim> &fe = cell->get_fe ();
typename DH<dim,spacedim>::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<dim> &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<face_dofs.size(); ++i)
if (fe.face_system_to_component_index(i).first ==
first_vector_component)
{
- // find corresponding other
- // components of vector
+ // find corresponding other components of vector
internal::VectorDoFTuple<dim> vector_dofs;
vector_dofs.dof_indices[0] = face_dofs[i];
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<dim> normal_vector
= (cell->face(face_no)->get_boundary()
.normal_vector (cell->face(face_no),
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,
}
}
- // 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)
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
<typename DH<dim,spacedim>::active_cell_iterator,
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),
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();
= 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();
normal[d] = 0;
normal /= normal.norm();
- // then construct constraints
- // from this:
+ // then construct constraints from this:
const internal::VectorDoFTuple<dim> &
dof_indices = same_dof_range[0]->first;
internal::add_constraint (dof_indices, normal,
}
- // 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. <DIM,
- // contributions. Here we
- // do not want to average
- // the normal
- // vectors. Since we have
- // DIM contributions, let's
- // assume (and verify) that
- // they are in fact all
- // linearly independent; in
- // that case, all vector
- // components are
- // constrained and we need
- // to set them to zero
+ // 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. <DIM,
+ // contributions. Here we do not want to average the normal
+ // vectors. Since we have DIM contributions, let's assume (and
+ // verify) that they are in fact all linearly independent; in that
+ // case, all vector components are constrained and we need to set
+ // them to zero
case dim:
{
- // assert that indeed
- // only a single cell has
- // contributed
+ // assert that indeed only a single cell has contributed
Assert (cell_to_normals_map.size() == 1,
ExcInternalError());
- // check linear
- // independence by
- // computing the
- // determinant of the
- // matrix created from
- // all the normal
- // vectors. if they are
- // linearly independent,
- // then the determinant
- // is nonzero. if they
- // are orthogonal, then
- // the matrix is in fact
- // equal to 1 (since they
- // are all unit vectors);
- // make sure the
- // determinant is larger
- // than 1e-3 to avoid
- // cases where cells are
- // degenerate
+ // check linear independence by computing the determinant of the
+ // matrix created from all the normal vectors. if they are
+ // linearly independent, then the determinant is nonzero. if they
+ // are orthogonal, then the matrix is in fact equal to 1 (since
+ // they are all unit vectors); make sure the determinant is larger
+ // than 1e-3 to avoid cases where cells are degenerate
{
Tensor<2,dim> t;
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; i<dim; ++i)
if (!constraints.is_constrained (same_dof_range[0]
->first.dof_indices[i])
}
- // 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<typename DH<dim,spacedim>::active_cell_iterator, std::list<Tensor<1,dim> > >
CellContributions;
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<Tensor<1,dim> > tangential_vectors;
for (typename CellContributions::const_iterator
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:
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();
*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<Tensor<1,dim> >::const_iterator
t = tangential_vectors.begin();
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<dim> &
dof_indices = same_dof_range[0]->first;
internal::add_tangentiality_constraints (dof_indices,
std::vector<unsigned int> dof_indices;
FullMatrix<double> cell_matrix;
dealii::Vector<double> cell_rhs;
+ const ConstraintMatrix *constraints;
};
}
(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; i<dofs_per_cell; ++i)
- (*right_hand_side)(data.dof_indices[i]) += data.cell_rhs(i);
+ data.constraints->distribute_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);
}
}
}
const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
void create_mass_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof,
- q, matrix, coefficient);
+ q, matrix, coefficient, constraints);
}
SparseMatrix<number> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
SparseMatrix<number> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_mass_matrix(StaticMappingQ1<dim,spacedim>::mapping,
- dof, q, matrix, rhs, rhs_vector, coefficient);
+ dof, q, matrix, rhs, rhs_vector, coefficient,
+ constraints);
}
const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
- const hp::QCollection<dim> &q,
+ void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
- create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, matrix, coefficient);
+ create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+ dof, q, matrix, coefficient, constraints);
}
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
- const hp::DoFHandler<dim,spacedim> &dof,
+ void create_mass_matrix (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
template <int dim, typename number, int spacedim>
- void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
- const hp::QCollection<dim> &q,
+ void create_mass_matrix (const hp::DoFHandler<dim,spacedim> &dof,
+ const hp::QCollection<dim> &q,
SparseMatrix<number> &matrix,
- const Function<spacedim> &rhs,
+ const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_mass_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
- matrix, rhs, rhs_vector, coefficient);
+ matrix, rhs, rhs_vector, coefficient, constraints);
}
const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<double> *,
void create_laplace_matrix (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q, matrix, coefficient);
}
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<double> *,
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_laplace_matrix(StaticMappingQ1<dim,spacedim>::mapping, dof, q,
matrix, rhs, rhs_vector, coefficient);
const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<double> *,
void create_laplace_matrix (const hp::DoFHandler<dim,spacedim> &dof,
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q, matrix, coefficient);
}
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
Assert (matrix.m() == dof.n_dofs(),
ExcDimensionMismatch (matrix.m(), dof.n_dofs()));
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<double> *,
SparseMatrix<double> &matrix,
const Function<spacedim> &rhs,
Vector<double> &rhs_vector,
- const Function<spacedim> *const coefficient)
+ const Function<spacedim> *const coefficient,
+ const ConstraintMatrix &constraints)
{
create_laplace_matrix(hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof, q,
matrix, rhs, rhs_vector, coefficient);
const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
(const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping,
SparseMatrix<double> &matrix,
const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,double,deal_II_space_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
SparseMatrix<double> &matrix,
const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
(const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping,
SparseMatrix<float> &matrix,
const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension,float,deal_II_space_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension> &dof,
SparseMatrix<float> &matrix,
const Function<deal_II_space_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_space_dimension> * const coefficient);
+ const Function<deal_II_space_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_boundary_mass_matrix<deal_II_dimension,deal_II_space_dimension>
const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::MappingCollection<deal_II_dimension> &mapping,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension> &dof,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<float> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::MappingCollection<deal_II_dimension> &mapping,
SparseMatrix<float> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_mass_matrix<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension> &dof,
SparseMatrix<float> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
#if deal_II_dimension != 1
(const DoFHandler<deal_II_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const Mapping<deal_II_dimension> &mapping,
const DoFHandler<deal_II_dimension> &dof,
const Quadrature<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const Mapping<deal_II_dimension> &mapping,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const DoFHandler<deal_II_dimension> &dof,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
// hp versions of create_laplace_matrix
template
(const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const hp::MappingCollection<deal_II_dimension> &mapping,
const hp::DoFHandler<deal_II_dimension> &dof,
const hp::QCollection<deal_II_dimension> &q,
SparseMatrix<double> &matrix,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const hp::MappingCollection<deal_II_dimension> &mapping,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
template
void MatrixCreator::create_laplace_matrix<deal_II_dimension>
(const hp::DoFHandler<deal_II_dimension> &dof,
SparseMatrix<double> &matrix,
const Function<deal_II_dimension> &rhs,
Vector<double> &rhs_vector,
- const Function<deal_II_dimension> * const coefficient);
+ const Function<deal_II_dimension> * const coefficient,
+ const ConstraintMatrix &constraints);
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(2);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_01/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_02/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ZeroFunction<dim>(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<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_02b/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(2);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_03/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_04/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(2);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_01/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_02/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ZeroFunction<dim>(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<double> 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<dim> *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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_02b/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(2);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_03/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element(FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1);
+ DoFHandler<dim> 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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> 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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_04/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element;
+ element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+ FE_Q<dim>(3), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+ FE_Q<dim>(4), 1));
+
+ hp::DoFHandler<dim> dof(tr);
+
+ for (typename hp::DoFHandler<dim>::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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+ MatrixTools::
+ create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(quadrature), matrix_ref,
+ rhs_function, rhs_ref);
+ constraints.condense(matrix_ref, rhs_ref);
+
+ const Function<dim> *const dummy = 0;
+ MatrixTools::create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_02/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element;
+ element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+ FE_Q<dim>(3), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+ FE_Q<dim>(4), 1));
+
+ hp::DoFHandler<dim> dof(tr);
+
+ for (typename hp::DoFHandler<dim>::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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ MatrixTools::
+ create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(quadrature), matrix_ref);
+ constraints.condense(matrix_ref);
+
+ const Function<dim> *const dummy = 0;
+ MatrixTools::create_laplace_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_laplace_matrix_constraints_02b/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element;
+ element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+ FE_Q<dim>(3), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+ FE_Q<dim>(4), 1));
+
+ hp::DoFHandler<dim> dof(tr);
+
+ for (typename hp::DoFHandler<dim>::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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ Functions::ExpFunction<dim> rhs_function;
+
+ Vector<double> rhs (dof.n_dofs()), rhs_ref(dof.n_dofs());
+
+ MatrixTools::
+ create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(quadrature), matrix_ref,
+ rhs_function, rhs_ref);
+ constraints.condense(matrix_ref, rhs_ref);
+
+ const Function<dim> *const dummy = 0;
+ MatrixTools::create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+
+ rhs -= rhs_ref;
+ deallog << "RHS vector error l2: " << rhs.l2_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_02/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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 <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tr;
+ if (dim==2)
+ GridGenerator::hyper_ball(tr, Point<dim>(), 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<dim> element;
+ element.push_back (FESystem<dim> (FE_Q<dim>(1), 1,
+ FE_Q<dim>(2), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(2), 1,
+ FE_Q<dim>(3), 1));
+ element.push_back (FESystem<dim> (FE_Q<dim>(3), 1,
+ FE_Q<dim>(4), 1));
+
+ hp::DoFHandler<dim> dof(tr);
+
+ for (typename hp::DoFHandler<dim>::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<dim> mapping (3);
+ QGauss<dim> quadrature(6);
+
+ ConstraintMatrix constraints;
+ DoFTools::make_hanging_node_constraints (dof, constraints);
+ VectorTools::interpolate_boundary_values(dof, 0,
+ ConstantFunction<dim>(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<double> matrix, matrix_ref;
+ matrix.reinit (sparsity);
+ matrix_ref.reinit (sparsity);
+
+ MatrixTools::
+ create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(quadrature), matrix_ref);
+ constraints.condense(matrix_ref);
+
+ const Function<dim> *const dummy = 0;
+ MatrixTools::create_mass_matrix (hp::MappingCollection<dim>(mapping), dof,
+ hp::QCollection<dim>(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<matrix.m(); ++i)
+ if (constraints.is_constrained(i)==true)
+ matrix.diag_element(i) = 0;
+ deallog << "Matrix error Frobenius: " << matrix.frobenius_norm() << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile ("create_mass_matrix_constraints_02b/output");
+ deallog << std::setprecision (3);
+ deallog.attach(logfile);
+ deallog.threshold_double(1e-10);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+DEAL:1d::Matrix error Frobenius: 0
+DEAL:2d::Matrix error Frobenius: 0
+DEAL:3d::Matrix error Frobenius: 0
projects / dealii-svn.git / commitdiff