From ee3d7e892a65e8184112d3444b0ba840d3c3ca45 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 1 Sep 2006 20:29:15 +0000 Subject: [PATCH] Move the Domination enum out of the FiniteElementBase class so that we can actually define operator| acting on it. Use this stuff in the computation of constraints. git-svn-id: https://svn.dealii.org/trunk@13794 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 2 +- deal.II/deal.II/include/fe/fe_base.h | 229 ++++-- deal.II/deal.II/include/fe/fe_dgp.h | 2 +- deal.II/deal.II/include/fe/fe_dgp_monomial.h | 2 +- .../deal.II/include/fe/fe_dgp_nonparametric.h | 2 +- deal.II/deal.II/include/fe/fe_dgq.h | 2 +- deal.II/deal.II/include/fe/fe_q.h | 2 +- deal.II/deal.II/include/fe/fe_system.h | 2 +- deal.II/deal.II/source/dofs/dof_tools.cc | 708 ++++++++++-------- deal.II/deal.II/source/fe/fe.cc | 4 +- deal.II/deal.II/source/fe/fe_dgp.cc | 8 +- deal.II/deal.II/source/fe/fe_dgp_monomial.cc | 8 +- .../deal.II/source/fe/fe_dgp_nonparametric.cc | 8 +- deal.II/deal.II/source/fe/fe_dgq.cc | 8 +- deal.II/deal.II/source/fe/fe_q.cc | 10 +- deal.II/deal.II/source/fe/fe_system.cc | 95 +-- 16 files changed, 576 insertions(+), 516 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 286334e481..399646cac7 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -1072,7 +1072,7 @@ class FiniteElement : public Subscriptor, * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; //@} diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 066c71fd52..c03c3a61b2 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -30,6 +30,120 @@ template class FESystem; + +/** + * A namespace solely for the purpose of defining the Domination enum as well + * as associated operators. + */ +namespace FiniteElementDomination +{ + /** + * An enum that describes the + * outcome of comparing two elements for + * mutual domination. If one element + * dominates another, then the + * restriction of the space described by + * the dominated element to a face of the + * cell is strictly larger than that of + * the dominating element. For example, + * in 2-d Q(2) elements dominate Q(4) + * elements, because the traces of Q(4) + * elements are quartic polynomials which + * is a space strictly larger than the + * quadratic polynomials (the restriction + * of the Q(2) element). In general, Q(k) + * dominates Q(k') if $k\le k'$. + * + * This enum is used in the + * FiniteElement::compare_fe_for_domination() + * function that is used in the context + * of hp finite element methods when + * determining what to do at faces where + * two different finite elements meet + * (see the hp paper for a more detailed + * description of the following). In that + * case, the degrees of freedom of one + * side need to be constrained to those + * on the other side. The determination + * which side is which is based on the + * outcome of a comparison for mutual + * domination: the dominated side is + * constrained to the dominating one. + * + * Note that there are situations where + * neither side dominates. The hp paper + * lists two case, with the simpler one + * being that a $Q_2\times Q_1$ + * vector-valued element (i.e. a + * FESystem(FE_Q(2),1,FE_Q(1),1)) + * meets a $Q_1\times Q_2$ element: here, + * for each of the two vector-components, + * we can define a domination + * relationship, but it is different for + * the two components. + * + * It is clear that the concept of + * domination doesn't matter for + * discontinuous elements. However, + * discontinuous elements may be part of + * vector-valued elements and may + * therefore be compared against each + * other for domination. They should + * return + * either_element_can_dominate + * in that case. Likewise, when comparing + * two identical finite elements, they + * should return this code; the reason is + * that we can not decide which element + * will dominate at the time we look at + * the first component of, for example, + * two $Q_2\times Q_1$ and $Q_2\times + * Q_2$ elements, and have to keep our + * options open until we get to the + * second base element. + */ + enum Domination + { + this_element_dominates, + other_element_dominates, + neither_element_dominates, + either_element_can_dominate + }; + + + /** + * A generalization of the binary + * or operator to a comparison + * relationship. The way this works is + * pretty much as when you would want to + * define a comparison relationship for + * vectors: either all elements of the + * first vector are smaller, equal, or + * larger than those of the second vector, + * or some are and some are not. + * + * This operator is pretty much the same: + * if both arguments are + * this_element_dominates or + * other_element_dominates, + * then the returned value is that + * value. On the other hand, if one of the + * values is + * either_element_can_dominate, + * then the returned value is that of the + * other argument. If either argument is + * neither_element_dominates, + * or if the two arguments are + * this_element_dominates and + * other_element_dominates, + * then the returned value is + * neither_element_dominates. + */ + Domination operator | (const Domination d1, + const Domination d2); +} + + /** * Dimension independent data for finite elements. See the derived * class FiniteElement class for information on its use. All @@ -155,81 +269,6 @@ class FiniteElementData */ H2 = 0x0e }; - - /** - * An enum that describes the - * outcome of comparing two elements for - * mutual domination. If one element - * dominates another, then the - * restriction of the space described by - * the dominated element to a face of the - * cell is strictly larger than that of - * the dominating element. For example, - * in 2-d Q(2) elements dominate Q(4) - * elements, because the traces of Q(4) - * elements are quartic polynomials which - * is a space strictly larger than the - * quadratic polynomials (the restriction - * of the Q(2) element). In general, Q(k) - * dominates Q(k') if $k\le k'$. - * - * This enum is used in the - * FiniteElement::compare_fe_for_domination() - * function that is used in the context - * of hp finite element methods when - * determining what to do at faces where - * two different finite elements meet - * (see the hp paper for a more detailed - * description of the following). In that - * case, the degrees of freedom of one - * side need to be constrained to those - * on the other side. The determination - * which side is which is based on the - * outcome of a comparison for mutual - * domination: the dominated side is - * constrained to the dominating one. - * - * Note that there are situations where - * neither side dominates. The hp paper - * lists two case, with the simpler one - * being that a $Q_2\times Q_1$ - * vector-valued element (i.e. a - * FESystem(FE_Q(2),1,FE_Q(1),1)) - * meets a $Q_1\times Q_2$ element: here, - * for each of the two vector-components, - * we can define a domination - * relationship, but it is different for - * the two components. - * - * It is clear that the concept of - * domination doesn't matter for - * discontinuous elements. However, - * discontinuous elements may be part of - * vector-valued elements and may - * therefore be compared against each - * other for domination. They should - * return - * either_element_can_dominate - * in that case. Likewise, when comparing - * two identical finite elements, they - * should return this code; the reason is - * that we can not decide which element - * will dominate at the time we look at - * the first component of, for example, - * two $Q_2\times Q_1$ and $Q_2\times - * Q_2$ elements, and have to keep our - * options open until we get to the - * second base element. - */ - enum Domination - { - this_element_dominates, - other_element_dominates, - neither_element_dominates, - either_element_can_dominate - }; - - /** * Number of degrees of freedom on @@ -550,6 +589,46 @@ class FiniteElementData // --------- inline and template functions --------------- + + +namespace FiniteElementDomination +{ + inline + Domination operator | (const Domination d1, + const Domination d2) + { + switch (d1) + { + case this_element_dominates: + if ((d2 == this_element_dominates) || + (d2 == either_element_can_dominate)) + return this_element_dominates; + else + return neither_element_dominates; + + case other_element_dominates: + if ((d2 == other_element_dominates) || + (d2 == either_element_can_dominate)) + return other_element_dominates; + else + return neither_element_dominates; + + case neither_element_dominates: + return neither_element_dominates; + + case either_element_can_dominate: + return d2; + + default: + // shouldn't get here + Assert (false, ExcInternalError()); + } + + return neither_element_dominates; + } +} + + template inline unsigned int diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h index 97a9dbd340..d790e9689d 100644 --- a/deal.II/deal.II/include/fe/fe_dgp.h +++ b/deal.II/deal.II/include/fe/fe_dgp.h @@ -127,7 +127,7 @@ class FE_DGP : public FE_Poly,dim> * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; /** diff --git a/deal.II/deal.II/include/fe/fe_dgp_monomial.h b/deal.II/deal.II/include/fe/fe_dgp_monomial.h index 7d2a27deab..e27bb6826c 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_monomial.h +++ b/deal.II/deal.II/include/fe/fe_dgp_monomial.h @@ -128,7 +128,7 @@ class FE_DGPMonomial : public FE_Poly,dim> * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; /** diff --git a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h index 9d6031478f..01a79fb166 100644 --- a/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h +++ b/deal.II/deal.II/include/fe/fe_dgp_nonparametric.h @@ -283,7 +283,7 @@ class FE_DGPNonparametric : public FiniteElement * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; /** diff --git a/deal.II/deal.II/include/fe/fe_dgq.h b/deal.II/deal.II/include/fe/fe_dgq.h index a9ecfb6dff..2eeca8b1ed 100644 --- a/deal.II/deal.II/include/fe/fe_dgq.h +++ b/deal.II/deal.II/include/fe/fe_dgq.h @@ -199,7 +199,7 @@ class FE_DGQ : public FE_Poly,dim> * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; /** diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index c90ee992a3..dd7f60e28f 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -433,7 +433,7 @@ class FE_Q : public FE_Poly,dim> * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; //@} diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index cb0eb81ac8..e807507044 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -535,7 +535,7 @@ class FESystem : public FiniteElement * particular the hp paper. */ virtual - typename FiniteElementData::Domination + FiniteElementDomination::Domination compare_for_domination (const FiniteElement &fe_other) const; //@} diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 0eef3d0407..b59925e8d5 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -2078,8 +2078,8 @@ namespace internal FullMatrix face_constraints; - // loop over all lines; only on - // lines there can be constraints. + // loop over all faces; only on + // face there can be constraints. // We do so by looping over all // active cells and checking // whether any of the faces are @@ -2128,335 +2128,386 @@ namespace internal Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, ExcInternalError()); - -//TODO: The proper way would be to ask the elements themselves which -//among them should be the master element. Have a poor-man's -//implementation here that simply decides this based on the number of -//DoFs per face - - // Store minimum degree element. - // For FE_Q it is the one with the - // lowest number of DoFs on the face. - unsigned int min_dofs_per_face = cell->get_fe ().dofs_per_face; - unsigned int min_degree_subface = 0; - bool mother_face_is_master = true; - + // first find out whether we + // can constrain each of the + // subfaces to the mother + // face. in the lingo of the hp + // paper, this would be the + // simple case + FiniteElementDomination::Domination + mother_face_dominates = FiniteElementDomination::either_element_can_dominate; for (unsigned int c=0; c::subfaces_per_face; ++c) + mother_face_dominates = mother_face_dominates | + (cell->get_fe().compare_for_domination + (cell->neighbor_child_on_subface (face, c)->get_fe())); + + switch (mother_face_dominates) { - typename DH::active_cell_iterator - neighbor_child - = cell->neighbor_child_on_subface (face, c); - - // Check if the element on one - // of the subfaces has a lower - // polynomial degree than the - // one of the other elements. - if (neighbor_child->get_fe ().dofs_per_face < min_dofs_per_face) - { - min_dofs_per_face = neighbor_child->get_fe ().dofs_per_face; - min_degree_subface = c; - mother_face_is_master = false; - } - } - // Case 1: The coarse element has - // the lowest polynomial degree. - // Therefore it will play the role - // of the master elements, to which - // the other elements will be constrained. - if (mother_face_is_master == true) - { - const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face; - dofs_on_mother.resize (n_dofs_on_mother); + case FiniteElementDomination::this_element_dominates: + { + // Case 1 (the simple case): + // The coarse element dominates + // the elements on the subfaces + const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face; + dofs_on_mother.resize (n_dofs_on_mother); - cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ()); + cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ()); - // Now create constraint matrix for - // the subfaces and assemble it. - for (unsigned int c=0; c::subfaces_per_face; ++c) - { - typename DH::active_cell_iterator neighbor_child - = cell->neighbor_child_on_subface (face, c); - - const unsigned int n_dofs_on_children - = neighbor_child->get_fe().dofs_per_face; - dofs_on_children.resize (n_dofs_on_children); - - const unsigned int subface_fe_index - = neighbor_child->active_fe_index(); - - // some sanity checks - // -- particularly - // useful if you start - // to think about faces - // with - // face_orientation==false - // and whether we - // really really have - // the right face... - Assert (neighbor_child->n_active_fe_indices() == 1, - ExcInternalError()); - Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, - ExcInternalError()); - Assert (cell->face(face)->child(c)->fe_index_is_active(subface_fe_index) - == true, - ExcInternalError()); + // Now create constraint matrix for + // the subfaces and assemble it. + for (unsigned int c=0; c::subfaces_per_face; ++c) + { + typename DH::active_cell_iterator neighbor_child + = cell->neighbor_child_on_subface (face, c); + + const unsigned int n_dofs_on_children + = neighbor_child->get_fe().dofs_per_face; + dofs_on_children.resize (n_dofs_on_children); + + const unsigned int subface_fe_index + = neighbor_child->active_fe_index(); + + // some sanity checks + // -- particularly + // useful if you start + // to think about faces + // with + // face_orientation==false + // and whether we + // really really have + // the right face... + Assert (neighbor_child->n_active_fe_indices() == 1, + ExcInternalError()); + Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, + ExcInternalError()); + Assert (cell->face(face)->child(c)->fe_index_is_active(subface_fe_index) + == true, + ExcInternalError()); - // Same procedure as for the - // mother cell. Extract the face - // DoFs from the cell DoFs. - cell->face(face)->child(c) - ->get_dof_indices (dofs_on_children, - subface_fe_index); + // Same procedure as for the + // mother cell. Extract the face + // DoFs from the cell DoFs. + cell->face(face)->child(c) + ->get_dof_indices (dofs_on_children, + subface_fe_index); - // Now create the - // element constraint - // for this subface. - // - // As a side remark, - // one may wonder the - // following: - // neighbor_child is - // clearly computed - // correctly, - // i.e. taking into - // account - // face_orientation - // (just look at the - // implementation of - // that - // function). however, - // we don't care about - // this here, when we - // ask for - // subface_interpolation - // on subface c. the - // question rather is: - // do we have to - // translate 'c' here - // as well? - // - // the answer is in - // fact 'no'. if one - // does that, results - // are wrong: - // constraints are - // added twice for the - // same pair of nodes - // but with differing - // weights. in - // addition, one can - // look at the - // deal.II/project_*_03 - // tests that look at - // exactly this case: - // there, we have a - // mesh with at least - // one - // face_orientation==false - // and hanging nodes, - // and the results of - // those tests show - // that the result of - // projection verifies - // the approximation - // properties of a - // finite element onto - // that mesh - face_constraints.reinit (n_dofs_on_mother, - n_dofs_on_children); - cell->get_fe() - .get_subface_interpolation_matrix (cell->get_dof_handler() - .get_fe()[subface_fe_index], - c, face_constraints); - - // Add constraints to global constraint - // matrix. - filter_constraints (dofs_on_mother, - dofs_on_children, - face_constraints, - constraints); - } - } - else - // Case 2: One of the finer elements - // has the lowest polynomial degree. - // First the coarse element will be - // constrained to that element. After - // that the other fine elements will - // be constrained to the coarse element. - { - // Assert (false, ExcNotImplemented ()); - - typename DH::active_cell_iterator neighbor_child - = cell->neighbor_child_on_subface (face, min_degree_subface); - const unsigned int n_dofs_on_children = neighbor_child->get_fe().dofs_per_face; - dofs_on_children.resize (n_dofs_on_children); - + // Now create the + // element constraint + // for this subface. + // + // As a side remark, + // one may wonder the + // following: + // neighbor_child is + // clearly computed + // correctly, + // i.e. taking into + // account + // face_orientation + // (just look at the + // implementation of + // that + // function). however, + // we don't care about + // this here, when we + // ask for + // subface_interpolation + // on subface c. the + // question rather is: + // do we have to + // translate 'c' here + // as well? + // + // the answer is in + // fact 'no'. if one + // does that, results + // are wrong: + // constraints are + // added twice for the + // same pair of nodes + // but with differing + // weights. in + // addition, one can + // look at the + // deal.II/project_*_03 + // tests that look at + // exactly this case: + // there, we have a + // mesh with at least + // one + // face_orientation==false + // and hanging nodes, + // and the results of + // those tests show + // that the result of + // projection verifies + // the approximation + // properties of a + // finite element onto + // that mesh + face_constraints.reinit (n_dofs_on_mother, + n_dofs_on_children); + cell->get_fe() + .get_subface_interpolation_matrix (cell->get_dof_handler() + .get_fe()[subface_fe_index], + c, face_constraints); + + // Add constraints to global constraint + // matrix. + filter_constraints (dofs_on_mother, + dofs_on_children, + face_constraints, + constraints); + } + + break; + } - // Get DoFs on child cell with - // lowest polynomial degree. - // All other DoFs will be constrained - // to the DoFs of this face. - const unsigned int subface_fe_index - = neighbor_child->active_fe_index(); + case FiniteElementDomination::other_element_dominates: + case FiniteElementDomination::neither_element_dominates: + { + // Case 2 (the "complex" + // case): at least one + // (the neither_... case) + // of the finer elements + // or all of them (the + // other_... case) is + // dominating. First the + // coarse element will be + // constrained to that + // element. After that + // the other fine + // elements will be + // constrained to the + // coarse element. + + // 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) + { + FiniteElementDomination::Domination + domination = FiniteElementDomination::either_element_can_dominate; + for (unsigned int sf=0; sf::subfaces_per_face; ++sf) + if (sf != dominating_subface_no) + domination = domination | + cell->neighbor_child_on_subface (face, dominating_subface_no) + ->get_fe().compare_for_domination + (cell->neighbor_child_on_subface (face, sf)->get_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 typename DH::active_cell_iterator neighbor_child + = cell->neighbor_child_on_subface (face, dominating_subface_no); + const unsigned int n_dofs_on_children = neighbor_child->get_fe().dofs_per_face; + dofs_on_children.resize (n_dofs_on_children); + + + // Get DoFs on child cell with + // lowest polynomial degree. + // All other DoFs will be constrained + // to the DoFs of this face. + const unsigned int subface_fe_index + = neighbor_child->active_fe_index(); - // Same procedure as for the - // mother cell. Extract the face - // DoFs from the cell DoFs. - cell->face(face)->child(min_degree_subface) - ->get_dof_indices (dofs_on_children, - subface_fe_index); - - - // The idea is to introduce - // a "virtual" intermediate coarse - // level face with the lowest - // polynomial degree. Then it is - // easy to constrain each of the - // connected faces to this intermediate - // coarse level face. As the DoFs on - // this intermediate coarse level face - // do not exist, they have to determined - // through the inverse of the constraint matrix - // from the lowest order subface to - // this intermediate coarse level face. - // - // Considering the following case: - // +---+----+ - // | | Q3 | - // |Q3 +----+ - // | | Q2 | - // +---+----+ - // - // The intermediate layer would be - // of order 2: - // +------+ * +---------+ - // + | | | F_1, Q3 | - // +Q3, C | * +---------+ - // + | | | F_2, Q2 | - // +------+ * +---------+ - // - // In this case, there are 3 DoFs on the - // intermediate layer. Assuming for the - // moment that these do exist, all DoFs - // on the connected faces can be - // expressed in terms of these DoFs. We - // have: - // C = A_1 * I - // F_1 = A_2 * I - // F_2 = A_3 * I - // where C, F_1, F_2 denote the DoFs - // on the faces of the elements and - // I denotes the DoFs on the intermediate - // face. A_1 to A_3 denote the corresponding - // face or subface interpolation matrices, - // describing the DoFs on one of the faces - // in terms of the DoFs on the intermediate - // layer. - // - // As the DoFs in I are only "virtual" - // they have to be expressed in terms - // of existing DoFs. In this case only - // A_3 is invertible. Therefore all - // other DoFs have to be constrained - // to the DoFs in F_2. - // This leads to - // I = A_3^-1 F_2 - // and - // C = A_1 * A_3^-1 F_2 - // F_1 = A_2 * A_3^-1 F_2 - // - // Therefore the constraint matrices - // in this case are: - // A_1 * A_3^-1 - // A_2 * A_3^-1 - // In 3D and for other configurations, - // the basic scheme is completely identical. + // Same procedure as for the + // mother cell. Extract the face + // DoFs from the cell DoFs. + cell->face(face)->child(dominating_subface_no) + ->get_dof_indices (dofs_on_children, + subface_fe_index); + + + // The idea is to introduce + // a "virtual" intermediate coarse + // level face with the lowest + // polynomial degree. Then it is + // easy to constrain each of the + // connected faces to this intermediate + // coarse level face. As the DoFs on + // this intermediate coarse level face + // do not exist, they have to determined + // through the inverse of the constraint matrix + // from the lowest order subface to + // this intermediate coarse level face. + // + // Considering the following case: + // +---+----+ + // | | Q3 | + // |Q3 +----+ + // | | Q2 | + // +---+----+ + // + // The intermediate layer would be + // of order 2: + // +------+ * +---------+ + // + | | | F_1, Q3 | + // +Q3, C | * +---------+ + // + | | | F_2, Q2 | + // +------+ * +---------+ + // + // In this case, there are 3 DoFs on the + // intermediate layer. Assuming for the + // moment that these do exist, all DoFs + // on the connected faces can be + // expressed in terms of these DoFs. We + // have: + // C = A_1 * I + // F_1 = A_2 * I + // F_2 = A_3 * I + // where C, F_1, F_2 denote the DoFs + // on the faces of the elements and + // I denotes the DoFs on the intermediate + // face. A_1 to A_3 denote the corresponding + // face or subface interpolation matrices, + // describing the DoFs on one of the faces + // in terms of the DoFs on the intermediate + // layer. + // + // As the DoFs in I are only "virtual" + // they have to be expressed in terms + // of existing DoFs. In this case only + // A_3 is invertible. Therefore all + // other DoFs have to be constrained + // to the DoFs in F_2. + // This leads to + // I = A_3^-1 F_2 + // and + // C = A_1 * A_3^-1 F_2 + // F_1 = A_2 * A_3^-1 F_2 + // + // Therefore the constraint matrices + // in this case are: + // A_1 * A_3^-1 + // A_2 * A_3^-1 + // In 3D and for other configurations, + // the basic scheme is completely identical. - // Now create the element - // constraint for this subface. - FullMatrix fc_sface_ipol (n_dofs_on_children, - n_dofs_on_children); - FullMatrix fc_ipol_sface (n_dofs_on_children, - n_dofs_on_children); - neighbor_child->get_fe().get_subface_interpolation_matrix (neighbor_child->get_fe (), - min_degree_subface, - fc_sface_ipol); - // Invert it, to get a mapping from the DoFs of the - // "Master-subface" to the intermediate layer. - fc_ipol_sface.invert (fc_sface_ipol); - - // Create constraint matrix for the mother face. - const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face; - dofs_on_mother.resize (n_dofs_on_mother); - cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ()); - - FullMatrix fc_mother_ipol (n_dofs_on_children, - n_dofs_on_mother); - FullMatrix fc_mother_sface (n_dofs_on_children, - n_dofs_on_mother); - neighbor_child->get_fe ().get_face_interpolation_matrix (cell->get_fe(), - fc_mother_ipol); - fc_ipol_sface.mmult (fc_mother_sface, fc_mother_ipol); - - // Add constraints to global constraint - // matrix. - filter_constraints (dofs_on_children, - dofs_on_mother, - fc_mother_sface, - constraints); - - // Now create constraint matrices for - // the subfaces and assemble them - for (unsigned int c=0; c::subfaces_per_face; ++c) - { - // As the "Master-subface" does not need constraints, skip it. - if (c != min_degree_subface) - { - typename DH::active_cell_iterator neighbor_child_slave - = cell->neighbor_child_on_subface (face, c); - const unsigned int n_dofs_on_mother = neighbor_child_slave->get_fe().dofs_per_face; - dofs_on_mother.resize (n_dofs_on_mother); - - // Find face number on the finer - // neighboring cell, which is - // shared the face with the - // face of the coarser cell. - const unsigned int neighbor2= - cell->neighbor_of_neighbor(face); - Assert (neighbor_child_slave->face(neighbor2) == cell->face(face)->child(c), - ExcInternalError()); + // Now create the element + // constraint for this subface. + FullMatrix fc_sface_ipol (n_dofs_on_children, + n_dofs_on_children); + FullMatrix fc_ipol_sface (n_dofs_on_children, + n_dofs_on_children); + neighbor_child->get_fe().get_subface_interpolation_matrix (neighbor_child->get_fe (), + dominating_subface_no, + fc_sface_ipol); + // Invert it, to get a mapping from the DoFs of the + // "Master-subface" to the intermediate layer. + fc_ipol_sface.invert (fc_sface_ipol); + + // Create constraint matrix for the mother face. + const unsigned int n_dofs_on_mother = cell->get_fe().dofs_per_face; + dofs_on_mother.resize (n_dofs_on_mother); + cell->face(face)->get_dof_indices (dofs_on_mother, cell->active_fe_index ()); + + FullMatrix fc_mother_ipol (n_dofs_on_children, + n_dofs_on_mother); + FullMatrix fc_mother_sface (n_dofs_on_children, + n_dofs_on_mother); + neighbor_child->get_fe ().get_face_interpolation_matrix (cell->get_fe(), + fc_mother_ipol); + fc_ipol_sface.mmult (fc_mother_sface, fc_mother_ipol); + + // Add constraints to global constraint + // matrix. + filter_constraints (dofs_on_children, + dofs_on_mother, + fc_mother_sface, + constraints); + + // Now create constraint matrices for + // the subfaces and assemble them + for (unsigned int c=0; c::subfaces_per_face; ++c) + { + // As the "Master-subface" does not need constraints, skip it. + if (c != dominating_subface_no) + { + typename DH::active_cell_iterator neighbor_child_slave + = cell->neighbor_child_on_subface (face, c); + const unsigned int n_dofs_on_mother = neighbor_child_slave->get_fe().dofs_per_face; + dofs_on_mother.resize (n_dofs_on_mother); + + // Find face number on the finer + // neighboring cell, which is + // shared the face with the + // face of the coarser cell. + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face); + Assert (neighbor_child_slave->face(neighbor2) == cell->face(face)->child(c), + ExcInternalError()); - // Same procedure as for the - // mother cell. Extract the face - // DoFs from the cell DoFs. - const unsigned int subface_fe_index - = neighbor_child_slave->active_fe_index(); + // Same procedure as for the + // mother cell. Extract the face + // DoFs from the cell DoFs. + const unsigned int subface_fe_index + = neighbor_child_slave->active_fe_index(); - cell->face(face)->child(c) - ->get_dof_indices (dofs_on_mother, - subface_fe_index); + cell->face(face)->child(c) + ->get_dof_indices (dofs_on_mother, + subface_fe_index); - // Now create the element - // constraint for this subface. - FullMatrix fc_child_sface_ipol (n_dofs_on_children, - n_dofs_on_mother); - FullMatrix fc_child_sface_sface (n_dofs_on_children, - n_dofs_on_mother); - neighbor_child->get_fe ().get_subface_interpolation_matrix - (neighbor_child_slave->get_fe(), - c, fc_child_sface_ipol); - - fc_ipol_sface.mmult (fc_child_sface_sface, fc_child_sface_ipol); - - // Add constraints to global constraint - // matrix. - filter_constraints (dofs_on_children, - dofs_on_mother, - fc_child_sface_sface, - constraints); - } - } + // Now create the element + // constraint for this subface. + FullMatrix fc_child_sface_ipol (n_dofs_on_children, + n_dofs_on_mother); + FullMatrix fc_child_sface_sface (n_dofs_on_children, + n_dofs_on_mother); + neighbor_child->get_fe ().get_subface_interpolation_matrix + (neighbor_child_slave->get_fe(), + c, fc_child_sface_ipol); + + fc_ipol_sface.mmult (fc_child_sface_sface, fc_child_sface_ipol); + + // Add constraints to global constraint + // matrix. + filter_constraints (dofs_on_children, + dofs_on_mother, + fc_child_sface_sface, + constraints); + } + } + + break; + } + + case FiniteElementDomination::either_element_can_dominate: + { + // hm, it isn't quite + // clear what exactly we + // would have to do + // here. sit tight until + // someone trips over the + // following statement + // and see what exactly + // is going on + Assert (false, ExcNotImplemented()); + } + + default: + // we shouldn't get here + Assert (false, ExcInternalError()); } } else @@ -2493,7 +2544,7 @@ namespace internal // constrain switch (cell->get_fe().compare_for_domination (neighbor->get_fe ())) { - case FiniteElementData::this_element_dominates: + case FiniteElementDomination::this_element_dominates: { // Get DoFs on // dominating and @@ -2525,7 +2576,7 @@ namespace internal break; } - case FiniteElementData::other_element_dominates: + case FiniteElementDomination::other_element_dominates: { // we don't do anything // here since we will @@ -2538,20 +2589,27 @@ namespace internal break; } - case FiniteElementData::either_element_can_dominate: + case FiniteElementDomination::either_element_can_dominate: { // it appears as if // neither element has // any constraints on - // its neighbor + // its neighbor. break; } - case FiniteElementData::neither_element_dominates: + case FiniteElementDomination::neither_element_dominates: { // we don't presently // know what exactly to - // do here + // do here. it isn't quite + // clear what exactly we + // would have to do + // here. sit tight until + // someone trips over the + // following statement + // and see what exactly + // is going on Assert (false, ExcNotImplemented()); break; } diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index f5f620729a..3ca65daf19 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -519,12 +519,12 @@ hp_quad_dof_identities (const FiniteElement &) const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FiniteElement:: compare_for_domination (const FiniteElement &) const { Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index 3e99b0892d..7d264e6734 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -156,18 +156,18 @@ FE_DGP::hp_constraints_are_implemented () const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FE_DGP::compare_for_domination (const FiniteElement &fe_other) const { // check whether both are discontinuous // elements and both could dominate, see // the description of - // FiniteElementData::Domination + // FiniteElementDomination::Domination if (dynamic_cast*>(&fe_other) != 0) - return FiniteElementData::either_element_can_dominate; + return FiniteElementDomination::either_element_can_dominate; Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc index 347ea12afc..b0c3245d71 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc @@ -317,19 +317,19 @@ FE_DGPMonomial::hp_constraints_are_implemented () const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FE_DGPMonomial:: compare_for_domination (const FiniteElement &fe_other) const { // check whether both are discontinuous // elements and both could dominate, see // the description of - // FiniteElementData::Domination + // FiniteElementDomination::Domination if (dynamic_cast*>(&fe_other) != 0) - return FiniteElementData::either_element_can_dominate; + return FiniteElementDomination::either_element_can_dominate; Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc index 73579bdfa2..1d98ec84ad 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -486,19 +486,19 @@ FE_DGPNonparametric::hp_constraints_are_implemented () const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FE_DGPNonparametric:: compare_for_domination (const FiniteElement &fe_other) const { // check whether both are discontinuous // elements and both could dominate, see // the description of - // FiniteElementData::Domination + // FiniteElementDomination::Domination if (dynamic_cast*>(&fe_other) != 0) - return FiniteElementData::either_element_can_dominate; + return FiniteElementDomination::either_element_can_dominate; Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index f69c236bc3..8dd3b05ded 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -493,18 +493,18 @@ FE_DGQ::hp_constraints_are_implemented () const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FE_DGQ::compare_for_domination (const FiniteElement &fe_other) const { // check whether both are discontinuous // elements and both could dominate, see // the description of - // FiniteElementData::Domination + // FiniteElementDomination::Domination if (dynamic_cast*>(&fe_other) != 0) - return FiniteElementData::either_element_can_dominate; + return FiniteElementDomination::either_element_can_dominate; Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index d32ff982ac..bd87fa0636 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -679,7 +679,7 @@ hp_quad_dof_identities (const FiniteElement &fe_other) const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FE_Q:: compare_for_domination (const FiniteElement &fe_other) const { @@ -687,15 +687,15 @@ compare_for_domination (const FiniteElement &fe_other) const = dynamic_cast*>(&fe_other)) { if (this->degree < fe_q_other->degree) - return FiniteElementData::this_element_dominates; + return FiniteElementDomination::this_element_dominates; else if (this->degree < fe_q_other->degree) - return FiniteElementData::either_element_can_dominate; + return FiniteElementDomination::either_element_can_dominate; else - return FiniteElementData::other_element_dominates; + return FiniteElementDomination::other_element_dominates; } Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index d9c04abe0e..5a877700eb 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -2273,7 +2273,7 @@ FESystem::hp_quad_dof_identities (const FiniteElement &fe_other) const template -typename FiniteElementData::Domination +FiniteElementDomination::Domination FESystem:: compare_for_domination (const FiniteElement &fe_other) const { @@ -2288,8 +2288,8 @@ compare_for_domination (const FiniteElement &fe_other) const Assert (this->n_base_elements() == fe_sys_other->n_base_elements(), ExcNotImplemented()); - typename FiniteElementData::Domination - domination = FiniteElementData::either_element_can_dominate; + FiniteElementDomination::Domination + domination = FiniteElementDomination::either_element_can_dominate; // loop over all base elements and do // some sanity checks @@ -2303,101 +2303,24 @@ compare_for_domination (const FiniteElement &fe_other) const ExcNotImplemented()); // for this pair of base elements, - // check who dominates - const typename FiniteElementData::Domination - base_domination - = (this->base_element(b) - .compare_for_domination (fe_sys_other->base_element(b))); - - // now see what that means with - // regard to the previous state - switch (domination) - { - case FiniteElementData::either_element_can_dominate: - { - // we haven't made a decision - // yet, simply copy what this - // pair of bases have to say - domination = base_domination; - break; - } - - case FiniteElementData::this_element_dominates: - { - // the present element - // previously dominated. this - // will still be the case if - // either the present base - // dominates or if the two - // bases don't - // care. otherwise, there is - // a tie which we will not be - // able to escape from no - // matter what the other - // pairs of bases are going - // to say - switch (base_domination) - { - case FiniteElementData::either_element_can_dominate: - case FiniteElementData::this_element_dominates: - { - break; - } - - case FiniteElementData::other_element_dominates: - case FiniteElementData::neither_element_dominates: - { - return FiniteElementData::neither_element_dominates; - } - - default: - // shouldn't get - // here - Assert (false, ExcInternalError()); - } - break; - } - - case FiniteElementData::other_element_dominates: - { - // the opposite case - switch (base_domination) - { - case FiniteElementData::either_element_can_dominate: - case FiniteElementData::other_element_dominates: - { - break; - } - - case FiniteElementData::this_element_dominates: - case FiniteElementData::neither_element_dominates: - { - return FiniteElementData::neither_element_dominates; - } - - default: - Assert (false, ExcInternalError()); - } - break; - } - - default: - Assert (false, ExcInternalError()); - } + // check who dominates and combine + // with previous result + domination = domination | (this->base_element(b) + .compare_for_domination (fe_sys_other->base_element(b))); } // if we've gotten here, then we've // either found a winner or either // element is fine being dominated Assert (domination != - FiniteElementData::neither_element_dominates, + FiniteElementDomination::neither_element_dominates, ExcInternalError()); return domination; } Assert (false, ExcNotImplemented()); - return FiniteElementData::neither_element_dominates; + return FiniteElementDomination::neither_element_dominates; } -- 2.39.5