From: bangerth Date: Wed, 20 Sep 2006 04:58:48 +0000 (+0000) Subject: Cache a few more things. Clean up a bit. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d6b1afa16760212793ebe6d2abfe76dbe8c4d1e;p=dealii-svn.git Cache a few more things. Clean up a bit. git-svn-id: https://svn.dealii.org/trunk@13933 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 091e7719ea..6a377b4dbb 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1658,6 +1658,78 @@ namespace internal { return 1; } + + + /** + * For a given face belonging + * to an active cell that + * borders to a more refined + * cell, return the fe_index of + * the most dominating finite + * element used on any of the + * face's subfaces. + */ + template + unsigned int + get_most_dominating_subface_fe_index (const face_iterator &face) + { + unsigned int dominating_subface_no = 0; + for (; dominating_subface_no::subfaces_per_face; + ++dominating_subface_no) + { + // each of the subfaces + // can have only a single + // fe_index associated + // with them, since there + // is no cell on the + // other side + Assert (face->child(dominating_subface_no) + ->n_active_fe_indices() + == 1, + ExcInternalError()); + + const FiniteElement & + this_subface_fe = (face->child(dominating_subface_no) + ->get_fe (face->child(dominating_subface_no) + ->nth_active_fe_index(0))); + + FiniteElementDomination::Domination + domination = FiniteElementDomination::either_element_can_dominate; + for (unsigned int sf=0; sf::subfaces_per_face; ++sf) + if (sf != dominating_subface_no) + { + const FiniteElement & + that_subface_fe = (face->child(sf) + ->get_fe (face->child(sf) + ->nth_active_fe_index(0))); + + domination = domination & + this_subface_fe.compare_for_face_domination(that_subface_fe); + } + + // see if the element + // on this subface is + // able to dominate + // the ones on all + // other subfaces, + // and if so take it + if ((domination == FiniteElementDomination::this_element_dominates) + || + (domination == FiniteElementDomination::either_element_can_dominate)) + break; + } + + // check that we have + // found one such subface + Assert (dominating_subface_no != GeometryInfo::subfaces_per_face, + ExcNotImplemented()); + + // return the finite element + // index used on it. note + // that only a single fe can + // be active on such subfaces + return face->child (dominating_subface_no)->nth_active_fe_index(0); + } @@ -2176,13 +2248,19 @@ namespace internal make_hp_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints) { -//TODO[WB]: we should call get_{sub,}face_interpolation_matrix only once per pair of finite elements -- since they return only stuff that counts in local dof numbers, their results are always the same. we should then just store the result in a lazy table -//TODO[WB]: we should short-circuit some of the more expensive things dealing with domination of finite elements for non-hp DoFHandlers. +//TODO[WB]: in the complex case, we should also store the inverses of the matrices in a cache table const unsigned int dim = DH::dimension; std::vector dofs_on_mother; std::vector dofs_on_children; + // a matrix to be used for + // constraints below. declared + // here and simply resized down + // below to avoid permanent + // re-allocation of memory + FullMatrix constraint_matrix; + // caches for the face and // subface interpolation // matrices between different @@ -2422,75 +2500,33 @@ namespace internal Assert (DoFHandlerSupportsDifferentP::value == true, ExcInternalError()); - // we first have to find - // one of the children - // for which the finite - // element is able to - // generate a space that - // all the other ones can - // be constrained to - unsigned int dominating_subface_no = 0; - for (; dominating_subface_no::subfaces_per_face; - ++dominating_subface_no) - { - // each of the - // subfaces can have - // only a single - // fe_index - // associated with - // them, since there - // is no cell on the - // other side - Assert (cell->face(face)->child(dominating_subface_no) - ->n_active_fe_indices() - == 1, - ExcInternalError()); - - const FiniteElement & - this_subface_fe = (cell->face(face)->child(dominating_subface_no) - ->get_fe (cell->face(face)->child(dominating_subface_no) - ->nth_active_fe_index(0))); - - FiniteElementDomination::Domination - domination = FiniteElementDomination::either_element_can_dominate; - for (unsigned int sf=0; sf::subfaces_per_face; ++sf) - if (sf != dominating_subface_no) - { - const FiniteElement & - that_subface_fe = (cell->face(face)->child(sf) - ->get_fe (cell->face(face)->child(sf) - ->nth_active_fe_index(0))); - - domination = domination & - this_subface_fe.compare_for_face_domination(that_subface_fe); - } - - // see if the element - // on this subface is - // able to dominate - // the ones on all - // other subfaces, - // and if so take it - if ((domination == FiniteElementDomination::this_element_dominates) - || - (domination == FiniteElementDomination::either_element_can_dominate)) - break; - } - - // check that we have - // found one such subface - Assert (dominating_subface_no != GeometryInfo::subfaces_per_face, - ExcNotImplemented()); - - const unsigned int dominating_fe_index - = cell->face(face)->child (dominating_subface_no)->nth_active_fe_index(0); - + // we first have + // to find the + // finite element + // that is able + // to generate a + // space that all + // the other ones + // can be + // constrained to + const unsigned int dominating_fe_index + = get_most_dominating_subface_fe_index (cell->face(face)); + const FiniteElement &dominating_fe - = cell->face(face)->child (dominating_subface_no)->get_fe (dominating_fe_index); + = dof_handler.get_fe()[dominating_fe_index]; - // check also that it is - // able to constrain the - // mother face + // check also + // that it is + // able to + // constrain the + // mother + // face. it + // should be, or + // we wouldn't + // have gotten + // into the + // branch for the + // 'complex' case Assert ((dominating_fe.compare_for_face_domination (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) == FiniteElementDomination::this_element_dominates) @@ -2508,10 +2544,16 @@ namespace internal Assert (dominating_fe.dofs_per_face <= cell->get_fe().dofs_per_face, ExcInternalError()); - FullMatrix restrict_mother_to_virtual (cell->get_fe().dofs_per_face, - dominating_fe.dofs_per_face); - dominating_fe.get_face_interpolation_matrix (cell->get_fe(), - restrict_mother_to_virtual); + + assert_existence_of_face_matrix + (dominating_fe, + cell->get_fe(), + face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]); + + const FullMatrix &restrict_mother_to_virtual + = *(face_interpolation_matrices + [dominating_fe_index][cell->active_fe_index()]); // split this matrix into // master and slave @@ -2536,11 +2578,11 @@ namespace internal restrict_mother_to_virtual_slave(i,j) = restrict_mother_to_virtual(i+dominating_fe.dofs_per_face,j); - FullMatrix face_constraints (cell->get_fe().dofs_per_face - - dominating_fe.dofs_per_face, - dominating_fe.dofs_per_face); + constraint_matrix.reinit (cell->get_fe().dofs_per_face - + dominating_fe.dofs_per_face, + dominating_fe.dofs_per_face); restrict_mother_to_virtual_slave - .mmult (face_constraints, + .mmult (constraint_matrix, restrict_mother_to_virtual_master_inv); std::vector face_dofs (cell->get_fe().dofs_per_face); @@ -2557,14 +2599,6 @@ namespace internal for (unsigned int i=0; iget_fe().dofs_per_face - dominating_fe.dofs_per_face; ++i) slave_dofs[i] = face_dofs[i + dominating_fe.dofs_per_face]; - - for (unsigned int i=0; i