// computing constraints
cout << " Computing constraints..." << endl;
ConstraintMatrix constraints;
- dof->make_constraint_matrix (constraints);
+ dof->make_hanging_node_constraints (constraints);
+ constraints.close ();
constraints.condense (sparsity);
cout << " Writing condensed sparsity pattern..." << endl;
#include <numerics/assembler.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
+#include <numerics/multigrid.h>
+#include <numerics/mg_smoother.h>
#include <lac/vector.h>
#include <lac/sparsematrix.h>
// make up sparsity pattern and
// compress with constraints
constraints.clear ();
- dof->DoFHandler<dim>::make_constraint_matrix (constraints);
+ dof->DoFHandler<dim>::make_hanging_node_constraints (constraints);
+ constraints.close ();
dof->DoFHandler<dim>::make_sparsity_pattern (system_sparsity);
constraints.condense (system_sparsity);
+ MGTransferPrebuilt p;
+ p.build_matrices (*dof);
+// MGSmoother smoother(*dof);
+
+
// reinite system matrix
system_matrix.reinit (system_sparsity);
// reinit solution vector, preset
solution.reinit (dof->DoFHandler<dim>::n_dofs());
// create assembler
- AssemblerData<dim> data (*dof,
- true, true, //assemble matrix and rhs
- system_matrix,
- right_hand_side,
- quadrature,
- update_flags);
+ Assembler<dim>::AssemblerData data (*dof,
+ true, true, //assemble matrix and rhs
+ system_matrix,
+ right_hand_side,
+ quadrature,
+ update_flags);
active_assemble_iterator assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
cout << ">" << endl;
cout << " Making grid... ";
- tria->create_hyper_ball ();
- HyperBallBoundary<dim> boundary_description;
- tria->set_boundary (&boundary_description);
+ tria->create_hypercube ();
+ tria->refine_global (level+1);
tria->begin_active()->set_refine_flag();
(++(++(tria->begin_active())))->set_refine_flag();
tria->execute_coarsening_and_refinement ();
- tria->refine_global (level);
cout << tria->n_active_cells() << " active cells." << endl;
rhs = new RHSPoly<dim>();
cout << endl;
+ const unsigned int n_dofs = dof->DoFHandler<dim>::n_dofs();
+ dof->clear ();
+ tria->set_boundary (0);
delete fe;
delete quadrature;
delete boundary_quadrature;
- return dof->DoFHandler<dim>::n_dofs();
+ return n_dofs;
};
void renumber_dofs (const vector<int> &new_numbers);
/**
- * Make up the constraint matrix which
- * is used to condensate the global
+ * Make up the constraints which
+ * is result from the use of hanging
+ * nodes. The object into which these
+ * are inserted is later
+ * used to condensate the global
* system matrices and to prolong
* the solution vectors from the true
* degrees of freedom also to the
*
* Since this method does not make sense in
* one dimension, the function returns
- * immediately after clearing the
- * constraint matrix.
- * For more than one dimension, the matrix
- * is cleared before usage. The constraint
- * matrix is closed anyway, no matter of the
- * dimension.
+ * immediately. The object is not cleared
+ * before use, so you should make sure
+ * it containts only constraints you still
+ * want; otherwise call the #clear#
+ * function.
*
* To condense a given sparsity pattern,
* use #ConstraintMatrix::condense#.
+ * Before doing so, you need to close
+ * the constraint object, which must be
+ * done after all constraints are entered.
+ * This function does not close the object
+ * since you may want to enter other
+ * constraints later on yourself.
*
* This function uses the user flags for
* the faces. If you need the user flags,
* store them beforehand.
*/
- void make_constraint_matrix (ConstraintMatrix &) const;
+ void make_hanging_node_constraints (ConstraintMatrix &) const;
/**
* $\ldots, u_1^m, u_2^m, u_1^{m+1}, u_2^{m+1},\ldots$ with $m$ denoting the
* $m$th basis function, or $\ldots, u_1^m, u_1^{m+1}, u_1^{m+2}, \ldots,
* u_2^m, u_2^{m+1}, u_2^{m+2}, \ldots$, respectively). Likewise, the
- * constraint matrix returned by #DoFHandler::make_constraint_matrix ()# is then
+ * constraint matrix returned by #DoFHandler::make_hanging_node_constraints ()#
+ * is then
* to be understood as a block matrix.
*
* The storage format of the degrees of freedom indices (short: DoF indices) is
-/**
- * Structure to be passed upon
- * construction of an assembler object
- * through the iterator object. See
- * \Ref{TriaRawIterator} for a discussion
- * of this mechanism.
- */
-template <int dim>
-struct AssemblerData {
- /**
- * Constructor.
- */
- AssemblerData (const DoFHandler<dim> &dof,
- const bool assemble_matrix,
- const bool assemble_rhs,
- SparseMatrix<double> &matrix,
- Vector<double> &rhs_vector,
- const Quadrature<dim> &quadrature,
- const UpdateFlags &update_flags);
-
- /**
- * Pointer to the dof handler object
- * to be used to iterate on.
- */
- const DoFHandler<dim> &dof;
-
- /**
- * Flags whether to assemble the matrix
- * and right hand sides.
- */
- const bool assemble_matrix, assemble_rhs;
-
- /**
- * Pointer to the matrix to be assembled
- * by this object. Elements are summed
- * up by the assembler, so you may want
- * to clear this object (set all entries
- * to zero) before use.
- */
- SparseMatrix<double> &matrix;
-
- /**
- * Pointer to the vector to be assembled
- * by this object. Elements are summed
- * up by the assembler, so you may want
- * to clear this object (set all entries
- * to zero) before use.
- */
- Vector<double> &rhs_vector;
-
- /**
- * Pointer to a quadrature object to be
- * used for this assemblage process.
- */
- const Quadrature<dim> &quadrature;
-
- /**
- * Store which of the fields of the
- * FEValues object need to be reinitialized
- * on each cell.
- */
- const UpdateFlags update_flags;
-};
-
-
-
-
/**
* An #Assembler# is a specialized version of a #DoFCellAccessor# which adds
* functionality to assemble global matrices and vectors from cell base ones.
template <int dim>
class Assembler : public DoFCellAccessor<dim> {
public:
+
+ /**
+ * Structure to be passed upon
+ * construction of an assembler object
+ * through the iterator object. See
+ * \Ref{TriaRawIterator} for a discussion
+ * of this mechanism.
+ */
+ struct AssemblerData {
+ /**
+ * Constructor.
+ */
+ AssemblerData (const DoFHandler<dim> &dof,
+ const bool assemble_matrix,
+ const bool assemble_rhs,
+ SparseMatrix<double> &matrix,
+ Vector<double> &rhs_vector,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags &update_flags);
+
+ /**
+ * Pointer to the dof handler object
+ * to be used to iterate on.
+ */
+ const DoFHandler<dim> &dof;
+
+ /**
+ * Flags whether to assemble the matrix
+ * and right hand sides.
+ */
+ const bool assemble_matrix, assemble_rhs;
+
+ /**
+ * Pointer to the matrix to be assembled
+ * by this object. Elements are summed
+ * up by the assembler, so you may want
+ * to clear this object (set all entries
+ * to zero) before use.
+ */
+ SparseMatrix<double> &matrix;
+
+ /**
+ * Pointer to the vector to be assembled
+ * by this object. Elements are summed
+ * up by the assembler, so you may want
+ * to clear this object (set all entries
+ * to zero) before use.
+ */
+ Vector<double> &rhs_vector;
+
+ /**
+ * Pointer to a quadrature object to be
+ * used for this assemblage process.
+ */
+ const Quadrature<dim> &quadrature;
+
+ /**
+ * Store which of the fields of the
+ * FEValues object need to be reinitialized
+ * on each cell.
+ */
+ const UpdateFlags update_flags;
+ };
+
+
+
/**
* Declare the data type that this accessor
* class expects to get passed from the
* iterator classes.
*/
- typedef AssemblerData<dim> AccessorData;
+ typedef AssemblerData AccessorData;
/**
* Default constructor, unused thus not
if (use_constraints)
{
ConstraintMatrix constraints;
- make_constraint_matrix (constraints);
+ make_hanging_node_constraints (constraints);
+ constraints.close ();
constraints.condense (sparsity);
};
#if deal_II_dimension == 1
template <>
-void DoFHandler<1>::make_constraint_matrix (ConstraintMatrix &cm) const {
- cm.clear ();
- cm.close ();
-};
+void DoFHandler<1>::make_hanging_node_constraints (ConstraintMatrix &) const {};
#endif
#if deal_II_dimension == 2
template <>
-void DoFHandler<2>::make_constraint_matrix (ConstraintMatrix &constraints) const {
+void DoFHandler<2>::make_hanging_node_constraints (ConstraintMatrix &constraints) const {
const unsigned int dim = 2;
- constraints.clear ();
-
// first mark all faces which are subject
// to constraints. We do so by looping
// over all active cells and checking
selected_fe->constraints()(row,i));
};
};
-
- constraints.close ();
};
#endif
#if deal_II_dimension == 3
template <>
-void DoFHandler<3>::make_constraint_matrix (ConstraintMatrix &constraints) const {
+void DoFHandler<3>::make_hanging_node_constraints (ConstraintMatrix &constraints) const {
const unsigned int dim = 3;
- constraints.clear ();
-
// first mark all faces which are subject
// to constraints. We do so by looping
// over all active cells and checking
selected_fe->constraints()(row,i));
};
};
-
- constraints.close ();
};
#endif
template <int dim>
-AssemblerData<dim>::AssemblerData (const DoFHandler<dim> &dof,
- const bool assemble_matrix,
- const bool assemble_rhs,
- SparseMatrix<double> &matrix,
- Vector<double> &rhs_vector,
- const Quadrature<dim> &quadrature,
- const UpdateFlags &update_flags) :
+Assembler<dim>::AssemblerData::AssemblerData (const DoFHandler<dim> &dof,
+ const bool assemble_matrix,
+ const bool assemble_rhs,
+ SparseMatrix<double> &matrix,
+ Vector<double> &rhs_vector,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags &update_flags) :
dof(dof),
assemble_matrix(assemble_matrix),
assemble_rhs(assemble_rhs),
template <int dim>
-Assembler<dim>::Assembler (Triangulation<dim> *tria,
- const int level,
- const int index,
- const AssemblerData<dim> *local_data) :
+Assembler<dim>::Assembler (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const AssemblerData *local_data) :
DoFCellAccessor<dim> (tria,level,index, &local_data->dof),
cell_matrix (dof_handler->get_fe().total_dofs),
cell_vector (Vector<double>(dof_handler->get_fe().total_dofs)),
// explicit instantiations
template class Equation<deal_II_dimension>;
template class Assembler<deal_II_dimension>;
-template class AssemblerData<deal_II_dimension>;
template class TriaRawIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
template class TriaIterator<deal_II_dimension,Assembler<deal_II_dimension> >;
// make up sparsity pattern and
// compress with constraints
constraints.clear ();
- dof_handler->make_constraint_matrix (constraints);
+ dof_handler->make_hanging_node_constraints (constraints);
+ constraints.close ();
dof_handler->make_sparsity_pattern (system_sparsity);
constraints.condense (system_sparsity);
solution.reinit (dof_handler->n_dofs());
// create assembler
- AssemblerData<dim> data (*dof_handler,
- true, true, //assemble matrix and rhs
- system_matrix,
- right_hand_side,
- quadrature,
- update_flags);
+ Assembler<dim>::AssemblerData data (*dof_handler,
+ true, true, //assemble matrix and rhs
+ system_matrix,
+ right_hand_side,
+ quadrature,
+ update_flags);
active_assemble_iterator assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
UpdateFlags update_flags = update_JxW_values;
if (a != 0)
update_flags = UpdateFlags (update_flags | update_q_points);
- const AssemblerData<dim> data (dof,
- true, false, // assemble matrix but not rhs
- matrix, dummy,
- q, update_flags);
+ const Assembler<dim>::AssemblerData data (dof,
+ true, false, // assemble matrix but not rhs
+ matrix, dummy,
+ q, update_flags);
TriaActiveIterator<dim, Assembler<dim> >
assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
dof.get_tria().begin_active()->level(),
const Function<dim> * const a) {
UpdateFlags update_flags = UpdateFlags(update_q_points |
update_JxW_values);
- const AssemblerData<dim> data (dof,
- true, true,
- matrix, rhs_vector,
- q, update_flags);
+ const Assembler<dim>::AssemblerData data (dof,
+ true, true,
+ matrix, rhs_vector,
+ q, update_flags);
TriaActiveIterator<dim, Assembler<dim> >
assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
dof.get_tria().begin_active()->level(),
update_JxW_values);
if (a != 0)
update_flags = UpdateFlags(update_flags | update_q_points);
- const AssemblerData<dim> data (dof,
- true, false, // assemble matrix but not rhs
- matrix, dummy,
- q, update_flags);
+ const Assembler<dim>::AssemblerData data (dof,
+ true, false, // assemble matrix but not rhs
+ matrix, dummy,
+ q, update_flags);
TriaActiveIterator<dim, Assembler<dim> >
assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
dof.get_tria().begin_active()->level(),
UpdateFlags update_flags = UpdateFlags(update_q_points |
update_gradients |
update_JxW_values);
- const AssemblerData<dim> data (dof,
- true, true,
- matrix, rhs_vector,
- q, update_flags);
+ const Assembler<dim>::AssemblerData data (dof,
+ true, true,
+ matrix, rhs_vector,
+ q, update_flags);
TriaActiveIterator<dim, Assembler<dim> >
assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
dof.get_tria().begin_active()->level(),
UpdateFlags update_flags = UpdateFlags(update_q_points |
update_JxW_values);
SparseMatrix<double> dummy;
- const AssemblerData<dim> data (dof,
- false, true,
- dummy, rhs_vector,
- quadrature, update_flags);
+ const Assembler<dim>::AssemblerData data (dof,
+ false, true,
+ dummy, rhs_vector,
+ quadrature, update_flags);
TriaActiveIterator<dim, Assembler<dim> >
assembler (const_cast<Triangulation<dim>*>(&dof.get_tria()),
dof.get_tria().begin_active()->level(),
VectorTools<dim>::project_boundary_values (const DoFHandler<dim> &dof,
const FunctionMap &boundary_functions,
const Quadrature<dim-1> &q,
- map<int,double> &boundary_values) {
+ map<int,double> &boundary_values) {
vector<int> dof_to_boundary_mapping;
dof.map_dof_to_boundary_indices (boundary_functions, dof_to_boundary_mapping);
// computing constraints
cout << " Computing constraints..." << endl;
ConstraintMatrix constraints;
- dof->make_constraint_matrix (constraints);
+ dof->make_hanging_node_constraints (constraints);
+ constraints.close ();
constraints.condense (sparsity);
cout << " Writing condensed sparsity pattern..." << endl;
#include <numerics/assembler.h>
#include <numerics/vectors.h>
#include <numerics/matrices.h>
+#include <numerics/multigrid.h>
+#include <numerics/mg_smoother.h>
#include <lac/vector.h>
#include <lac/sparsematrix.h>
// make up sparsity pattern and
// compress with constraints
constraints.clear ();
- dof->DoFHandler<dim>::make_constraint_matrix (constraints);
+ dof->DoFHandler<dim>::make_hanging_node_constraints (constraints);
+ constraints.close ();
dof->DoFHandler<dim>::make_sparsity_pattern (system_sparsity);
constraints.condense (system_sparsity);
+ MGTransferPrebuilt p;
+ p.build_matrices (*dof);
+// MGSmoother smoother(*dof);
+
+
// reinite system matrix
system_matrix.reinit (system_sparsity);
// reinit solution vector, preset
solution.reinit (dof->DoFHandler<dim>::n_dofs());
// create assembler
- AssemblerData<dim> data (*dof,
- true, true, //assemble matrix and rhs
- system_matrix,
- right_hand_side,
- quadrature,
- update_flags);
+ Assembler<dim>::AssemblerData data (*dof,
+ true, true, //assemble matrix and rhs
+ system_matrix,
+ right_hand_side,
+ quadrature,
+ update_flags);
active_assemble_iterator assembler (tria,
tria->begin_active()->level(),
tria->begin_active()->index(),
cout << ">" << endl;
cout << " Making grid... ";
- tria->create_hyper_ball ();
- HyperBallBoundary<dim> boundary_description;
- tria->set_boundary (&boundary_description);
+ tria->create_hypercube ();
+ tria->refine_global (level+1);
tria->begin_active()->set_refine_flag();
(++(++(tria->begin_active())))->set_refine_flag();
tria->execute_coarsening_and_refinement ();
- tria->refine_global (level);
cout << tria->n_active_cells() << " active cells." << endl;
rhs = new RHSPoly<dim>();
cout << endl;
+ const unsigned int n_dofs = dof->DoFHandler<dim>::n_dofs();
+ dof->clear ();
+ tria->set_boundary (0);
delete fe;
delete quadrature;
delete boundary_quadrature;
- return dof->DoFHandler<dim>::n_dofs();
+ return n_dofs;
};