#if deal_II_dimension == 1
static
void
- do_make_hanging_node_constraints (const ::DoFHandler<1> &,
- ConstraintMatrix &,
- ::internal::int2type<1>)
+ make_hp_hanging_node_constraints (const ::DoFHandler<1> &,
+ ConstraintMatrix &)
{
// nothing to do for regular
// dof handlers in 1d
}
+
+ static
+ void
+ make_oldstyle_hanging_node_constraints (const ::DoFHandler<1> &,
+ ConstraintMatrix &,
+ ::internal::int2type<1>)
+ {
+ // nothing to do for regular
+ // dof handlers in 1d
+ }
+
+
+ static
+ void
+ make_hp_hanging_node_constraints (const ::MGDoFHandler<1> &,
+ ConstraintMatrix &)
+ {
+ // nothing to do for regular
+ // dof handlers in 1d
+ }
+
+
+
+ static
+ void
+ make_oldstyle_hanging_node_constraints (const ::MGDoFHandler<1> &,
+ ConstraintMatrix &,
+ ::internal::int2type<1>)
+ {
+ // nothing to do for regular
+ // dof handlers in 1d
+ }
+
+
+ static
+ void
+ make_hp_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+ ConstraintMatrix &/*constraints*/)
+ {
+ // we may have to compute
+ // constraints for
+ // vertices. gotta think about
+ // that a bit more
+//TODO[WB]: think about what to do here...
+ }
+
+
+
static
void
- do_make_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
- ConstraintMatrix &/*constraints*/,
- ::internal::int2type<1>)
+ make_oldstyle_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+ ConstraintMatrix &/*constraints*/,
+ ::internal::int2type<1>)
{
// we may have to compute
// constraints for
#endif
+
#if deal_II_dimension == 2
template <class DH>
static
void
- do_make_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints,
- ::internal::int2type<2>)
+ make_oldstyle_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints,
+ ::internal::int2type<2>)
{
- // Decide whether to use the
- // new or old make_hanging_node_constraints
- // function. If all the FiniteElement
- // or all elements in a FECollection support
- // the new face constraint matrix, the
- // new code will be used.
- // Otherwise, the old implementation is used
- // for the moment.
- if (dof_handler.get_fe().hp_constraints_are_implemented ())
- {
- do_make_hp_hanging_node_constraints (dof_handler,
- constraints);
- return;
- }
-
const unsigned int dim = 2;
std::vector<unsigned int> dofs_on_mother;
ExcInternalError());
}
}
+#endif
+
+
+#if deal_II_dimension == 3
+ template <class DH>
+ static
+ void
+ make_oldstyle_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints,
+ ::internal::int2type<3>)
+ {
+ const unsigned int dim = 3;
+
+ 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);
+
+ // 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());
+
+ 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);
+
+ // 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);
+
+ // 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++]
+ = 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);
+
+ // 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);
+
+ // 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
-#if deal_II_dimension > 1
+
template <class DH>
static
void
- do_make_hp_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints)
+ make_hp_hanging_node_constraints (const DH &dof_handler,
+ ConstraintMatrix &constraints)
{
const unsigned int dim = DH::dimension;
}
}
}
-#endif
-
-
-
-#if deal_II_dimension == 3
- template <class DH>
- static
- void
- do_make_hanging_node_constraints (const DH &dof_handler,
- ConstraintMatrix &constraints,
- ::internal::int2type<3>)
- {
- // Decide whether to use the
- // new or old make_hanging_node_constraints
- // function. If all the FiniteElement
- // or all elements in a FECollection support
- // the new face constraint matrix, the
- // new code will be used.
- // Otherwise, the old implementation is used
- // for the moment.
- if (dof_handler.get_fe().hp_constraints_are_implemented ())
- {
- do_make_hp_hanging_node_constraints (dof_handler,
- constraints);
- return;
- }
-
- const unsigned int dim = 3;
-
- 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);
-
- // 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());
-
- 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);
-
- // 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);
-
- // 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++]
- = 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);
-
- // 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);
-
- // 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
-
}
}
DoFTools::make_hanging_node_constraints (const DH &dof_handler,
ConstraintMatrix &constraints)
{
- // distribute this one template
- // function to the appropriate
- // functions in 1d, 2d, and 3d
- internal::DoFTools
- ::do_make_hanging_node_constraints (dof_handler,
- constraints,
- ::internal::int2type<DH::dimension>());
+ // Decide whether to use the
+ // new or old make_hanging_node_constraints
+ // function. If all the FiniteElement
+ // or all elements in a FECollection support
+ // the new face constraint matrix, the
+ // new code will be used.
+ // Otherwise, the old implementation is used
+ // for the moment.
+ if (dof_handler.get_fe().hp_constraints_are_implemented ())
+ internal::DoFTools::
+ make_hp_hanging_node_constraints (dof_handler,
+ constraints);
+ else
+ internal::DoFTools::
+ make_oldstyle_hanging_node_constraints (dof_handler,
+ constraints,
+ ::internal::int2type<DH::dimension>());
}