From: Peter Munch Date: Sat, 1 Aug 2020 12:45:54 +0000 (+0200) Subject: Use face_no in FE X-Git-Tag: v9.3.0-rc1~1191^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5080b6916e227801a5b0e83983e26141846d387a;p=dealii.git Use face_no in FE --- diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 9763abf10d..88b5787810 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -3129,9 +3129,12 @@ template inline std::pair FiniteElement::face_system_to_component_index( const unsigned int index, - const unsigned int) const + const unsigned int face_no) const { - AssertIndexRange(index, face_system_to_component_table[0].size()); + AssertIndexRange( + index, + face_system_to_component_table[this->n_unique_faces() == 1 ? 0 : face_no] + .size()); // in debug mode, check whether the // function is primitive, since @@ -3146,11 +3149,13 @@ FiniteElement::face_system_to_component_index( // // in 1d, the face index is equal // to the cell index - Assert(is_primitive(this->face_to_cell_index(index, 0)), + Assert(is_primitive(this->face_to_cell_index(index, face_no)), (typename FiniteElement::ExcShapeFunctionNotPrimitive( index))); - return face_system_to_component_table[0][index]; + return face_system_to_component_table[this->n_unique_faces() == 1 ? + 0 : + face_no][index]; } @@ -3170,10 +3175,14 @@ template inline std::pair, unsigned int> FiniteElement::face_system_to_base_index( const unsigned int index, - const unsigned int) const + const unsigned int face_no) const { - AssertIndexRange(index, face_system_to_base_table[0].size()); - return face_system_to_base_table[0][index]; + AssertIndexRange( + index, + face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no] + .size()); + return face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no] + [index]; } @@ -3301,7 +3310,7 @@ FiniteElement::get_associated_geometry_primitive( // are enumerated on the reference cell if (cell_dof_index < this->get_first_line_index()) return GeometryPrimitive::vertex; - else if (cell_dof_index < this->get_first_quad_index()) + else if (cell_dof_index < this->get_first_quad_index(0)) return GeometryPrimitive::line; else if (cell_dof_index < this->get_first_hex_index()) return GeometryPrimitive::quad; diff --git a/include/deal.II/fe/fe_poly_face.templates.h b/include/deal.II/fe/fe_poly_face.templates.h index fb94219cc1..06bc3ed006 100644 --- a/include/deal.II/fe/fe_poly_face.templates.h +++ b/include/deal.II/fe/fe_poly_face.templates.h @@ -141,16 +141,16 @@ FE_PolyFace::fill_fe_face_values( case 3: { // Fill data for quad shape functions - if (this->n_dofs_per_quad() != 0) + if (this->n_dofs_per_quad(face_no) != 0) { const unsigned int foffset = - this->get_first_quad_index() + - this->n_dofs_per_quad() * face_no; - for (unsigned int k = 0; k < this->n_dofs_per_quad(); ++k) + this->get_first_quad_index(face_no); + for (unsigned int k = 0; k < this->n_dofs_per_quad(face_no); + ++k) output_data.shape_values(foffset + k, i) = fe_data - .shape_values[k + this->get_first_face_quad_index()] - [i]; + .shape_values[k + this->get_first_face_quad_index( + face_no)][i]; } } DEAL_II_FALLTHROUGH; @@ -173,9 +173,10 @@ FE_PolyFace::fill_fe_face_values( line) * this->n_dofs_per_line() + k, - i) = fe_data.shape_values - [k + (line * this->n_dofs_per_line()) + - this->get_first_face_line_index()][i]; + i) = + fe_data.shape_values + [k + (line * this->n_dofs_per_line()) + + this->get_first_face_line_index(face_no)][i]; } } } diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 4d6b8e6899..96fb938532 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -965,7 +965,7 @@ public: */ virtual std::vector> hp_quad_dof_identities(const FiniteElement &fe_other, - const unsigned int = 0) const override; + const unsigned int face_no = 0) const override; /** * @copydoc FiniteElement::compare_for_domination() @@ -1162,7 +1162,8 @@ private: */ template std::vector> - hp_object_dof_identities(const FiniteElement &fe_other) const; + hp_object_dof_identities(const FiniteElement &fe_other, + const unsigned int face_no = 0) const; /** * Usually: Fields of cell-independent data. diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index 929d536ce4..a3c60daa3e 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -1246,7 +1246,8 @@ namespace FETools std::vector> & face_system_to_component_table, const FiniteElement &finite_element, - const bool do_tensor_product = true); + const bool do_tensor_product = true, + const unsigned int face_no = 0 /*TODO*/); } // namespace Compositing diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index 34c7d92789..cf8ac7d153 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -106,12 +106,17 @@ namespace FETools for (unsigned int i = 0; i < fes.size(); i++) if (multiplicities[i] > 0) { + // TODO: the implementation makes the assumption that all faces have + // the same number of dofs -> don't construct DPO but + // PrecomputedFiniteElementData + AssertDimension(fes[i]->n_unique_quads(), 1); + multiplied_dofs_per_vertex += fes[i]->n_dofs_per_vertex() * multiplicities[i]; multiplied_dofs_per_line += fes[i]->n_dofs_per_line() * multiplicities[i]; multiplied_dofs_per_quad += - fes[i]->n_dofs_per_quad() * multiplicities[i]; + fes[i]->n_dofs_per_quad(0) * multiplicities[i]; multiplied_dofs_per_hex += fes[i]->n_dofs_per_hex() * multiplicities[i]; @@ -287,12 +292,11 @@ namespace FETools for (unsigned int base = 0; base < fes.size(); ++base) for (unsigned int m = 0; m < multiplicities[base]; ++m) for (unsigned int local_index = 0; - local_index < fes[base]->n_dofs_per_quad(); + local_index < fes[base]->n_dofs_per_quad(quad_number); ++local_index, ++total_index) { const unsigned int index_in_base = - (fes[base]->n_dofs_per_quad() * quad_number + local_index + - fes[base]->get_first_quad_index()); + local_index + fes[base]->get_first_quad_index(quad_number); Assert(index_in_base < fes[base]->n_dofs_per_cell(), ExcInternalError()); @@ -509,12 +513,11 @@ namespace FETools comp_start += fes[base]->n_components() * do_tensor_product) for (unsigned int local_index = 0; - local_index < fes[base]->n_dofs_per_quad(); + local_index < fes[base]->n_dofs_per_quad(quad_number); ++local_index, ++total_index) { const unsigned int index_in_base = - (fes[base]->n_dofs_per_quad() * quad_number + local_index + - fes[base]->get_first_quad_index()); + local_index + fes[base]->get_first_quad_index(quad_number); Assert(comp_start + fes[base]->n_components() <= retval[total_index].size(), @@ -773,13 +776,13 @@ namespace FETools fe.base_element(base).n_components() * do_tensor_product) for (unsigned int local_index = 0; - local_index < fe.base_element(base).n_dofs_per_quad(); + local_index < + fe.base_element(base).n_dofs_per_quad(quad_number); ++local_index, ++total_index) { const unsigned int index_in_base = - (fe.base_element(base).n_dofs_per_quad() * quad_number + - local_index + - fe.base_element(base).get_first_quad_index()); + local_index + + fe.base_element(base).get_first_quad_index(quad_number); system_to_base_table[total_index] = std::make_pair(std::make_pair(base, m), index_in_base); @@ -855,7 +858,8 @@ namespace FETools std::vector> & face_system_to_component_table, const FiniteElement &fe, - const bool do_tensor_product) + const bool do_tensor_product, + const unsigned int face_no) { // Initialize index tables. do this in the same way as done for the cell // tables, except that we now loop over the objects of faces @@ -868,7 +872,8 @@ namespace FETools unsigned int total_index = 0; for (unsigned int vertex_number = 0; vertex_number < - ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0) + ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), + face_no) .n_vertices(); ++vertex_number) { @@ -906,12 +911,14 @@ namespace FETools { const unsigned int comp_in_base = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .first; const unsigned int comp = comp_start + comp_in_base; const unsigned int face_index_in_comp = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .second; face_system_to_component_table[total_index] = std::make_pair(comp, face_index_in_comp); @@ -925,7 +932,8 @@ namespace FETools // 2. Lines for (unsigned int line_number = 0; line_number < - ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0) + ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), + face_no) .n_lines(); ++line_number) { @@ -947,7 +955,7 @@ namespace FETools fe.base_element(base).get_first_line_index()); const unsigned int face_index_in_base = - (fe.base_element(base).get_first_face_line_index() + + (fe.base_element(base).get_first_face_line_index(face_no) + fe.base_element(base).n_dofs_per_line() * line_number + local_index); @@ -958,12 +966,14 @@ namespace FETools { const unsigned int comp_in_base = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .first; const unsigned int comp = comp_start + comp_in_base; const unsigned int face_index_in_comp = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .second; face_system_to_component_table[total_index] = std::make_pair(comp, face_index_in_comp); @@ -985,16 +995,16 @@ namespace FETools fe.base_element(base).n_components() * do_tensor_product) for (unsigned int local_index = 0; - local_index < fe.base_element(base).n_dofs_per_quad(); + local_index < fe.base_element(base).n_dofs_per_quad(face_no); ++local_index, ++total_index) { // do everything alike for this type of object const unsigned int index_in_base = (local_index + - fe.base_element(base).get_first_quad_index()); + fe.base_element(base).get_first_quad_index(face_no)); const unsigned int face_index_in_base = - (fe.base_element(base).get_first_face_quad_index() + + (fe.base_element(base).get_first_face_quad_index(face_no) + local_index); face_system_to_base_table[total_index] = @@ -1004,12 +1014,14 @@ namespace FETools { const unsigned int comp_in_base = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .first; const unsigned int comp = comp_start + comp_in_base; const unsigned int face_index_in_comp = fe.base_element(base) - .face_system_to_component_index(face_index_in_base) + .face_system_to_component_index(face_index_in_base, + face_no) .second; face_system_to_component_table[total_index] = std::make_pair(comp, face_index_in_comp); @@ -1019,7 +1031,7 @@ namespace FETools non_primitive_index; } } - Assert(total_index == fe.n_dofs_per_face(), ExcInternalError()); + Assert(total_index == fe.n_dofs_per_face(face_no), ExcInternalError()); Assert(total_index == face_system_to_component_table.size(), ExcInternalError()); Assert(total_index == face_system_to_base_table.size(), @@ -1934,11 +1946,13 @@ namespace FETools const unsigned int face_fine, const double threshold) { + const unsigned int face_no = face_coarse; + Assert(face_coarse == 0, ExcNotImplemented()); Assert(face_fine == 0, ExcNotImplemented()); const unsigned int nc = GeometryInfo::max_children_per_face; - const unsigned int n = fe.n_dofs_per_face(); + const unsigned int n = fe.n_dofs_per_face(face_no); const unsigned int nd = fe.n_components(); const unsigned int degree = fe.degree; @@ -1961,9 +1975,9 @@ namespace FETools { unsigned int face_dof = 0; for (unsigned int i = 0; - i < - ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0) - .n_vertices(); + i < ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), + face_no) + .n_vertices(); ++i) { const unsigned int offset_c = @@ -1980,10 +1994,9 @@ namespace FETools } } - for (unsigned int i = 1; - i <= - ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0) - .n_lines(); + for (unsigned int i = 1; i <= ReferenceCell::internal::Info::get_face( + fe.reference_cell_type(), face_no) + .n_lines(); ++i) { const unsigned int offset_c = @@ -2004,18 +2017,16 @@ namespace FETools if (dim == 3) { - const unsigned int offset_c = - fe.get_first_quad_index() + face_coarse * fe.n_dofs_per_quad(); - const unsigned int offset_f = - fe.get_first_quad_index() + face_fine * fe.n_dofs_per_quad(); - for (unsigned int j = 0; j < fe.n_dofs_per_quad(); ++j) + const unsigned int offset_c = fe.get_first_quad_index(face_coarse); + const unsigned int offset_f = fe.get_first_quad_index(face_fine); + for (unsigned int j = 0; j < fe.n_dofs_per_quad(face_no); ++j) { face_c_dofs[face_dof] = offset_c + j; face_f_dofs[face_dof] = offset_f + j; ++face_dof; } } - Assert(face_dof == fe.n_dofs_per_face(), ExcInternalError()); + Assert(face_dof == fe.n_dofs_per_face(face_no), ExcInternalError()); } // Set up meshes, one with a single diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h index 9f622655ad..b85649a2bf 100644 --- a/include/deal.II/fe/fe_tools_extrapolate.templates.h +++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h @@ -1412,9 +1412,8 @@ namespace FETools compute_all_non_local_data(dof2, u2_relevant); // exclude dofs on more refined ghosted cells - const FiniteElement &fe = dof2.get_fe(); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); - if (dofs_per_face > 0) + const FiniteElement &fe = dof2.get_fe(); + if (fe.max_dofs_per_face() > 0) { const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); std::vector indices(dofs_per_cell); @@ -1432,7 +1431,8 @@ namespace FETools const typename DoFHandler::cell_iterator neighbor = cell->neighbor(face); if (neighbor->level() != cell->level()) - for (unsigned int i = 0; i < dofs_per_face; ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); + ++i) { const types::global_dof_index index = indices[fe.face_to_cell_index(i, face)]; diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 53a850375a..abd2bd4655 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -107,20 +107,27 @@ FiniteElement::FiniteElement( for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) system_to_component_table[j] = std::pair(0, j); - face_system_to_component_table.resize(1); - face_system_to_component_table[0].resize(this->n_dofs_per_face()); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) - face_system_to_component_table[0][j] = - std::pair(0, j); + face_system_to_component_table.resize(this->n_unique_faces()); + for (unsigned int f = 0; f < this->n_unique_faces(); ++f) + { + face_system_to_component_table[f].resize(this->n_dofs_per_face(f)); + for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j) + face_system_to_component_table[f][j] = + std::pair(0, j); + } } for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j) system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j); - face_system_to_base_table.resize(1); - face_system_to_base_table[0].resize(this->n_dofs_per_face()); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) - face_system_to_base_table[0][j] = std::make_pair(std::make_pair(0U, 0U), j); + face_system_to_base_table.resize(this->n_unique_faces()); + for (unsigned int f = 0; f < this->n_unique_faces(); ++f) + { + face_system_to_base_table[f].resize(this->n_dofs_per_face(f)); + for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j) + face_system_to_base_table[f][j] = + std::make_pair(std::make_pair(0U, 0U), j); + } // Fill with default value; may be changed by constructor of derived class. base_to_block_indices.reinit(1, 1); @@ -141,17 +148,26 @@ FiniteElement::FiniteElement( FullMatrix()); } + if (dim == 3) { - adjust_quad_dof_index_for_face_orientation_table.resize(1); + adjust_quad_dof_index_for_face_orientation_table.resize( + this->n_unique_quads()); - adjust_quad_dof_index_for_face_orientation_table[0] = - Table<2, int>(this->n_dofs_per_quad(), 8); - adjust_quad_dof_index_for_face_orientation_table[0].fill(0); + for (unsigned int f = 0; f < this->n_unique_quads(); ++f) + { + adjust_quad_dof_index_for_face_orientation_table[f] = Table<2, int>( + this->n_dofs_per_quad(f), + ReferenceCell::internal::Info::get_cell(this->reference_cell_type()) + .face_reference_cell_type(f) == ReferenceCell::Type::Quad ? + 8 : + 6); + adjust_quad_dof_index_for_face_orientation_table[f].fill(0); + } } - unit_face_support_points.resize(1); - generalized_face_support_points.resize(1); + unit_face_support_points.resize(this->n_unique_faces()); + generalized_face_support_points.resize(this->n_unique_faces()); } @@ -558,7 +574,7 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, const auto &refence_cell = ReferenceCell::internal::Info::get_cell(this->reference_cell_type()); - AssertIndexRange(face_index, this->n_dofs_per_face()); + AssertIndexRange(face_index, this->n_dofs_per_face(face)); AssertIndexRange(face, refence_cell.n_faces()); // TODO: we could presumably solve the 3d case below using the @@ -576,7 +592,7 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, // other than standard orientation if ((face_orientation != true) || (face_flip != false) || (face_rotation != false)) - Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad() <= 1), + Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1), ExcMessage( "The function in this base class can not handle this case. " "Rather, the derived class you are using must provide " @@ -585,7 +601,7 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, // we need to distinguish between DoFs on vertices, lines and in 3d quads. // do so in a sequence of if-else statements - if (face_index < this->get_first_face_line_index()) + if (face_index < this->get_first_face_line_index(face)) // DoF is on a vertex { // get the number of the vertex on the face that corresponds to this DoF, @@ -604,12 +620,13 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, this->n_dofs_per_vertex() + dof_index_on_vertex); } - else if (face_index < this->get_first_face_quad_index()) + else if (face_index < this->get_first_face_quad_index(face)) // DoF is on a face { // do the same kind of translation as before. we need to only consider // DoFs on the lines, i.e., ignoring those on the vertices - const unsigned int index = face_index - this->get_first_face_line_index(); + const unsigned int index = + face_index - this->get_first_face_line_index(face); const unsigned int face_line = index / this->n_dofs_per_line(); const unsigned int dof_index_on_line = index % this->n_dofs_per_line(); @@ -629,10 +646,10 @@ FiniteElement::face_to_cell_index(const unsigned int face_index, Assert(dim >= 3, ExcInternalError()); // ignore vertex and line dofs - const unsigned int index = face_index - this->get_first_face_quad_index(); + const unsigned int index = + face_index - this->get_first_face_quad_index(face); - return (this->get_first_quad_index() + face * this->n_dofs_per_quad() + - index); + return (this->get_first_quad_index(face) + index); } } @@ -642,13 +659,11 @@ template unsigned int FiniteElement::adjust_quad_dof_index_for_face_orientation( const unsigned int index, - const unsigned int face_no, + const unsigned int face, const bool face_orientation, const bool face_flip, const bool face_rotation) const { - (void)face_no; - // general template for 1D and 2D: not // implemented. in fact, the function // shouldn't even be called unless we are @@ -666,12 +681,20 @@ FiniteElement::adjust_quad_dof_index_for_face_orientation( // in 3d), so we don't need the table, but // the function should also not have been // called - AssertIndexRange(index, this->n_dofs_per_quad()); - Assert(adjust_quad_dof_index_for_face_orientation_table[0].n_elements() == - 8 * this->n_dofs_per_quad(), + AssertIndexRange(index, this->n_dofs_per_quad(face)); + Assert(adjust_quad_dof_index_for_face_orientation_table + [this->n_unique_quads() == 1 ? 0 : face] + .n_elements() == + (ReferenceCell::internal::Info::get_cell(this->reference_cell_type()) + .face_reference_cell_type(face) == ReferenceCell::Type::Quad ? + 8 : + 6) * + this->n_dofs_per_quad(face), ExcInternalError()); - return index + adjust_quad_dof_index_for_face_orientation_table[0]( - index, 4 * face_orientation + 2 * face_flip + face_rotation); + return index + + adjust_quad_dof_index_for_face_orientation_table + [this->n_unique_quads() == 1 ? 0 : face]( + index, 4 * face_orientation + 2 * face_flip + face_rotation); } @@ -818,8 +841,14 @@ bool FiniteElement::constraints_are_implemented( const internal::SubfaceCase &subface_case) const { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + if (subface_case == internal::SubfaceCase::case_isotropic) - return (this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0); + return (this->n_dofs_per_face(face_no) == 0) || + (interface_constraints.m() != 0); else return false; } @@ -840,6 +869,12 @@ const FullMatrix & FiniteElement::constraints( const internal::SubfaceCase &subface_case) const { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + (void)face_no; + (void)subface_case; Assert(subface_case == internal::SubfaceCase::case_isotropic, ExcMessage("Constraints for this element are only implemented " @@ -847,7 +882,8 @@ FiniteElement::constraints( "(which is always the case in 2d, and in 3d requires " "that the neighboring cell of a coarse cell presents " "exactly four children on the common face).")); - Assert((this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0), + Assert((this->n_dofs_per_face(face_no) == 0) || + (interface_constraints.m() != 0), ExcMessage("The finite element for which you try to obtain " "hanging node constraints does not appear to " "implement them.")); @@ -866,17 +902,22 @@ template TableIndices<2> FiniteElement::interface_constraints_size() const { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + switch (dim) { case 1: return {0U, 0U}; case 2: return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(), - this->n_dofs_per_face()}; + this->n_dofs_per_face(face_no)}; case 3: return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() + - 4 * this->n_dofs_per_quad(), - this->n_dofs_per_face()}; + 4 * this->n_dofs_per_quad(face_no), + this->n_dofs_per_face(face_no)}; default: Assert(false, ExcNotImplemented()); } @@ -1069,38 +1110,45 @@ FiniteElement::unit_support_point(const unsigned int index) const template const std::vector> & FiniteElement::get_unit_face_support_points( - const unsigned int) const + const unsigned int face_no) const { // a finite element may define // support points, but only if // there are as many as there are // degrees of freedom on a face - Assert((unit_face_support_points[0].size() == 0) || - (unit_face_support_points[0].size() == this->n_dofs_per_face()), + Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() == 0) || + (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() == this->n_dofs_per_face(face_no)), ExcInternalError()); - return unit_face_support_points[0]; + return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]; } template bool -FiniteElement::has_face_support_points(const unsigned int) const +FiniteElement::has_face_support_points( + const unsigned int face_no) const { - return (unit_face_support_points[0].size() != 0); + return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() != 0); } template Point -FiniteElement::unit_face_support_point(const unsigned int index, - const unsigned int) const +FiniteElement::unit_face_support_point( + const unsigned int index, + const unsigned int face_no) const { - AssertIndexRange(index, this->n_dofs_per_face()); - Assert(unit_face_support_points[0].size() == this->n_dofs_per_face(), + AssertIndexRange(index, this->n_dofs_per_face(face_no)); + Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() == this->n_dofs_per_face(face_no), ExcFEHasNoSupportPoints()); - return unit_face_support_points[0][index]; + return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + [index]; } diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index c55552234b..0b11c719cf 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -88,17 +88,24 @@ FE_ABF::FE_ABF(const unsigned int deg) initialize_restriction(); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // TODO[TL]: for anisotropic refinement we will probably need a table of // submatrices with an array for each refine case std::vector> face_embeddings( 1 << (dim - 1), - FullMatrix(this->n_dofs_per_face(), this->n_dofs_per_face())); + FullMatrix(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no))); // TODO: Something goes wrong there. The error of the least squares fit // is to large ... // FETools::compute_face_embedding_matrices(*this, face_embeddings.data(), 0, // 0); - this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit((1 << (dim - 1)) * + this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (const auto &face_embedding : face_embeddings) for (unsigned int i = 0; i < face_embedding.m(); ++i) @@ -153,6 +160,11 @@ FE_ABF::initialize_support_points(const unsigned int deg) QGauss cell_quadrature(deg + 2); const unsigned int n_interior_points = cell_quadrature.size(); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + unsigned int n_face_points = (dim > 1) ? 1 : 0; // compute (deg+1)^(dim-1) for (unsigned int d = 1; d < dim; ++d) @@ -160,7 +172,7 @@ FE_ABF::initialize_support_points(const unsigned int deg) this->generalized_support_points.resize( GeometryInfo::faces_per_cell * n_face_points + n_interior_points); - this->generalized_face_support_points[0].resize(n_face_points); + this->generalized_face_support_points[face_no].resize(n_face_points); // These might be required when the faces contribution is computed @@ -194,7 +206,8 @@ FE_ABF::initialize_support_points(const unsigned int deg) for (unsigned int k = 0; k < n_face_points; ++k) { - this->generalized_face_support_points[0][k] = face_points.point(k); + this->generalized_face_support_points[face_no][k] = + face_points.point(k); // Compute its quadrature // contribution for each // moment. @@ -381,7 +394,8 @@ FE_ABF::initialize_restriction() for (unsigned int k = 0; k < n_face_points; ++k) for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell(); ++i_child) - for (unsigned int i_face = 0; i_face < this->n_dofs_per_face(); + for (unsigned int i_face = 0; + i_face < this->n_dofs_per_face(face); ++i_face) { // The quadrature @@ -389,13 +403,12 @@ FE_ABF::initialize_restriction() // subcell are NOT // transformed, so we // have to do it here. - this->restriction[iso][child](face * this->n_dofs_per_face() + - i_face, - i_child) += + this->restriction[iso][child]( + face * this->n_dofs_per_face(face) + i_face, i_child) += Utilities::fixed_power(.5) * q_sub.weight(k) * cached_values_face(i_child, k) * this->shape_value_component( - face * this->n_dofs_per_face() + i_face, + face * this->n_dofs_per_face(face) + i_face, q_sub.point(k), GeometryInfo::unit_normal_direction[face]); } @@ -419,9 +432,14 @@ FE_ABF::initialize_restriction() polynomials[dd] = new AnisotropicPolynomials(poly); } + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + QGauss q_cell(rt_order + 1); const unsigned int start_cell_dofs = - GeometryInfo::faces_per_cell * this->n_dofs_per_face(); + GeometryInfo::faces_per_cell * this->n_dofs_per_face(face_no); // Store shape values, since the // evaluation suffers if not @@ -561,14 +579,19 @@ FE_ABF::convert_generalized_support_point_values_to_dof_values( for (unsigned int k = 0; k < n_face_points; ++k) for (unsigned int i = 0; i < boundary_weights.size(1); ++i) { - nodal_values[i + face * this->n_dofs_per_face()] += + nodal_values[i + face * this->n_dofs_per_face(face)] += boundary_weights(k, i) * support_point_values[face * n_face_points + k][GeometryInfo< dim>::unit_normal_direction[face]]; } + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + const unsigned int start_cell_dofs = - GeometryInfo::faces_per_cell * this->n_dofs_per_face(); + GeometryInfo::faces_per_cell * this->n_dofs_per_face(face_no); const unsigned int start_cell_points = GeometryInfo::faces_per_cell * n_face_points; diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index dd1e4b469e..2efea8a775 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -80,12 +80,17 @@ FE_BDM::FE_BDM(const unsigned int deg) this->reinit_restriction_and_prolongation_matrices(true, true); FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.); + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face()); + face_embeddings[i].reinit(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.); - this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit((1 << (dim - 1)) * + this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (unsigned int d = 0; d < GeometryInfo::max_children_per_face; ++d) for (unsigned int i = 0; i < face_embeddings[d].m(); ++i) @@ -155,15 +160,15 @@ FE_BDM::convert_generalized_support_point_values_to_dof_values( // initialize_support_points() if (test_values_face.size() == 0) { - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i) nodal_values[dbase + i] = support_point_values[pbase + i] [GeometryInfo::unit_normal_direction[f]]; - pbase += this->n_dofs_per_face(); + pbase += this->n_dofs_per_face(f); } else { - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i) { double s = 0.; for (unsigned int k = 0; k < test_values_face.size(); ++k) @@ -175,11 +180,16 @@ FE_BDM::convert_generalized_support_point_values_to_dof_values( } pbase += test_values_face.size(); } - dbase += this->n_dofs_per_face(); + dbase += this->n_dofs_per_face(f); } + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + (void)face_no; + AssertDimension(dbase, - this->n_dofs_per_face() * GeometryInfo::faces_per_cell); + this->n_dofs_per_face(face_no) * + GeometryInfo::faces_per_cell); AssertDimension(pbase, this->generalized_support_points.size() - test_values_cell.size()); @@ -332,10 +342,15 @@ FE_BDM::initialize_support_points(const unsigned int deg) // considered later. In 2D, we can use point values. QGauss face_points(deg + 1); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // Copy the quadrature formula to the face points. - this->generalized_face_support_points[0].resize(face_points.size()); + this->generalized_face_support_points[face_no].resize(face_points.size()); for (unsigned int k = 0; k < face_points.size(); ++k) - this->generalized_face_support_points[0][k] = face_points.point(k); + this->generalized_face_support_points[face_no][k] = face_points.point(k); // In the interior, we only test with polynomials of degree up to // deg-2, thus we use deg points. Note that deg>=1 and the lowest @@ -356,12 +371,13 @@ FE_BDM::initialize_support_points(const unsigned int deg) k < face_points.size() * GeometryInfo::faces_per_cell; ++k) this->generalized_support_points[k] = faces.point( - k + QProjector::DataSetDescriptor::face(this->reference_cell_type(), - 0, - true, - false, - false, - this->n_dofs_per_face())); + k + + QProjector::DataSetDescriptor::face(this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face(face_no))); // Currently, for backward compatibility, we do not use moments, but // point values on faces in 2D. In 3D, this is impossible, since the diff --git a/source/fe/fe_bernstein.cc b/source/fe/fe_bernstein.cc index 8dee487809..86525d8659 100644 --- a/source/fe/fe_bernstein.cc +++ b/source/fe/fe_bernstein.cc @@ -112,11 +112,11 @@ FE_Bernstein::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // see if source is a Bernstein element if (const FE_Bernstein *source_fe = @@ -124,9 +124,9 @@ FE_Bernstein::get_subface_interpolation_matrix( { // have this test in here since a table of size 2x0 reports its size as // 0x0 - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); + this->n_dofs_per_face(face_no))); // Make sure that the element for which the DoFs should be constrained // is the one with the higher polynomial degree. Actually the procedure @@ -134,13 +134,13 @@ FE_Bernstein::get_subface_interpolation_matrix( // produced in that case might lead to problems in the hp procedures, // which use this method. Assert( - this->n_dofs_per_face() <= source_fe->n_dofs_per_face(), + this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); const Quadrature quad_face_support( FE_Q(QIterated<1>(QTrapez<1>(), source_fe->degree)) - .get_unit_face_support_points()); + .get_unit_face_support_points(face_no)); // Rule of thumb for FP accuracy, that can be expected for a given // polynomial degree. This value is used to cut off values close to @@ -163,10 +163,10 @@ FE_Bernstein::get_subface_interpolation_matrix( 0, subface); - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) { const Point &p = subface_quadrature.point(i); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = this->shape_value(this->face_to_cell_index(j, 0), p); @@ -185,11 +185,11 @@ FE_Bernstein::get_subface_interpolation_matrix( // make sure that the row sum of each of the matrices is 1 at this // point. this must be so since the shape functions sum up to 1 - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) { double sum = 0.; - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) sum += interpolation_matrix(j, i); Assert(std::fabs(sum - 1) < eps, ExcInternalError()); @@ -236,7 +236,7 @@ FE_Bernstein::hp_vertex_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc index 6f1cc3d2c5..f470b2e008 100644 --- a/source/fe/fe_enriched.cc +++ b/source/fe/fe_enriched.cc @@ -440,7 +440,6 @@ FE_Enriched::initialize( // If the system is not primitive, these have not been initialized by // FiniteElement this->system_to_component_table.resize(this->n_dofs_per_cell()); - this->face_system_to_component_table[0].resize(this->n_dofs_per_face()); FETools::Compositing::build_cell_tables(this->system_to_base_table, this->system_to_component_table, @@ -448,11 +447,21 @@ FE_Enriched::initialize( *this, false); - FETools::Compositing::build_face_tables( - this->face_system_to_base_table[0], - this->face_system_to_component_table[0], - *this, - false); + this->face_system_to_component_table.resize(this->n_unique_faces()); + + for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no) + { + this->face_system_to_component_table[0].resize( + this->n_dofs_per_face(face_no)); + + + FETools::Compositing::build_face_tables( + this->face_system_to_base_table[face_no], + this->face_system_to_component_table[face_no], + *this, + false, + face_no); + } } // restriction and prolongation matrices are built on demand @@ -468,8 +477,8 @@ FE_Enriched::initialize( // this FE sits on the boundary or not. Thus for moment just copy support // points from fe system: { - this->unit_support_points = fe_system->unit_support_points; - this->unit_face_support_points[0] = fe_system->unit_face_support_points[0]; + this->unit_support_points = fe_system->unit_support_points; + this->unit_face_support_points = fe_system->unit_face_support_points; } // take adjust_quad_dof_index_for_face_orientation_table from FESystem: diff --git a/source/fe/fe_face.cc b/source/fe/fe_face.cc index b9534a584e..d673916566 100644 --- a/source/fe/fe_face.cc +++ b/source/fe/fe_face.cc @@ -163,16 +163,16 @@ FE_FaceQ::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this function is similar to the respective method in FE_Q - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // see if source is a FaceQ element if (const FE_FaceQ *source_fe = @@ -184,13 +184,13 @@ FE_FaceQ::get_subface_interpolation_matrix( // produced in that case might lead to problems in the hp procedures, // which use this method. Assert( - this->n_dofs_per_face() <= source_fe->n_dofs_per_face(), + this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); // generate a quadrature with the unit face support points. const Quadrature face_quadrature( - source_fe->get_unit_face_support_points()); + source_fe->get_unit_face_support_points(face_no)); // Rule of thumb for FP accuracy, that can be expected for a given // polynomial degree. This value is used to cut off values close to @@ -199,7 +199,7 @@ FE_FaceQ::get_subface_interpolation_matrix( // compute the interpolation matrix by simply taking the value at the // support points. - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) { const Point p = subface == numbers::invalid_unsigned_int ? @@ -207,7 +207,7 @@ FE_FaceQ::get_subface_interpolation_matrix( GeometryInfo::child_to_cell_coordinates( face_quadrature.point(i), subface); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = this->poly_space.compute_value(j, p); @@ -225,11 +225,11 @@ FE_FaceQ::get_subface_interpolation_matrix( // make sure that the row sum of each of the matrices is 1 at this // point. this must be so since the shape functions sum up to 1 - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) { double sum = 0.; - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) sum += interpolation_matrix(j, i); Assert(std::fabs(sum - 1) < eps, ExcInternalError()); @@ -254,7 +254,7 @@ FE_FaceQ::has_support_on_face( const unsigned int shape_index, const unsigned int face_index) const { - return (face_index == (shape_index / this->n_dofs_per_face())); + return (face_index == (shape_index / this->n_dofs_per_face(face_index))); } @@ -340,7 +340,8 @@ FE_FaceQ::hp_line_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && + fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of @@ -418,7 +419,8 @@ FE_FaceQ::hp_quad_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && + fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of @@ -574,15 +576,17 @@ FE_FaceQ<1, spacedim>::get_subface_interpolation_matrix( const FiniteElement<1, spacedim> &x_source_fe, const unsigned int /*subface*/, FullMatrix &interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { (void)x_source_fe; - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + (void)face_no; + + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); interpolation_matrix(0, 0) = 1.; } @@ -798,7 +802,7 @@ FE_FaceP::has_support_on_face( const unsigned int shape_index, const unsigned int face_index) const { - return (face_index == (shape_index / this->n_dofs_per_face())); + return (face_index == (shape_index / this->n_dofs_per_face(face_index))); } @@ -888,16 +892,16 @@ FE_FaceP::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this function is similar to the respective method in FE_Q - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // see if source is a FaceP element if (const FE_FaceP *source_fe = @@ -909,7 +913,7 @@ FE_FaceP::get_subface_interpolation_matrix( // produced in that case might lead to problems in the hp procedures, // which use this method. Assert( - this->n_dofs_per_face() <= source_fe->n_dofs_per_face(), + this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); @@ -924,7 +928,7 @@ FE_FaceP::get_subface_interpolation_matrix( const double eps = 2e-13 * (this->degree + 1) * (dim - 1); FullMatrix mass(face_quadrature.size(), - source_fe->n_dofs_per_face()); + source_fe->n_dofs_per_face(face_no)); for (unsigned int k = 0; k < face_quadrature.size(); ++k) { @@ -934,18 +938,18 @@ FE_FaceP::get_subface_interpolation_matrix( GeometryInfo::child_to_cell_coordinates( face_quadrature.point(k), subface); - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) mass(k, j) = source_fe->poly_space.compute_value(j, p); } Householder H(mass); Vector v_in(face_quadrature.size()); - Vector v_out(source_fe->n_dofs_per_face()); + Vector v_out(source_fe->n_dofs_per_face(face_no)); // compute the interpolation matrix by evaluating on the fine side and // then solving the least squares problem - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) { for (unsigned int k = 0; k < face_quadrature.size(); ++k) { @@ -960,7 +964,7 @@ FE_FaceP::get_subface_interpolation_matrix( (void)result; Assert(result < 1e-12, FETools::ExcLeastSquaresError(result)); - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) { double matrix_entry = v_out(j); @@ -995,7 +999,7 @@ FE_FaceP::get_constant_modes() const { Table<2, bool> constant_modes(1, this->n_dofs_per_cell()); for (unsigned int face : GeometryInfo::face_indices()) - constant_modes(0, face * this->n_dofs_per_face()) = true; + constant_modes(0, face * this->n_dofs_per_face(face)) = true; return std::pair, std::vector>( constant_modes, std::vector(1, 0)); } diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 6f46551ccb..caa4f1f39d 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -108,8 +108,14 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) #endif FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face()); + face_embeddings[i].reinit(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); FETools::compute_face_embedding_matrices( *this, @@ -128,14 +134,15 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) case 2: { - this->interface_constraints.reinit(2 * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face; ++i) - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) - this->interface_constraints(i * this->n_dofs_per_face() + j, + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) + this->interface_constraints(i * this->n_dofs_per_face(face_no) + + j, k) = face_embeddings[i](j, k); break; @@ -143,16 +150,16 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) case 3: { - this->interface_constraints.reinit(4 * (this->n_dofs_per_face() - - this->degree), - this->n_dofs_per_face()); + this->interface_constraints.reinit( + 4 * (this->n_dofs_per_face(face_no) - this->degree), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = this->degree; j < 2 * this->degree; ++j, ++target_row) - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) this->interface_constraints(target_row, k) = face_embeddings[2 * i](j, k); @@ -160,7 +167,7 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) for (unsigned int j = 3 * this->degree; j < GeometryInfo<3>::lines_per_face * this->degree; ++j, ++target_row) - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) this->interface_constraints(target_row, k) = face_embeddings[i](j, k); @@ -169,7 +176,8 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) for (unsigned int k = i * this->degree; k < (i + 1) * this->degree; ++k, ++target_row) - for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l) + for (unsigned int l = 0; l < this->n_dofs_per_face(face_no); + ++l) this->interface_constraints(target_row, l) = face_embeddings[i + 2 * j](k, l); @@ -178,7 +186,8 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) for (unsigned int k = (i + 2) * this->degree; k < (i + 3) * this->degree; ++k, ++target_row) - for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l) + for (unsigned int l = 0; l < this->n_dofs_per_face(face_no); + ++l) this->interface_constraints(target_row, l) = face_embeddings[2 * i + j](k, l); @@ -186,9 +195,9 @@ FE_Nedelec::FE_Nedelec(const unsigned int order) ++i) for (unsigned int j = GeometryInfo<3>::lines_per_face * this->degree; - j < this->n_dofs_per_face(); + j < this->n_dofs_per_face(face_no); ++j, ++target_row) - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) this->interface_constraints(target_row, k) = face_embeddings[i](j, k); @@ -254,6 +263,11 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) { const int dim = 2; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // Create polynomial basis. const std::vector> &lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis(order + 1); @@ -273,11 +287,11 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) QProjector::project_to_all_faces(this->reference_cell_type(), reference_edge_quadrature); - this->generalized_face_support_points[0].resize(n_edge_points); + this->generalized_face_support_points[face_no].resize(n_edge_points); // Create face support points. for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point) - this->generalized_face_support_points[0][q_point] = + this->generalized_face_support_points[face_no][q_point] = reference_edge_quadrature.point(q_point); if (order > 0) @@ -310,7 +324,7 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order) boundary_weights(q_point, i) = reference_edge_quadrature.weight(q_point) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[0][q_point](0)); + this->generalized_face_support_points[face_no][q_point](0)); } for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) @@ -347,6 +361,11 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) { const int dim = 3; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // Create polynomial basis. const std::vector> &lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis(order + 1); @@ -380,8 +399,8 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) boundary_weights.reinit(n_edge_points + n_face_points, 2 * (order + 1) * order); - this->generalized_face_support_points[0].resize(4 * n_edge_points + - n_face_points); + this->generalized_face_support_points[face_no].resize(4 * n_edge_points + + n_face_points); this->generalized_support_points.resize(n_boundary_points + n_interior_points); @@ -391,8 +410,9 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; ++line) - this->generalized_face_support_points[0][line * n_edge_points + - q_point] = + this + ->generalized_face_support_points[face_no][line * n_edge_points + + q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( ReferenceCell::get_hypercube(dim - 1), @@ -421,14 +441,14 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) boundary_weights(q_point, i) = reference_edge_quadrature.weight(q_point) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[0][q_point](1)); + this->generalized_face_support_points[face_no][q_point](1)); } // Create support points on faces. for (unsigned int q_point = 0; q_point < n_face_points; ++q_point) { - this - ->generalized_face_support_points[0][q_point + 4 * n_edge_points] = + this->generalized_face_support_points[face_no] + [q_point + 4 * n_edge_points] = reference_face_quadrature.point(q_point); for (unsigned int i = 0; i <= order; ++i) @@ -437,24 +457,20 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) boundary_weights(q_point + n_edge_points, 2 * (i * order + j)) = reference_face_quadrature.weight(q_point) * lobatto_polynomials_grad[i].value( - this->generalized_face_support_points[0][q_point + - 4 * n_edge_points]( - 0)) * + this->generalized_face_support_points + [face_no][q_point + 4 * n_edge_points](0)) * lobatto_polynomials[j + 2].value( - this->generalized_face_support_points[0][q_point + - 4 * n_edge_points]( - 1)); + this->generalized_face_support_points + [face_no][q_point + 4 * n_edge_points](1)); boundary_weights(q_point + n_edge_points, 2 * (i * order + j) + 1) = reference_face_quadrature.weight(q_point) * lobatto_polynomials_grad[i].value( - this->generalized_face_support_points[0][q_point + - 4 * n_edge_points]( - 1)) * + this->generalized_face_support_points + [face_no][q_point + 4 * n_edge_points](1)) * lobatto_polynomials[j + 2].value( - this->generalized_face_support_points[0][q_point + - 4 * n_edge_points]( - 0)); + this->generalized_face_support_points + [face_no][q_point + 4 * n_edge_points](0)); } } @@ -486,7 +502,7 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) else { - this->generalized_face_support_points[0].resize(4 * n_edge_points); + this->generalized_face_support_points[face_no].resize(4 * n_edge_points); this->generalized_support_points.resize( GeometryInfo::lines_per_cell * n_edge_points); @@ -495,8 +511,9 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; ++line) - this->generalized_face_support_points[0][line * n_edge_points + - q_point] = + this + ->generalized_face_support_points[face_no][line * n_edge_points + + q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( ReferenceCell::get_hypercube(dim - 1), @@ -2453,20 +2470,21 @@ void FE_Nedelec::get_face_interpolation_matrix( const FiniteElement &source, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { + (void)face_no; // this is only implemented, if the // source FE is also a // Nedelec element AssertThrow((source.get_name().find("FE_Nedelec<") == 0) || (dynamic_cast *>(&source) != nullptr), (typename FiniteElement::ExcInterpolationNotImplemented())); - Assert(interpolation_matrix.m() == source.n_dofs_per_face(), + Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - source.n_dofs_per_face())); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + source.n_dofs_per_face(face_no))); + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); + this->n_dofs_per_face(face_no))); // ok, source is a Nedelec element, so // we will be able to do the work @@ -2484,7 +2502,7 @@ FE_Nedelec::get_face_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); interpolation_matrix = 0; @@ -2560,7 +2578,7 @@ FE_Nedelec::get_subface_interpolation_matrix( const FiniteElement &source, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is only implemented, if the // source FE is also a @@ -2568,12 +2586,12 @@ FE_Nedelec::get_subface_interpolation_matrix( AssertThrow((source.get_name().find("FE_Nedelec<") == 0) || (dynamic_cast *>(&source) != nullptr), typename FiniteElement::ExcInterpolationNotImplemented()); - Assert(interpolation_matrix.m() == source.n_dofs_per_face(), + Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - source.n_dofs_per_face())); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + source.n_dofs_per_face(face_no))); + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); + this->n_dofs_per_face(face_no))); // ok, source is a Nedelec element, so // we will be able to do the work @@ -2591,7 +2609,7 @@ FE_Nedelec::get_subface_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); interpolation_matrix = 0.0; // Perform projection-based interpolation @@ -2605,7 +2623,8 @@ FE_Nedelec::get_subface_interpolation_matrix( { case 2: { - for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof) + for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no); + ++dof) for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { @@ -2653,7 +2672,8 @@ FE_Nedelec::get_subface_interpolation_matrix( Vector solution(source_fe.degree - 1); Vector system_rhs(source_fe.degree - 1); - for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof) + for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no); + ++dof) { system_rhs = 0.0; @@ -2700,7 +2720,8 @@ FE_Nedelec::get_subface_interpolation_matrix( {0.0, 1.0}, {1.0, 1.0}}; - for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof) + for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no); + ++dof) for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { @@ -2769,7 +2790,8 @@ FE_Nedelec::get_subface_interpolation_matrix( GeometryInfo::lines_per_face); Vector tmp(GeometryInfo::lines_per_face); - for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof) + for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no); + ++dof) { system_rhs = 0.0; @@ -2892,7 +2914,8 @@ FE_Nedelec::get_subface_interpolation_matrix( system_rhs.reinit(system_matrix_inv.m(), 2); tmp.reinit(2); - for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof) + for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no); + ++dof) { system_rhs = 0.0; @@ -3112,6 +3135,11 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( const std::vector> &support_point_values, std::vector & nodal_values) const { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + const unsigned int deg = this->degree - 1; Assert(support_point_values.size() == this->generalized_support_points.size(), ExcDimensionMismatch(support_point_values.size(), @@ -3186,7 +3214,8 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( system_matrix(i, j) += boundary_weights(q_point, j) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[0][q_point](0)); + this->generalized_face_support_points[face_no][q_point]( + 0)); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3465,7 +3494,8 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( system_matrix(i, j) += boundary_weights(q_point, j) * lobatto_polynomials_grad[i + 1].value( - this->generalized_face_support_points[0][q_point](1)); + this->generalized_face_support_points[face_no][q_point]( + 1)); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3544,10 +3574,10 @@ FE_Nedelec::convert_generalized_support_point_values_to_dof_values( 2 * (k * (this->degree - 1) + l)) * legendre_polynomials[i].value( this->generalized_face_support_points - [0][q_point + 4 * n_edge_points](0)) * + [face_no][q_point + 4 * n_edge_points](0)) * lobatto_polynomials[j + 2].value( this->generalized_face_support_points - [0][q_point + 4 * n_edge_points](1)); + [face_no][q_point + 4 * n_edge_points](1)); system_matrix_inv.reinit(system_matrix.m(), system_matrix.m()); system_matrix_inv.invert(system_matrix); diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc index 7a9d9bb646..f2a6e5d049 100644 --- a/source/fe/fe_nedelec_sz.cc +++ b/source/fe/fe_nedelec_sz.cc @@ -147,7 +147,9 @@ FE_NedelecSZ::get_data( const unsigned int faces_per_cell = GeometryInfo::faces_per_cell; const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell; - const unsigned int n_face_dofs = this->n_dofs_per_quad() * faces_per_cell; + + // we assume that all quads have the same numer of dofs + const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell; const UpdateFlags flags(data.update_each); const unsigned int n_q_points = quadrature.size(); @@ -1695,7 +1697,8 @@ FE_NedelecSZ::fill_face_values( // Loop through quad points: for (unsigned int m = 0; m < faces_per_cell; ++m) { - const unsigned int shift_m(m * this->n_dofs_per_quad()); + // we assume that all quads have the same numer of dofs + const unsigned int shift_m(m * this->n_dofs_per_quad(0)); // Calculate the offsets for each face-based shape function: // // Type-1 (gradients) diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index d6f1794bf3..6df444121b 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -83,11 +83,11 @@ namespace internal const unsigned int nn = cell->neighbor_face_no(f); if (nn < GeometryInfo::faces_per_cell / 2) - for (unsigned int j = 0; j < fe.n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j) { const unsigned int cell_j = fe.face_to_cell_index(j, f); - Assert(f * fe.n_dofs_per_face() + j < face_sign.size(), + Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(), ExcInternalError()); Assert(mapping_kind.size() == 1 || cell_j < mapping_kind.size(), @@ -98,7 +98,7 @@ namespace internal if ((mapping_kind.size() > 1 ? mapping_kind[cell_j] : mapping_kind[0]) == mapping_raviart_thomas) - face_sign[f * fe.n_dofs_per_face() + j] = -1.0; + face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0; } } } diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 38c301c2ba..3d6cf65dd4 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -624,11 +624,11 @@ FE_Q_Base::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // see if source is a Q element if (const FE_Q_Base *source_fe = @@ -637,9 +637,9 @@ FE_Q_Base::get_subface_interpolation_matrix( { // have this test in here since a table of size 2x0 reports its size as // 0x0 - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); + this->n_dofs_per_face(face_no))); // Make sure that the element for which the DoFs should be constrained // is the one with the higher polynomial degree. Actually the procedure @@ -647,13 +647,13 @@ FE_Q_Base::get_subface_interpolation_matrix( // produced in that case might lead to problems in the hp procedures, // which use this method. Assert( - this->n_dofs_per_face() <= source_fe->n_dofs_per_face(), + this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); // generate a point on this cell and evaluate the shape functions there const Quadrature quad_face_support( - source_fe->get_unit_face_support_points()); + source_fe->get_unit_face_support_points(face_no)); // Rule of thumb for FP accuracy, that can be expected for a given // polynomial degree. This value is used to cut off values close to @@ -674,11 +674,11 @@ FE_Q_Base::get_subface_interpolation_matrix( quad_face_support, 0, subface); - for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i) { const Point &p = subface_quadrature.point(i); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = this->shape_value(this->face_to_cell_index(j, 0), p); @@ -697,11 +697,11 @@ FE_Q_Base::get_subface_interpolation_matrix( // make sure that the row sum of each of the matrices is 1 at this // point. this must be so since the shape functions sum up to 1 - for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j) { double sum = 0.; - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) sum += interpolation_matrix(j, i); Assert(std::fabs(sum - 1) < eps, ExcInternalError()); @@ -750,7 +750,7 @@ FE_Q_Base::hp_vertex_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of @@ -815,7 +815,7 @@ FE_Q_Base::hp_line_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of @@ -886,7 +886,7 @@ FE_Q_Base::hp_quad_dof_identities( // equivalencies to be recorded return std::vector>(); } - else if (fe_other.n_dofs_per_face() == 0) + else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0) { // if the other element has no elements on faces at all, // then it would be impossible to enforce any kind of @@ -946,7 +946,12 @@ FE_Q_Base::initialize_unit_face_support_points( if (dim == 1) return; - this->unit_face_support_points[0].resize( + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + + this->unit_face_support_points[face_no].resize( Utilities::fixed_power(q_degree + 1)); // find renumbering of faces and assign from values of quadrature @@ -962,9 +967,9 @@ FE_Q_Base::initialize_unit_face_support_points( // The only thing we have to do is reorder the points from tensor // product order to the order in which we enumerate DoFs on cells - this->unit_face_support_points[0].resize(support_quadrature.size()); + this->unit_face_support_points[face_no].resize(support_quadrature.size()); for (unsigned int k = 0; k < support_quadrature.size(); ++k) - this->unit_face_support_points[0][face_index_map[k]] = + this->unit_face_support_points[face_no][face_index_map[k]] = support_quadrature.point(k); } @@ -979,12 +984,17 @@ FE_Q_Base:: if (dim < 3) return; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + Assert(this->adjust_quad_dof_index_for_face_orientation_table[0] - .n_elements() == 8 * this->n_dofs_per_quad(), + .n_elements() == 8 * this->n_dofs_per_quad(face_no), ExcInternalError()); const unsigned int n = q_degree - 1; - Assert(n * n == this->n_dofs_per_quad(), ExcInternalError()); + Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError()); // the dofs on a face are connected to a n x n matrix. for example, for // degree==4 we have the following dofs on a quad @@ -1004,34 +1014,42 @@ FE_Q_Base:: // rotated and mirrored numbers. - for (unsigned int local = 0; local < this->n_dofs_per_quad(); ++local) + for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local) // face support points are in lexicographic ordering with x running // fastest. invert that (y running fastest) { unsigned int i = local % n, j = local / n; // face_orientation=false, face_flip=false, face_rotation=false - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 0) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 0) = j + i * n - local; // face_orientation=false, face_flip=false, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 1) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 1) = i + (n - 1 - j) * n - local; // face_orientation=false, face_flip=true, face_rotation=false - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 2) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 2) = (n - 1 - j) + (n - 1 - i) * n - local; // face_orientation=false, face_flip=true, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 3) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 3) = (n - 1 - i) + j * n - local; // face_orientation=true, face_flip=false, face_rotation=false - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 4) = 0; + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 4) = 0; // face_orientation=true, face_flip=false, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 5) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 5) = j + (n - 1 - i) * n - local; // face_orientation=true, face_flip=true, face_rotation=false - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 6) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 6) = (n - 1 - i) + (n - 1 - j) * n - local; // face_orientation=true, face_flip=true, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table[0](local, 7) = + this->adjust_quad_dof_index_for_face_orientation_table[face_no](local, + 7) = (n - 1 - j) + i * n - local; } @@ -1052,7 +1070,7 @@ FE_Q_Base::face_to_cell_index( const bool face_flip, const bool face_rotation) const { - AssertIndexRange(face_index, this->n_dofs_per_face()); + AssertIndexRange(face_index, this->n_dofs_per_face(face)); AssertIndexRange(face, GeometryInfo::faces_per_cell); // TODO: we could presumably solve the 3d case below using the @@ -1066,7 +1084,7 @@ FE_Q_Base::face_to_cell_index( // we need to distinguish between DoFs on vertices, lines and in 3d quads. // do so in a sequence of if-else statements - if (face_index < this->get_first_face_line_index()) + if (face_index < this->get_first_face_line_index(face)) // DoF is on a vertex { // get the number of the vertex on the face that corresponds to this DoF, @@ -1082,12 +1100,13 @@ FE_Q_Base::face_to_cell_index( this->n_dofs_per_vertex() + dof_index_on_vertex); } - else if (face_index < this->get_first_face_quad_index()) + else if (face_index < this->get_first_face_quad_index(face)) // DoF is on a face { // do the same kind of translation as before. we need to only consider // DoFs on the lines, i.e., ignoring those on the vertices - const unsigned int index = face_index - this->get_first_face_line_index(); + const unsigned int index = + face_index - this->get_first_face_line_index(face); const unsigned int face_line = index / this->n_dofs_per_line(); const unsigned int dof_index_on_line = index % this->n_dofs_per_line(); @@ -1142,17 +1161,17 @@ FE_Q_Base::face_to_cell_index( Assert(dim >= 3, ExcInternalError()); // ignore vertex and line dofs - const unsigned int index = face_index - this->get_first_face_quad_index(); + const unsigned int index = + face_index - this->get_first_face_quad_index(face); // the same is true here as above for the 3d case -- someone will // just have to draw a bunch of pictures. in the meantime, // we can implement the Q2 case in which it is simple - Assert((this->n_dofs_per_quad() <= 1) || + Assert((this->n_dofs_per_quad(face) <= 1) || ((face_orientation == true) && (face_flip == false) && (face_rotation == false)), ExcNotImplemented()); - return (this->get_first_quad_index() + face * this->n_dofs_per_quad() + - index); + return (this->get_first_quad_index(face) + index); } } @@ -1544,7 +1563,8 @@ FE_Q_Base::has_support_on_face( // first, special-case interior shape functions, since they have no support // no-where on the boundary - if (((dim == 2) && (shape_index >= this->get_first_quad_index())) || + if (((dim == 2) && + (shape_index >= this->get_first_quad_index(0 /*first quad*/))) || ((dim == 3) && (shape_index >= this->get_first_hex_index()))) return false; @@ -1565,7 +1585,7 @@ FE_Q_Base::has_support_on_face( return false; } - else if (shape_index < this->get_first_quad_index()) + else if (shape_index < this->get_first_quad_index(0 /*first quad*/)) // ok, dof is on a line { const unsigned int line_index = @@ -1596,7 +1616,8 @@ FE_Q_Base::has_support_on_face( // dof is on a quad { const unsigned int quad_index = - (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad(); + (shape_index - this->get_first_quad_index(0)) / + this->n_dofs_per_quad(face_index); // this won't work Assert(static_cast(quad_index) < static_cast(GeometryInfo::quads_per_cell), ExcInternalError()); diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 3af55b0c28..670ec26d5e 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -307,8 +307,13 @@ FE_Q_Hierarchical::hp_quad_dof_identities( // one is an FE_Nothing. if (dynamic_cast *>(&fe_other) != nullptr) { - const unsigned int this_dpq = this->n_dofs_per_quad(); - const unsigned int other_dpq = fe_other.n_dofs_per_quad(); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + + const unsigned int this_dpq = this->n_dofs_per_quad(face_no); + const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no); // we deal with hierarchical 1d polynomials where dofs are enumerated // increasingly. Thus we return a vector of pairs for the first N-1, where @@ -927,7 +932,7 @@ void FE_Q_Hierarchical::get_face_interpolation_matrix( const FiniteElement &x_source_fe, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is only implemented, if the // source FE is also a @@ -938,12 +943,12 @@ FE_Q_Hierarchical::get_face_interpolation_matrix( nullptr), (typename FiniteElement::ExcInterpolationNotImplemented())); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // ok, source is a Q_Hierarchical element, so // we will be able to do the work @@ -962,7 +967,7 @@ FE_Q_Hierarchical::get_face_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); interpolation_matrix = 0; @@ -976,7 +981,7 @@ FE_Q_Hierarchical::get_face_interpolation_matrix( // element, which corresponds to 1 on diagonal of the matrix. // DoFs which correspond to higher polynomials // are zeroed (zero rows in the matrix). - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) interpolation_matrix(i, i) = 1; break; @@ -1016,7 +1021,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is only implemented, if the // source FE is also a @@ -1027,12 +1032,12 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( nullptr), (typename FiniteElement::ExcInterpolationNotImplemented())); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // ok, source is a Q_Hierarchical element, so // we will be able to do the work @@ -1050,7 +1055,7 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), (typename FiniteElement::ExcInterpolationNotImplemented())); switch (dim) @@ -1065,7 +1070,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(1, 0) = 0.5; interpolation_matrix(1, 1) = 0.5; - for (unsigned int dof = 2; dof < this->n_dofs_per_face();) + for (unsigned int dof = 2; + dof < this->n_dofs_per_face(face_no);) { interpolation_matrix(1, dof) = -1.0; dof = dof + 2; @@ -1073,14 +1079,16 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( int factorial_i = 1; - for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 2; i < this->n_dofs_per_face(face_no); + ++i) { interpolation_matrix(i, i) = std::pow(0.5, i); factorial_i *= i; int factorial_j = factorial_i; int factorial_ij = 1; - for (unsigned int j = i + 1; j < this->n_dofs_per_face(); + for (unsigned int j = i + 1; + j < this->n_dofs_per_face(face_no); ++j) { factorial_ij *= j - i; @@ -1106,7 +1114,8 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( interpolation_matrix(0, 0) = 0.5; interpolation_matrix(0, 1) = 0.5; - for (unsigned int dof = 2; dof < this->n_dofs_per_face();) + for (unsigned int dof = 2; + dof < this->n_dofs_per_face(face_no);) { interpolation_matrix(0, dof) = -1.0; dof = dof + 2; @@ -1116,14 +1125,16 @@ FE_Q_Hierarchical::get_subface_interpolation_matrix( int factorial_i = 1; - for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 2; i < this->n_dofs_per_face(face_no); + ++i) { interpolation_matrix(i, i) = std::pow(0.5, i); factorial_i *= i; int factorial_j = factorial_i; int factorial_ij = 1; - for (unsigned int j = i + 1; j < this->n_dofs_per_face(); + for (unsigned int j = i + 1; + j < this->n_dofs_per_face(face_no); ++j) { factorial_ij *= j - i; @@ -1903,12 +1914,17 @@ FE_Q_Hierarchical::initialize_generalized_face_support_points() { const unsigned int codim = dim - 1; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // number of points: (degree+1)^codim unsigned int n = this->degree + 1; for (unsigned int i = 1; i < codim; ++i) n *= this->degree + 1; - this->generalized_face_support_points[0].resize(n); + this->generalized_face_support_points[face_no].resize(n); Point p; @@ -1941,7 +1957,8 @@ FE_Q_Hierarchical::initialize_generalized_face_support_points() else p(2) = .5; } - this->generalized_face_support_points[0][face_renumber[k++]] = p; + this->generalized_face_support_points[face_no][face_renumber[k++]] = + p; } } @@ -2029,7 +2046,7 @@ FE_Q_Hierarchical::hierarchic_to_fe_q_hierarchical_numbering( for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i) h2l[next_index++] = n + 2 + i; // inside quad - Assert(fe.n_dofs_per_quad() == + Assert(fe.n_dofs_per_quad(0 /*only one quad in 2D*/) == fe.n_dofs_per_line() * fe.n_dofs_per_line(), ExcInternalError()); for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i) @@ -2086,8 +2103,14 @@ FE_Q_Hierarchical::hierarchic_to_fe_q_hierarchical_numbering( for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i) h2l[next_index++] = (2 + i) * n2 + n + 1; + // TODO: the implementation makes the assumption that all faces have + // the same number of dofs + AssertDimension(fe.n_unique_faces(), 1); + const unsigned int face_no = 0; + (void)face_no; + // inside quads - Assert(fe.n_dofs_per_quad() == + Assert(fe.n_dofs_per_quad(face_no) == fe.n_dofs_per_line() * fe.n_dofs_per_line(), ExcInternalError()); // left face @@ -2117,7 +2140,7 @@ FE_Q_Hierarchical::hierarchic_to_fe_q_hierarchical_numbering( // inside hex Assert(fe.n_dofs_per_hex() == - fe.n_dofs_per_quad() * fe.n_dofs_per_line(), + fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(), ExcInternalError()); for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i) for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j) @@ -2192,7 +2215,8 @@ FE_Q_Hierarchical::has_support_on_face(const unsigned int shape_index, // shape functions, since they // have no support no-where on // the boundary - if (((dim == 2) && (shape_index >= this->get_first_quad_index())) || + if (((dim == 2) && (shape_index >= + this->get_first_quad_index(0 /*only one quad in 2D*/))) || ((dim == 3) && (shape_index >= this->get_first_hex_index()))) return false; @@ -2214,7 +2238,7 @@ FE_Q_Hierarchical::has_support_on_face(const unsigned int shape_index, return true; return false; } - else if (shape_index < this->get_first_quad_index()) + else if (shape_index < this->get_first_quad_index(0)) // ok, dof is on a line { const unsigned int line_index = @@ -2231,7 +2255,8 @@ FE_Q_Hierarchical::has_support_on_face(const unsigned int shape_index, // dof is on a quad { const unsigned int quad_index = - (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad(); + (shape_index - this->get_first_quad_index(0 /*first quad*/)) / + this->n_dofs_per_quad(face_index); Assert(static_cast(quad_index) < static_cast(GeometryInfo::quads_per_cell), ExcInternalError()); diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 5b4334ae6c..78d3271a9a 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -87,17 +87,24 @@ FE_RaviartThomas::FE_RaviartThomas(const unsigned int deg) FETools::compute_embedding_matrices(*this, this->prolongation); initialize_restriction(); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // TODO[TL]: for anisotropic refinement we will probably need a table of // submatrices with an array for each refine case FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face()); + face_embeddings[i].reinit(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0); - this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit((1 << (dim - 1)) * + this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (unsigned int d = 0; d < GeometryInfo::max_children_per_face; ++d) for (unsigned int i = 0; i < face_embeddings[d].m(); ++i) @@ -152,6 +159,11 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) QGauss cell_quadrature(deg + 1); const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0; + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + unsigned int n_face_points = (dim > 1) ? 1 : 0; // compute (deg+1)^(dim-1) for (unsigned int d = 1; d < dim; ++d) @@ -160,7 +172,7 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) this->generalized_support_points.resize( GeometryInfo::faces_per_cell * n_face_points + n_interior_points); - this->generalized_face_support_points[0].resize(n_face_points); + this->generalized_face_support_points[face_no].resize(n_face_points); // Number of the point being entered unsigned int current = 0; @@ -173,12 +185,10 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) boundary_weights.reinit(n_face_points, legendre.n()); - // Assert (face_points.size() == this->n_dofs_per_face(), - // ExcInternalError()); - for (unsigned int k = 0; k < n_face_points; ++k) { - this->generalized_face_support_points[0][k] = face_points.point(k); + this->generalized_face_support_points[face_no][k] = + face_points.point(k); // Compute its quadrature // contribution for each // moment. @@ -323,7 +333,8 @@ FE_RaviartThomas::initialize_restriction() for (unsigned int k = 0; k < n_face_points; ++k) for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell(); ++i_child) - for (unsigned int i_face = 0; i_face < this->n_dofs_per_face(); + for (unsigned int i_face = 0; + i_face < this->n_dofs_per_face(face); ++i_face) { // The quadrature @@ -331,13 +342,12 @@ FE_RaviartThomas::initialize_restriction() // subcell are NOT // transformed, so we // have to do it here. - this->restriction[iso][child](face * this->n_dofs_per_face() + - i_face, - i_child) += + this->restriction[iso][child]( + face * this->n_dofs_per_face(face) + i_face, i_child) += Utilities::fixed_power(.5) * q_sub.weight(k) * cached_values_on_face(i_child, k) * this->shape_value_component( - face * this->n_dofs_per_face() + i_face, + face * this->n_dofs_per_face(face) + i_face, q_sub.point(k), GeometryInfo::unit_normal_direction[face]); } @@ -362,9 +372,14 @@ FE_RaviartThomas::initialize_restriction() polynomials[dd] = std::make_unique>(poly); } + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + QGauss q_cell(this->degree); const unsigned int start_cell_dofs = - GeometryInfo::faces_per_cell * this->n_dofs_per_face(); + GeometryInfo::faces_per_cell * this->n_dofs_per_face(face_no); // Store shape values, since the // evaluation suffers if not @@ -512,14 +527,19 @@ FE_RaviartThomas::convert_generalized_support_point_values_to_dof_values( for (unsigned int k = 0; k < n_face_points; ++k) for (unsigned int i = 0; i < boundary_weights.size(1); ++i) { - nodal_values[i + face * this->n_dofs_per_face()] += + nodal_values[i + face * this->n_dofs_per_face(face)] += boundary_weights(k, i) * support_point_values[face * n_face_points + k]( GeometryInfo::unit_normal_direction[face]); } + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + const unsigned int start_cell_dofs = - GeometryInfo::faces_per_cell * this->n_dofs_per_face(); + GeometryInfo::faces_per_cell * this->n_dofs_per_face(face_no); const unsigned int start_cell_points = GeometryInfo::faces_per_cell * n_face_points; diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 90ca1195c6..4ca716558b 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -82,19 +82,27 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal(const unsigned int deg) for (unsigned int i = 0; i < nc; ++i) this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs); } + + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // Fill prolongation matrices with embedding operators FETools::compute_embedding_matrices(*this, this->prolongation); // TODO[TL]: for anisotropic refinement we will probably need a table of // submatrices with an array for each refine case FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face()); + face_embeddings[i].reinit(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0); - this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit((1 << (dim - 1)) * + this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (unsigned int d = 0; d < GeometryInfo::max_children_per_face; ++d) for (unsigned int i = 0; i < face_embeddings[d].m(); ++i) @@ -147,8 +155,14 @@ template void FE_RaviartThomasNodal::initialize_support_points(const unsigned int deg) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + this->generalized_support_points.resize(this->n_dofs_per_cell()); - this->generalized_face_support_points[0].resize(this->n_dofs_per_face()); + this->generalized_face_support_points[face_no].resize( + this->n_dofs_per_face(face_no)); // Number of the point being entered unsigned int current = 0; @@ -162,25 +176,28 @@ FE_RaviartThomasNodal::initialize_support_points(const unsigned int deg) if (dim > 1) { QGauss face_points(deg + 1); - Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError()); - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) - this->generalized_face_support_points[0][k] = face_points.point(k); + Assert(face_points.size() == this->n_dofs_per_face(face_no), + ExcInternalError()); + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) + this->generalized_face_support_points[face_no][k] = + face_points.point(k); Quadrature faces = QProjector::project_to_all_faces(this->reference_cell_type(), face_points); - for (unsigned int k = 0; - k < this->n_dofs_per_face() * GeometryInfo::faces_per_cell; + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) * + GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = faces.point( - k + - QProjector::DataSetDescriptor::face(this->reference_cell_type(), - 0, - true, - false, - false, - this->n_dofs_per_face())); - - current = this->n_dofs_per_face() * GeometryInfo::faces_per_cell; + this->generalized_support_points[k] = + faces.point(k + QProjector::DataSetDescriptor::face( + this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face(face_no))); + + current = + this->n_dofs_per_face(face_no) * GeometryInfo::faces_per_cell; } if (deg == 0) @@ -329,9 +346,9 @@ FE_RaviartThomasNodal:: unsigned int fbase = 0; unsigned int f = 0; for (; f < GeometryInfo::faces_per_cell; - ++f, fbase += this->n_dofs_per_face()) + ++f, fbase += this->n_dofs_per_face(f)) { - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i) { nodal_values[fbase + i] = support_point_values[fbase + i]( GeometryInfo::unit_normal_direction[f]); @@ -463,7 +480,7 @@ template std::vector> FE_RaviartThomasNodal::hp_quad_dof_identities( const FiniteElement &fe_other, - const unsigned int) const + const unsigned int face_no) const { // we can presently only compute // these identities if both FEs are @@ -479,8 +496,10 @@ FE_RaviartThomasNodal::hp_quad_dof_identities( // this works exactly like the line // case above - const unsigned int p = this->n_dofs_per_quad(); - const unsigned int q = fe_q_other->n_dofs_per_quad(); + const unsigned int p = this->n_dofs_per_quad(face_no); + + AssertDimension(fe_q_other->n_unique_faces(), 1); + const unsigned int q = fe_q_other->n_dofs_per_quad(0); std::vector> identities; @@ -575,7 +594,7 @@ void FE_RaviartThomasNodal::get_face_interpolation_matrix( const FiniteElement &x_source_fe, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is only implemented, if the // source FE is also a @@ -585,12 +604,12 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( &x_source_fe) != nullptr), typename FiniteElement::ExcInterpolationNotImplemented()); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // ok, source is a RaviartThomasNodal element, so // we will be able to do the work @@ -608,7 +627,7 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), typename FiniteElement::ExcInterpolationNotImplemented()); // generate a quadrature @@ -618,7 +637,7 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( // which returns the support // points on the face. Quadrature quad_face_support( - source_fe.generalized_face_support_points[0]); + source_fe.generalized_face_support_points[face_no]); // Rule of thumb for FP accuracy, // that can be expected for a @@ -635,11 +654,11 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( quad_face_support, 0); - for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i) { const Point &p = face_projection.point(i); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = this->shape_value_component(this->face_to_cell_index(j, 0), p, 0); @@ -664,11 +683,11 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix( // this point. this must be so // since the shape functions sum up // to 1 - for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j) { double sum = 0.; - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) sum += interpolation_matrix(j, i); Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1), @@ -683,7 +702,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is only implemented, if the // source FE is also a @@ -693,12 +712,12 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( &x_source_fe) != nullptr), typename FiniteElement::ExcInterpolationNotImplemented()); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // ok, source is a RaviartThomasNodal element, so // we will be able to do the work @@ -716,7 +735,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( // lead to problems in the // hp procedures, which use this // method. - Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(), + Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no), typename FiniteElement::ExcInterpolationNotImplemented()); // generate a quadrature @@ -726,7 +745,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( // which returns the support // points on the face. Quadrature quad_face_support( - source_fe.generalized_face_support_points[0]); + source_fe.generalized_face_support_points[face_no]); // Rule of thumb for FP accuracy, // that can be expected for a @@ -745,11 +764,11 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( 0, subface); - for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i) { const Point &p = subface_projection.point(i); - for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j) { double matrix_entry = this->shape_value_component(this->face_to_cell_index(j, 0), p, 0); @@ -774,11 +793,11 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix( // this point. this must be so // since the shape functions sum up // to 1 - for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j) + for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j) { double sum = 0.; - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) sum += interpolation_matrix(j, i); Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1), diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc index 9f0b2025db..145ca3a659 100644 --- a/source/fe/fe_rt_bubbles.cc +++ b/source/fe/fe_rt_bubbles.cc @@ -75,18 +75,26 @@ FE_RT_Bubbles::FE_RT_Bubbles(const unsigned int deg) for (unsigned int i = 0; i < nc; ++i) this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs); } + + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // Fill prolongation matrices with embedding operators // set tolerance to 1, as embedding error accumulate quickly FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.0); FullMatrix face_embeddings[GeometryInfo::max_children_per_face]; for (unsigned int i = 0; i < GeometryInfo::max_children_per_face; ++i) - face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face()); + face_embeddings[i].reinit(this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0); - this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(), - this->n_dofs_per_face()); + this->interface_constraints.reinit((1 << (dim - 1)) * + this->n_dofs_per_face(face_no), + this->n_dofs_per_face(face_no)); unsigned int target_row = 0; for (unsigned int d = 0; d < GeometryInfo::max_children_per_face; ++d) for (unsigned int i = 0; i < face_embeddings[d].m(); ++i) @@ -131,8 +139,14 @@ template void FE_RT_Bubbles::initialize_support_points(const unsigned int deg) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + this->generalized_support_points.resize(this->n_dofs_per_cell()); - this->generalized_face_support_points[0].resize(this->n_dofs_per_face()); + this->generalized_face_support_points[face_no].resize( + this->n_dofs_per_face(face_no)); // Index of the point being entered unsigned int current = 0; @@ -143,25 +157,28 @@ FE_RT_Bubbles::initialize_support_points(const unsigned int deg) if (dim > 1) { QGaussLobatto face_points(deg + 1); - Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError()); - for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k) - this->generalized_face_support_points[0][k] = face_points.point(k); + Assert(face_points.size() == this->n_dofs_per_face(face_no), + ExcInternalError()); + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k) + this->generalized_face_support_points[face_no][k] = + face_points.point(k); Quadrature faces = QProjector::project_to_all_faces(this->reference_cell_type(), face_points); - for (unsigned int k = 0; - k < this->n_dofs_per_face() * GeometryInfo::faces_per_cell; + for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) * + GeometryInfo::faces_per_cell; ++k) - this->generalized_support_points[k] = faces.point( - k + - QProjector::DataSetDescriptor::face(this->reference_cell_type(), - 0, - true, - false, - false, - this->n_dofs_per_face())); - - current = this->n_dofs_per_face() * GeometryInfo::faces_per_cell; + this->generalized_support_points[k] = + faces.point(k + QProjector::DataSetDescriptor::face( + this->reference_cell_type(), + 0, + true, + false, + false, + this->n_dofs_per_face(face_no))); + + current = + this->n_dofs_per_face(face_no) * GeometryInfo::faces_per_cell; } if (deg == 1) @@ -283,9 +300,9 @@ FE_RT_Bubbles::convert_generalized_support_point_values_to_dof_values( unsigned int fbase = 0; unsigned int f = 0; for (; f < GeometryInfo::faces_per_cell; - ++f, fbase += this->n_dofs_per_face()) + ++f, fbase += this->n_dofs_per_face(f)) { - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i) { nodal_values[fbase + i] = support_point_values[fbase + i]( GeometryInfo::unit_normal_direction[f]); diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index 78d40b2ef0..de8169987e 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -858,7 +858,7 @@ FESystem::face_to_cell_index(const unsigned int face_dof_index, // the DoFs within their own numbering. thus, translate to // the base element numbering and then back const std::pair, unsigned int> - face_base_index = this->face_system_to_base_index(face_dof_index); + face_base_index = this->face_system_to_base_index(face_dof_index, face); const unsigned int base_face_to_cell_index = this->base_element(face_base_index.first.first) @@ -1392,6 +1392,11 @@ template void FESystem::build_interface_constraints() { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(this->n_unique_faces(), 1); + const unsigned int face_no = 0; + // check whether all base elements implement their interface constraint // matrices. if this is not the case, then leave the interface costraints of // this composed element empty as well; however, the rest of the element is @@ -1422,7 +1427,7 @@ FESystem::build_interface_constraints() // data type, first value in pair is (base element,instance of base // element), second is index within this instance const std::pair, unsigned int> - n_index = this->face_system_to_base_table[0][n]; + n_index = this->face_system_to_base_table[face_no][n]; // likewise for the m index. this is more complicated due to the // strange ordering we have for the dofs on the refined faces. @@ -1461,7 +1466,7 @@ FESystem::build_interface_constraints() const unsigned int tmp1 = 2 * this->n_dofs_per_vertex() + index_in_line; m_index.first = - this->face_system_to_base_table[0][tmp1].first; + this->face_system_to_base_table[face_no][tmp1].first; // what we are still missing is the index of m within the // base elements interface_constraints table @@ -1472,12 +1477,13 @@ FESystem::build_interface_constraints() // dof, we can construct the rest: tmp2 will denote the // index of this shape function among the line shape // functions: - Assert(this->face_system_to_base_table[0][tmp1].second >= - 2 * base_element(m_index.first.first) - .n_dofs_per_vertex(), - ExcInternalError()); + Assert( + this->face_system_to_base_table[face_no][tmp1].second >= + 2 * + base_element(m_index.first.first).n_dofs_per_vertex(), + ExcInternalError()); const unsigned int tmp2 = - this->face_system_to_base_table[0][tmp1].second - + this->face_system_to_base_table[face_no][tmp1].second - 2 * base_element(m_index.first.first).n_dofs_per_vertex(); Assert(tmp2 < base_element(m_index.first.first) .n_dofs_per_line(), @@ -1517,14 +1523,15 @@ FESystem::build_interface_constraints() const unsigned int tmp1 = 4 * this->n_dofs_per_vertex() + index_in_line; m_index.first = - this->face_system_to_base_table[0][tmp1].first; + this->face_system_to_base_table[face_no][tmp1].first; - Assert(this->face_system_to_base_table[0][tmp1].second >= - 4 * base_element(m_index.first.first) - .n_dofs_per_vertex(), - ExcInternalError()); + Assert( + this->face_system_to_base_table[face_no][tmp1].second >= + 4 * + base_element(m_index.first.first).n_dofs_per_vertex(), + ExcInternalError()); const unsigned int tmp2 = - this->face_system_to_base_table[0][tmp1].second - + this->face_system_to_base_table[face_no][tmp1].second - 4 * base_element(m_index.first.first).n_dofs_per_vertex(); Assert(tmp2 < base_element(m_index.first.first) .n_dofs_per_line(), @@ -1543,42 +1550,45 @@ FESystem::build_interface_constraints() const unsigned int index_in_quad = (m - 5 * this->n_dofs_per_vertex() - 12 * this->n_dofs_per_line()) % - this->n_dofs_per_quad(); - Assert(index_in_quad < this->n_dofs_per_quad(), + this->n_dofs_per_quad(face_no); + Assert(index_in_quad < this->n_dofs_per_quad(face_no), ExcInternalError()); const unsigned int sub_quad = ((m - 5 * this->n_dofs_per_vertex() - 12 * this->n_dofs_per_line()) / - this->n_dofs_per_quad()); + this->n_dofs_per_quad(face_no)); Assert(sub_quad < 4, ExcInternalError()); const unsigned int tmp1 = 4 * this->n_dofs_per_vertex() + 4 * this->n_dofs_per_line() + index_in_quad; - Assert(tmp1 < this->face_system_to_base_table[0].size(), + Assert(tmp1 < + this->face_system_to_base_table[face_no].size(), ExcInternalError()); m_index.first = - this->face_system_to_base_table[0][tmp1].first; - - Assert(this->face_system_to_base_table[0][tmp1].second >= - 4 * base_element(m_index.first.first) - .n_dofs_per_vertex() + - 4 * base_element(m_index.first.first) - .n_dofs_per_line(), - ExcInternalError()); + this->face_system_to_base_table[face_no][tmp1].first; + + Assert( + this->face_system_to_base_table[face_no][tmp1].second >= + 4 * base_element(m_index.first.first) + .n_dofs_per_vertex() + + 4 * + base_element(m_index.first.first).n_dofs_per_line(), + ExcInternalError()); const unsigned int tmp2 = - this->face_system_to_base_table[0][tmp1].second - + this->face_system_to_base_table[face_no][tmp1].second - 4 * base_element(m_index.first.first).n_dofs_per_vertex() - 4 * base_element(m_index.first.first).n_dofs_per_line(); Assert(tmp2 < base_element(m_index.first.first) - .n_dofs_per_quad(), + .n_dofs_per_quad(face_no), ExcInternalError()); m_index.second = 5 * base_element(m_index.first.first).n_dofs_per_vertex() + 12 * base_element(m_index.first.first).n_dofs_per_line() + - base_element(m_index.first.first).n_dofs_per_quad() * + base_element(m_index.first.first) + .n_dofs_per_quad(face_no) * sub_quad + tmp2; } @@ -1648,18 +1658,25 @@ FESystem::initialize( // FiniteElement this->system_to_component_table.resize(this->n_dofs_per_cell()); - this->face_system_to_component_table.resize(1); - this->face_system_to_component_table[0].resize(this->n_dofs_per_face()); - FETools::Compositing::build_cell_tables(this->system_to_base_table, this->system_to_component_table, this->component_to_base_table, *this); - FETools::Compositing::build_face_tables( - this->face_system_to_base_table[0], - this->face_system_to_component_table[0], - *this); + this->face_system_to_component_table.resize(this->n_unique_faces()); + + for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no) + { + this->face_system_to_component_table[0].resize( + this->n_dofs_per_face(face_no)); + + FETools::Compositing::build_face_tables( + this->face_system_to_base_table[face_no], + this->face_system_to_component_table[face_no], + *this, + true, + face_no); + } } // now initialize interface constraints, support points, and other tables. @@ -1702,43 +1719,56 @@ FESystem::initialize( // initialize face support points (for dim==2,3). same procedure as above if (dim > 1) init_tasks += Threads::new_task([&]() { - // if one of the base elements has no support points, then it makes no - // sense to define support points for the composed element. In that case, - // return an empty array to demonstrate that fact (note that we ask - // whether the base element has no support points at all, not only none on - // the face!) - // - // on the other hand, if there is an element that simply has no degrees of - // freedom on the face at all, then we don't care whether it has support - // points or not. this is, for example, the case for the stable Stokes - // element Q(p)^dim \times DGP(p-1). - for (unsigned int base_el = 0; base_el < this->n_base_elements(); - ++base_el) - if (!base_element(base_el).has_support_points() && - (base_element(base_el).n_dofs_per_face() > 0)) - { - this->unit_face_support_points[0].resize(0); - return; - } - + for (unsigned int face_no = 0; face_no < this->n_unique_faces(); + ++face_no) + { + // if one of the base elements has no support points, then it makes + // no sense to define support points for the composed element. In + // that case, return an empty array to demonstrate that fact (note + // that we ask whether the base element has no support points at + // all, not only none on the face!) + // + // on the other hand, if there is an element that simply has no + // degrees of freedom on the face at all, then we don't care whether + // it has support points or not. this is, for example, the case for + // the stable Stokes element Q(p)^dim \times DGP(p-1). + bool flag_has_no_support_points = false; + + for (unsigned int base_el = 0; base_el < this->n_base_elements(); + ++base_el) + if (!base_element(base_el).has_support_points() && + (base_element(base_el).n_dofs_per_face(face_no) > 0)) + { + this->unit_face_support_points[face_no].resize(0); + flag_has_no_support_points = true; + break; + } - // generate unit face support points from unit support points of sub - // elements - this->unit_face_support_points[0].resize(this->n_dofs_per_face()); - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) - { - const unsigned int base_i = - this->face_system_to_base_table[0][i].first.first; - const unsigned int index_in_base = - this->face_system_to_base_table[0][i].second; + if (flag_has_no_support_points) + continue; - Assert(index_in_base < - base_element(base_i).unit_face_support_points[0].size(), - ExcInternalError()); + // generate unit face support points from unit support points of sub + // elements + this->unit_face_support_points[face_no].resize( + this->n_dofs_per_face(face_no)); - this->unit_face_support_points[0][i] = - base_element(base_i).unit_face_support_points[0][index_in_base]; + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) + { + const unsigned int base_i = + this->face_system_to_base_table[face_no][i].first.first; + const unsigned int index_in_base = + this->face_system_to_base_table[face_no][i].second; + + Assert( + index_in_base < + base_element(base_i).unit_face_support_points[face_no].size(), + ExcInternalError()); + + this->unit_face_support_points[face_no][i] = + base_element(base_i) + .unit_face_support_points[face_no][index_in_base]; + } } }); @@ -1813,36 +1843,40 @@ FESystem::initialize( // initialize quad dof index permutation in 3d and higher if (dim >= 3) init_tasks += Threads::new_task([&]() { - // the array into which we want to write should have the correct size - // already. - Assert(this->adjust_quad_dof_index_for_face_orientation_table[0] - .n_elements() == 8 * this->n_dofs_per_quad(), - ExcInternalError()); - - // to obtain the shifts for this composed element, copy the shift - // information of the base elements - unsigned int index = 0; - for (unsigned int b = 0; b < this->n_base_elements(); ++b) + for (unsigned int face_no = 0; face_no < this->n_unique_faces(); + ++face_no) { - const Table<2, int> &temp = - this->base_element(b) - .adjust_quad_dof_index_for_face_orientation_table[0]; - for (unsigned int c = 0; c < this->element_multiplicity(b); ++c) + // the array into which we want to write should have the correct size + // already. + Assert(this->adjust_quad_dof_index_for_face_orientation_table[face_no] + .n_elements() == 8 * this->n_dofs_per_quad(face_no), + ExcInternalError()); + + // to obtain the shifts for this composed element, copy the shift + // information of the base elements + unsigned int index = 0; + for (unsigned int b = 0; b < this->n_base_elements(); ++b) { - for (unsigned int i = 0; i < temp.size(0); ++i) - for (unsigned int j = 0; j < 8; ++j) - this->adjust_quad_dof_index_for_face_orientation_table[0]( - index + i, j) = temp(i, j); - index += temp.size(0); + const Table<2, int> &temp = + this->base_element(b) + .adjust_quad_dof_index_for_face_orientation_table[face_no]; + for (unsigned int c = 0; c < this->element_multiplicity(b); ++c) + { + for (unsigned int i = 0; i < temp.size(0); ++i) + for (unsigned int j = 0; j < 8; ++j) + this->adjust_quad_dof_index_for_face_orientation_table + [face_no](index + i, j) = temp(i, j); + index += temp.size(0); + } } + Assert(index == this->n_dofs_per_quad(face_no), ExcInternalError()); } - Assert(index == this->n_dofs_per_quad(), ExcInternalError()); // additionally compose the permutation information for lines Assert(this->adjust_line_dof_index_for_line_orientation_table.size() == this->n_dofs_per_line(), ExcInternalError()); - index = 0; + unsigned int index = 0; for (unsigned int b = 0; b < this->n_base_elements(); ++b) { const std::vector &temp2 = @@ -1887,12 +1921,12 @@ FESystem::get_face_interpolation_matrix( FullMatrix & interpolation_matrix, const unsigned int face_no) const { - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // since dofs for each base are independent, we only have to stack things up // from base element to base element @@ -1928,8 +1962,8 @@ FESystem::get_face_interpolation_matrix( ExcNotImplemented()); // get the interpolation from the bases - base_to_base_interpolation.reinit(base_other.n_dofs_per_face(), - base.n_dofs_per_face()); + base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no), + base.n_dofs_per_face(face_no)); base.get_face_interpolation_matrix(base_other, base_to_base_interpolation, face_no); @@ -1937,16 +1971,19 @@ FESystem::get_face_interpolation_matrix( // now translate entries. we'd like to have something like // face_base_to_system_index, but that doesn't exist. rather, all we // have is the reverse. well, use that then - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) - if (this->face_system_to_base_index(i).first == + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) + if (this->face_system_to_base_index(i, face_no).first == std::make_pair(base_index, multiplicity)) - for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face(); + for (unsigned int j = 0; + j < fe_other_system->n_dofs_per_face(face_no); ++j) - if (fe_other_system->face_system_to_base_index(j).first == + if (fe_other_system->face_system_to_base_index(j, face_no) + .first == std::make_pair(base_index_other, multiplicity_other)) interpolation_matrix(j, i) = base_to_base_interpolation( - fe_other_system->face_system_to_base_index(j).second, - this->face_system_to_base_index(i).second); + fe_other_system->face_system_to_base_index(j, face_no) + .second, + this->face_system_to_base_index(i, face_no).second); // advance to the next base element for this and the other fe_system; // see if we can simply advance the multiplicity by one, or if have to @@ -2006,12 +2043,12 @@ FESystem::get_subface_interpolation_matrix( (dynamic_cast *>(&x_source_fe) != nullptr), (typename FiniteElement::ExcInterpolationNotImplemented())); - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // since dofs for each base are independent, we only have to stack things up // from base element to base element @@ -2048,8 +2085,8 @@ FESystem::get_subface_interpolation_matrix( ExcNotImplemented()); // get the interpolation from the bases - base_to_base_interpolation.reinit(base_other.n_dofs_per_face(), - base.n_dofs_per_face()); + base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no), + base.n_dofs_per_face(face_no)); base.get_subface_interpolation_matrix(base_other, subface, base_to_base_interpolation, @@ -2058,16 +2095,19 @@ FESystem::get_subface_interpolation_matrix( // now translate entries. we'd like to have something like // face_base_to_system_index, but that doesn't exist. rather, all we // have is the reverse. well, use that then - for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i) - if (this->face_system_to_base_index(i).first == + for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i) + if (this->face_system_to_base_index(i, face_no).first == std::make_pair(base_index, multiplicity)) - for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face(); + for (unsigned int j = 0; + j < fe_other_system->n_dofs_per_face(face_no); ++j) - if (fe_other_system->face_system_to_base_index(j).first == + if (fe_other_system->face_system_to_base_index(j, face_no) + .first == std::make_pair(base_index_other, multiplicity_other)) interpolation_matrix(j, i) = base_to_base_interpolation( - fe_other_system->face_system_to_base_index(j).second, - this->face_system_to_base_index(i).second); + fe_other_system->face_system_to_base_index(j, face_no) + .second, + this->face_system_to_base_index(i, face_no).second); // advance to the next base element for this and the other fe_system; // see if we can simply advance the multiplicity by one, or if have to @@ -2117,7 +2157,8 @@ template template std::vector> FESystem::hp_object_dof_identities( - const FiniteElement &fe_other) const + const FiniteElement &fe_other, + const unsigned int face_no) const { // since dofs on each subobject (vertex, line, ...) are ordered such that // first come all from the first base element all multiplicities, then @@ -2167,7 +2208,8 @@ FESystem::hp_object_dof_identities( base_identities = base.hp_line_dof_identities(base_other); break; case 2: - base_identities = base.hp_quad_dof_identities(base_other); + base_identities = + base.hp_quad_dof_identities(base_other, face_no); break; default: Assert(false, ExcNotImplemented()); @@ -2247,9 +2289,9 @@ template std::vector> FESystem::hp_quad_dof_identities( const FiniteElement &fe_other, - const unsigned int) const + const unsigned int face_no) const { - return hp_object_dof_identities<2>(fe_other); + return hp_object_dof_identities<2>(fe_other, face_no); } @@ -2354,22 +2396,27 @@ FESystem::unit_face_support_point( const unsigned int index, const unsigned int face_no) const { - AssertIndexRange(index, this->n_dofs_per_face()); - Assert((this->unit_face_support_points[0].size() == - this->n_dofs_per_face()) || - (this->unit_face_support_points[0].size() == 0), - (typename FiniteElement::ExcFEHasNoSupportPoints())); + AssertIndexRange(index, this->n_dofs_per_face(face_no)); + Assert( + (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() == this->n_dofs_per_face(face_no)) || + (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() == 0), + (typename FiniteElement::ExcFEHasNoSupportPoints())); // let's see whether we have the information pre-computed - if (this->unit_face_support_points[0].size() != 0) - return this->unit_face_support_points[0][index]; + if (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + .size() != 0) + return this + ->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no] + [index]; else // no. ask the base element whether it would like to provide this // information return ( - base_element(this->face_system_to_base_index(index).first.first) - .unit_face_support_point(this->face_system_to_base_index(index).second, - face_no)); + base_element(this->face_system_to_base_index(index, face_no).first.first) + .unit_face_support_point( + this->face_system_to_base_index(index, face_no).second, face_no)); } diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index c77d24dc53..981a5bc527 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -102,7 +102,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) std::vector> &face_system_to_component_table, const FiniteElement &fe, - const bool do_tensor_product); + const bool do_tensor_product, + const unsigned int face_no); \} diff --git a/source/fe/fe_trace.cc b/source/fe/fe_trace.cc index 5190e0e655..584465e413 100644 --- a/source/fe/fe_trace.cc +++ b/source/fe/fe_trace.cc @@ -53,7 +53,10 @@ FE_TraceQ::FE_TraceQ(const unsigned int degree) FETools::hierarchic_to_lexicographic_numbering(degree)); // Initialize face support points - this->unit_face_support_points[0] = fe_q.get_unit_face_support_points(); + AssertDimension(this->n_unique_faces(), fe_q.n_unique_faces()); + for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no) + this->unit_face_support_points[face_no] = + fe_q.get_unit_face_support_points(face_no); // initialize unit support points (this makes it possible to assign initial // values to FE_TraceQ). Note that we simply take the points of fe_q but @@ -231,15 +234,15 @@ FE_TraceQ::get_subface_interpolation_matrix( const FiniteElement &x_source_fe, const unsigned int subface, FullMatrix & interpolation_matrix, - const unsigned int) const + const unsigned int face_no) const { // this is the code from FE_FaceQ - Assert(interpolation_matrix.n() == this->n_dofs_per_face(), + Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.n(), - this->n_dofs_per_face())); - Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(), + this->n_dofs_per_face(face_no))); + Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no), ExcDimensionMismatch(interpolation_matrix.m(), - x_source_fe.n_dofs_per_face())); + x_source_fe.n_dofs_per_face(face_no))); // see if source is a FaceQ element if (const FE_TraceQ *source_fe = @@ -247,7 +250,8 @@ FE_TraceQ::get_subface_interpolation_matrix( { fe_q.get_subface_interpolation_matrix(source_fe->fe_q, subface, - interpolation_matrix); + interpolation_matrix, + face_no); } else if (dynamic_cast *>(&x_source_fe) != nullptr) {