*/
/**
- * 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
+ * Compute the constraints
+ * resulting from the presence of
+ * hanging nodes. The object into
+ * which these are inserted is
+ * later used to condense the
+ * global system matrix and right
+ * hand side, and to extend the
+ * solution vectors from the true
* degrees of freedom also to the
- * constraint nodes.
+ * constraint nodes. This
+ * function is explained in
+ * detail in the @ref step_6
+ * "step-6" tutorial program and
+ * is used in almost all
+ * following programs as well.
*
- * Since this method does not make sense in
- * one dimension, the function returns
- * immediately. The object is not cleared
- * before use, so you should make sure
- * it containts only constraints you still
- * want; otherwise call the @p clear
- * function.
+ * This function does not clear
+ * the constraint matrix object
+ * before use, in order to allow
+ * adding constraints from
+ * different sources to the same
+ * object. You therefore need to
+ * make sure it contains only
+ * constraints you still want;
+ * otherwise call the
+ * ConstraintMatrix::clear()
+ * function. Likewise, this
+ * function does not close the
+ * object since you may want to
+ * enter other constraints later
+ * on yourself.
*
- * To condense a given sparsity pattern,
- * use ConstraintMatrix@p ::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.
- */
- static void
- make_hanging_node_constraints (const DoFHandler<1> &dof_handler,
- ConstraintMatrix &constraints);
-
- /**
- * Declaration of same function
- * for different space dimension.
- */
- static void
- make_hanging_node_constraints (const DoFHandler<2> &dof_handler,
- ConstraintMatrix &constraints);
-
- /**
- * Declaration of same function
- * for different space dimension.
- */
- static void
- make_hanging_node_constraints (const DoFHandler<3> &dof_handler,
- ConstraintMatrix &constraints);
-
- /**
- * Declaration of same function
- * for hp::DoFHandler
- */
- static void
- make_hanging_node_constraints (const hp::DoFHandler<1> &dof_handler,
- ConstraintMatrix &constraints);
-
- /**
- * Declaration of same function
- * for different space dimension.
- */
- static void
- make_hanging_node_constraints (const hp::DoFHandler<2> &dof_handler,
- ConstraintMatrix &constraints);
-
- /**
- * Declaration of same function
- * for different space dimension.
+ * In the hp-case, i.e. when the
+ * argument is of type
+ * hp::DoFHandler, we consider
+ * constraints due to different
+ * finite elements used on two
+ * sides of a face between cells
+ * as hanging nodes as well. In
+ * other words, for hp finite
+ * elements, this function
+ * computes all constraints due
+ * to differing mesh sizes (h) or
+ * polynomial degrees (p) between
+ * adjacent cells.
+ *
+ * The template argument (and by
+ * consequence the type of the
+ * first argument to this
+ * function) can be either a
+ * ::DoFHandler, hp::DoFHandler,
+ * or MGDoFHandler.
*/
+ template <class DH>
static void
- make_hanging_node_constraints (const hp::DoFHandler<3> &dof_handler,
- ConstraintMatrix &constraints);
-
+ make_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints);
//@}
/**
-#if deal_II_dimension == 1
-
-void DoFTools::make_hanging_node_constraints (
- const DoFHandler<1> &,
- ConstraintMatrix &)
+namespace internal
{
- // nothing to be done here
-}
+ namespace DoFTools
+ {
+ template <int> class int2type
+ {};
-void DoFTools::make_hanging_node_constraints (
- const hp::DoFHandler<1> &,
- ConstraintMatrix &)
-{
- // we may have to compute
- // constraints for vertices
- Assert (false, ExcNotImplemented());
-}
+#if deal_II_dimension == 1
+ void
+ do_make_hanging_node_constraints (const ::DoFHandler<1> &,
+ ConstraintMatrix &,
+ int2type<1>)
+ {
+ // nothing to do for regular
+ // dof handlers in 1d
+ }
-#endif
+ void
+ do_make_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+ ConstraintMatrix &/*constraints*/,
+ int2type<1>)
+ {
+ // we may have to compute
+ // constraints for
+ // vertices. gotta think about
+ // that a bit more
+//TODO[WB]: think about what to do here...
+ }
+#endif
#if deal_II_dimension == 2
-
-void DoFTools::make_hanging_node_constraints (
- const DoFHandler<2> &dof_handler,
- ConstraintMatrix &constraints)
-{
- const unsigned int dim = 2;
-
- const FiniteElement<dim> &fe = dof_handler.get_fe();
+ template <class DH>
+ void
+ do_make_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints,
+ int2type<2>)
+ {
+ const unsigned int dim = 2;
- // have space for the degrees of
- // freedom on mother and child
- // lines
- const unsigned int n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line,
- n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
-
- std::vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
- std::vector<unsigned int> dofs_on_children(n_dofs_on_children);
-
- Assert(n_dofs_on_mother == fe.constraints().n(),
- ExcDimensionMismatch(n_dofs_on_mother,
- fe.constraints().n()));
- Assert(n_dofs_on_children == fe.constraints().m(),
- ExcDimensionMismatch(n_dofs_on_children,
- fe.constraints().m()));
-
- // loop over all lines; only on
- // lines there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
- // constraints
- //
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
- DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ std::vector<unsigned int> dofs_on_mother;
+ std::vector<unsigned int> dofs_on_children;
+
+ // loop over all lines; only on
+ // lines there can be constraints.
+ // We do so by looping over all
+ // active cells and checking
+ // whether any of the faces are
+ // refined which can only be from
+ // the neighboring cell because
+ // this one is active. In that
+ // case, the face is subject to
+ // constraints
+ //
+ // note that even though we may
+ // visit a face twice if the
+ // neighboring cells are equally
+ // refined, we can only visit each
+ // face with hanging nodes once
+ typename DH::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->has_children())
- {
- const DoFHandler<dim>::line_iterator line = cell->face(face);
-
- // fill the dofs indices. Use same
- // enumeration scheme as in
- // @p{FiniteElement::constraints()}
- unsigned int next_index = 0;
- for (unsigned int vertex=0; vertex<2; ++vertex)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_mother[next_index++] = line->vertex_dof_index(vertex,dof);
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_mother[next_index++] = line->dof_index(dof);
- Assert (next_index == dofs_on_mother.size(),
- ExcInternalError());
-
- next_index = 0;
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++] = line->child(0)->vertex_dof_index(1,dof);
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++] = line->child(child)->dof_index(dof);
- Assert (next_index == dofs_on_children.size(),
- ExcInternalError());
-
- // for each row in the constraint
- // matrix for this line:
- for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
{
- constraints.add_line (dofs_on_children[row]);
- for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
- constraints.add_entry (dofs_on_children[row],
- dofs_on_mother[i],
- fe.constraints()(row,i));
- }
- }
-}
-
+ // so now we've found a
+ // face of an active
+ // cell that has
+ // children. that means
+ // that there are
+ // hanging nodes here.
+
+ // in any case, faces
+ // can have at most two
+ // active fe indices,
+ // but here the face
+ // can have only one
+ // (namely the same as
+ // that from the cell
+ // we're sitting on),
+ // and each of the
+ // children can have
+ // only one as
+ // well. check this
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcInternalError());
+ Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+ == true,
+ ExcInternalError());
+ for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+ Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+ ExcInternalError());
-void DoFTools::make_hanging_node_constraints (
- const hp::DoFHandler<2> &/*dof_handler*/,
- ConstraintMatrix &/*constraints*/)
-{
-//TODO[?]: Implement (required for continuous elements)
- Assert (false, ExcNotImplemented());
-}
+ // right now, all that
+ // is implemented is
+ // the case that both
+ // sides use the same
+ // fe
+ for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+ Assert (cell->face(face)->child(c)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+
+ // ok, start up the work
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int fe_index = cell->active_fe_index();
+
+ const unsigned int
+ n_dofs_on_mother = 2*fe.dofs_per_vertex + fe.dofs_per_line,
+ n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
+ dofs_on_mother.resize (n_dofs_on_mother);
+ dofs_on_children.resize (n_dofs_on_children);
+
+ Assert(n_dofs_on_mother == fe.constraints().n(),
+ ExcDimensionMismatch(n_dofs_on_mother,
+ fe.constraints().n()));
+ Assert(n_dofs_on_children == fe.constraints().m(),
+ ExcDimensionMismatch(n_dofs_on_children,
+ fe.constraints().m()));
+
+ const typename DH::line_iterator this_face = cell->face(face);
+
+ // fill the dofs indices. Use same
+ // enumeration scheme as in
+ // @p{FiniteElement::constraints()}
+ unsigned int next_index = 0;
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+ fe_index);
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+ Assert (next_index == dofs_on_mother.size(),
+ ExcInternalError());
+
+ next_index = 0;
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->vertex_dof_index(1,dof,fe_index);
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(child)->dof_index(dof, fe_index);
+ Assert (next_index == dofs_on_children.size(),
+ ExcInternalError());
+
+ // for each row in the constraint
+ // matrix for this line:
+ for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
+ {
+ constraints.add_line (dofs_on_children[row]);
+ for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+ constraints.add_entry (dofs_on_children[row],
+ dofs_on_mother[i],
+ fe.constraints()(row,i));
+ }
+ }
+ else
+ {
+ // this face has no
+ // children, but it
+ // could still be that
+ // it is shared by two
+ // cells that use a
+ // different fe index
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcInternalError());
+ }
+ }
#endif
#if deal_II_dimension == 3
-
-void DoFTools::make_hanging_node_constraints (
- const DoFHandler<3> &dof_handler,
- ConstraintMatrix &constraints)
-{
- const unsigned int dim = 3;
-
- const FiniteElement<dim> &fe = dof_handler.get_fe();
+ template <class DH>
+ void
+ do_make_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints,
+ int2type<3>)
+ {
+ const unsigned int dim = 3;
- // have space for the degrees of
- // freedom on mother and child
- // lines
- const unsigned int
- n_dofs_on_mother = (4*fe.dofs_per_vertex+
- 4*fe.dofs_per_line+
- fe.dofs_per_quad),
- n_dofs_on_children = (5*fe.dofs_per_vertex+
- 12*fe.dofs_per_line+
- 4*fe.dofs_per_quad);
-
- std::vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
- std::vector<unsigned int> dofs_on_children(n_dofs_on_children);
-
- Assert(n_dofs_on_mother == fe.constraints().n(),
- ExcDimensionMismatch(n_dofs_on_mother,
- fe.constraints().n()));
- Assert(n_dofs_on_children == fe.constraints().m(),
- ExcDimensionMismatch(n_dofs_on_children,
- fe.constraints().m()));
-
- // loop over all lines; only on
- // lines there can be constraints.
- // We do so by looping over all
- // active cells and checking
- // whether any of the faces are
- // refined which can only be from
- // the neighboring cell because
- // this one is active. In that
- // case, the face is subject to
- // constraints
- //
- // note that even though we may
- // visit a face twice if the
- // neighboring cells are equally
- // refined, we can only visit each
- // face with hanging nodes once
- DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
- endc = dof_handler.end();
- for (; cell!=endc; ++cell)
- for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- if (cell->face(f)->has_children())
- {
- const DoFHandler<dim>::face_iterator face = cell->face(f);
-
- // fill the dofs indices. Use same
- // enumeration scheme as in
- // @p{FiniteElement::constraints()}
- unsigned int next_index = 0;
- for (unsigned int vertex=0; vertex<4; ++vertex)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_mother[next_index++] = face->vertex_dof_index(vertex,dof);
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
- dofs_on_mother[next_index++] = face->line(line)->dof_index(dof);
- for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- dofs_on_mother[next_index++] = face->dof_index(dof);
- Assert (next_index == dofs_on_mother.size(),
- ExcInternalError());
+ std::vector<unsigned int> dofs_on_mother;
+ std::vector<unsigned int> dofs_on_children;
+
+ // loop over all quads; only on
+ // quads there can be constraints.
+ // We do so by looping over all
+ // active cells and checking
+ // whether any of the faces are
+ // refined which can only be from
+ // the neighboring cell because
+ // this one is active. In that
+ // case, the face is subject to
+ // constraints
+ //
+ // note that even though we may
+ // visit a face twice if the
+ // neighboring cells are equally
+ // refined, we can only visit each
+ // face with hanging nodes once
+ typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->has_children())
+ {
+ // so now we've found a
+ // face of an active
+ // cell that has
+ // children. that means
+ // that there are
+ // hanging nodes here.
+
+ // in any case, faces
+ // can have at most two
+ // active fe indices,
+ // but here the face
+ // can have only one
+ // (namely the same as
+ // that from the cell
+ // we're sitting on),
+ // and each of the
+ // children can have
+ // only one as
+ // well. check this
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcInternalError());
+ Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+ == true,
+ ExcInternalError());
+ for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+ Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+ ExcInternalError());
+
+ // right now, all that
+ // is implemented is
+ // the case that both
+ // sides use the same
+ // fe, and not only
+ // that but also that
+ // all lines bounding
+ // this face and the
+ // children have the
+ // same fe
+ for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+ {
+ Assert (cell->face(face)->child(c)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ for (unsigned int e=0; e<4; ++e)
+ {
+ Assert (cell->face(face)->child(c)->line(e)
+ ->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)->child(c)->line(e)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ }
+ }
+ for (unsigned int e=0; e<4; ++e)
+ {
+ Assert (cell->face(face)->line(e)
+ ->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)->line(e)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcNotImplemented());
+ }
+
+ // ok, start up the work
+ const FiniteElement<dim> &fe = cell->get_fe();
+ const unsigned int fe_index = cell->active_fe_index();
+
+ const unsigned int
+ n_dofs_on_mother = (4*fe.dofs_per_vertex+
+ 4*fe.dofs_per_line+
+ fe.dofs_per_quad),
+ n_dofs_on_children = (5*fe.dofs_per_vertex+
+ 12*fe.dofs_per_line+
+ 4*fe.dofs_per_quad);
+
+ dofs_on_mother.resize (n_dofs_on_mother);
+ dofs_on_children.resize (n_dofs_on_children);
+
+ Assert(n_dofs_on_mother == fe.constraints().n(),
+ ExcDimensionMismatch(n_dofs_on_mother,
+ fe.constraints().n()));
+ Assert(n_dofs_on_children == fe.constraints().m(),
+ ExcDimensionMismatch(n_dofs_on_children,
+ fe.constraints().m()));
+
+ const typename DH::face_iterator this_face = cell->face(face);
- next_index = 0;
-
- // assert some consistency
- // assumptions
- Assert ((face->child(0)->vertex_index(3) ==
- face->child(1)->vertex_index(2)) &&
- (face->child(0)->vertex_index(3) ==
- face->child(2)->vertex_index(1)) &&
- (face->child(0)->vertex_index(3) ==
- face->child(3)->vertex_index(0)),
- ExcInternalError());
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++]
- = face->child(0)->vertex_dof_index(3,dof);
+ // fill the dofs indices. Use same
+ // enumeration scheme as in
+ // @p{FiniteElement::constraints()}
+ unsigned int next_index = 0;
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+ fe_index);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_mother[next_index++]
+ = this_face->line(line)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+ dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+ Assert (next_index == dofs_on_mother.size(),
+ ExcInternalError());
- // dof numbers on the centers of
- // the lines bounding this face
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
- dofs_on_children[next_index++]
- = face->line(line)->child(0)->vertex_dof_index(1,dof);
+ next_index = 0;
+
+ // assert some consistency
+ // assumptions
+ Assert ((this_face->child(0)->vertex_index(3) ==
+ this_face->child(1)->vertex_index(2)) &&
+ (this_face->child(0)->vertex_index(3) ==
+ this_face->child(2)->vertex_index(1)) &&
+ (this_face->child(0)->vertex_index(3) ==
+ this_face->child(3)->vertex_index(0)),
+ ExcInternalError());
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->vertex_dof_index(3,dof);
- // next the dofs on the lines interior
- // to the face; the order of these
- // lines is laid down in the
- // FiniteElement class documentation
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = face->child(0)->line(1)->dof_index(dof);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = face->child(2)->line(1)->dof_index(dof);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = face->child(0)->line(3)->dof_index(dof);
- for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
- dofs_on_children[next_index++]
- = face->child(1)->line(3)->dof_index(dof);
+ // dof numbers on the centers of
+ // the lines bounding this face
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
- // dofs on the bordering lines
- for (unsigned int line=0; line<4; ++line)
- for (unsigned int child=0; child<2; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ // next the dofs on the lines interior
+ // to the face; the order of these
+ // lines is laid down in the
+ // FiniteElement class documentation
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(0)->line(1)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(2)->line(1)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
dofs_on_children[next_index++]
- = face->line(line)->child(child)->dof_index(dof);
+ = this_face->child(0)->line(3)->dof_index(dof, fe_index);
+ for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(1)->line(3)->dof_index(dof, fe_index);
- // finally, for the dofs interior
- // to the four child faces
- for (unsigned int child=0; child<4; ++child)
- for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- dofs_on_children[next_index++]
- = face->child(child)->dof_index(dof);
- Assert (next_index == dofs_on_children.size(),
- ExcInternalError());
+ // dofs on the bordering lines
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int child=0; child<2; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->line(line)->child(child)->dof_index(dof, fe_index);
- // for each row in the constraint
- // matrix for this line:
- for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
- {
- constraints.add_line (dofs_on_children[row]);
- for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
- constraints.add_entry (dofs_on_children[row],
- dofs_on_mother[i],
- fe.constraints()(row,i));
+ // finally, for the dofs interior
+ // to the four child faces
+ for (unsigned int child=0; child<4; ++child)
+ for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+ dofs_on_children[next_index++]
+ = this_face->child(child)->dof_index(dof, fe_index);
+ Assert (next_index == dofs_on_children.size(),
+ ExcInternalError());
+
+ // for each row in the constraint
+ // matrix for this line:
+ for (unsigned int row=0; row!=dofs_on_children.size(); ++row)
+ {
+ constraints.add_line (dofs_on_children[row]);
+ for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+ constraints.add_entry (dofs_on_children[row],
+ dofs_on_mother[i],
+ fe.constraints()(row,i));
+ }
}
- }
+ else
+ {
+ // this face has no
+ // children, but it
+ // could still be that
+ // it is shared by two
+ // cells that use a
+ // different fe index
+ Assert (cell->face(face)->n_active_fe_indices() == 1,
+ ExcNotImplemented());
+ Assert (cell->face(face)
+ ->fe_index_is_active(cell->active_fe_index()) == true,
+ ExcInternalError());
+ }
+ }
+#endif
+
+ }
}
-void DoFTools::make_hanging_node_constraints (
- const hp::DoFHandler<3> &/*dof_handler*/,
- ConstraintMatrix &/*constraints*/)
+
+template <class DH>
+void
+DoFTools::make_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints)
{
-//TODO:[?] Implement (required for continuous elements)
- Assert (false, ExcNotImplemented());
+ // distribute this one template
+ // function to the appropriate
+ // functions in 1d, 2d, and 3d
+ do_make_hanging_node_constraints (dof_handler,
+ constraints,
+ internal::DoFTools::int2type<DH::dimension>());
}
-#endif
-
template <class DH, typename Number>
#endif
+template
+void
+DoFTools::make_hanging_node_constraints (const DoFHandler<deal_II_dimension> &dof_handler,
+ ConstraintMatrix &constraints);
+template
+void
+DoFTools::make_hanging_node_constraints (const MGDoFHandler<deal_II_dimension> &dof_handler,
+ ConstraintMatrix &constraints);
+template
+void
+DoFTools::make_hanging_node_constraints (const hp::DoFHandler<deal_II_dimension> &dof_handler,
+ ConstraintMatrix &constraints);
+
+
template
void