From 725d6be8b0696afd71753f50d639d0290752f59f Mon Sep 17 00:00:00 2001 From: bangerth Date: Sun, 13 Aug 2006 00:51:14 +0000 Subject: [PATCH] Cleanup some of the code; remove one level of indirection and re-sort functions to make reading through them a bit easier git-svn-id: https://svn.dealii.org/trunk@13690 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/dofs/dof_tools.cc | 575 ++++++++++++----------- 1 file changed, 300 insertions(+), 275 deletions(-) diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 8147d1a7d9..cb60456dff 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1625,20 +1625,67 @@ namespace internal #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 @@ -1649,29 +1696,15 @@ namespace internal #endif + #if deal_II_dimension == 2 template 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 dofs_on_mother; @@ -1807,14 +1840,237 @@ namespace internal ExcInternalError()); } } +#endif + + +#if deal_II_dimension == 3 + template + static + void + make_oldstyle_hanging_node_constraints (const DH &dof_handler, + ConstraintMatrix &constraints, + ::internal::int2type<3>) + { + const unsigned int dim = 3; + + std::vector dofs_on_mother; + std::vector 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::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::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::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 &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; dofchild(0)->line(1)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(2)->line(1)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(0)->line(3)->dof_index(dof, fe_index); + for (unsigned int dof=0; dofchild(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 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; @@ -2142,247 +2398,6 @@ namespace internal } } } -#endif - - - -#if deal_II_dimension == 3 - template - 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 dofs_on_mother; - std::vector 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::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::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::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 &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; dofchild(0)->line(1)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(2)->line(1)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(0)->line(3)->dof_index(dof, fe_index); - for (unsigned int dof=0; dofchild(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 - } } @@ -2393,13 +2408,23 @@ void 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()); + // 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()); } -- 2.39.5