From 88cd461ff35f5a81c261c1af903d285db1a3e1bd Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 31 Jul 2020 23:55:57 +0200 Subject: [PATCH] Convert all face-related data structures of FiniteElement to vectors --- include/deal.II/fe/fe.h | 19 +++++---- source/fe/fe.cc | 45 ++++++++++++-------- source/fe/fe_abf.cc | 4 +- source/fe/fe_bdm.cc | 4 +- source/fe/fe_enriched.cc | 10 ++--- source/fe/fe_face.cc | 18 ++++---- source/fe/fe_nedelec.cc | 49 +++++++++++----------- source/fe/fe_p1nc.cc | 6 +-- source/fe/fe_q_base.cc | 26 ++++++------ source/fe/fe_q_hierarchical.cc | 4 +- source/fe/fe_raviart_thomas.cc | 4 +- source/fe/fe_raviart_thomas_nodal.cc | 8 ++-- source/fe/fe_rt_bubbles.cc | 4 +- source/fe/fe_system.cc | 62 +++++++++++++++------------- source/fe/fe_trace.cc | 2 +- source/simplex/fe_lib.cc | 28 ++++++------- 16 files changed, 156 insertions(+), 137 deletions(-) diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 586a3784ff..8a24356dbc 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -2436,7 +2436,7 @@ protected: * get_unit_face_support_points() function for a discussion of what * contributes a face support point. */ - std::vector> unit_face_support_points; + std::vector>> unit_face_support_points; /** * Support points used for interpolation functions of non-Lagrangian @@ -2448,7 +2448,7 @@ protected: * Face support points used for interpolation functions of non-Lagrangian * elements. */ - std::vector> generalized_face_support_points; + std::vector>> generalized_face_support_points; /** * For faces with non-standard face_orientation in 3D, the dofs on faces @@ -2465,7 +2465,7 @@ protected: * no permutation at all. Derived finite element classes have to * fill this Table with the correct values. */ - Table<2, int> adjust_quad_dof_index_for_face_orientation_table; + std::vector> adjust_quad_dof_index_for_face_orientation_table; /** * For lines with non-standard line_orientation in 3D, the dofs on lines @@ -2496,7 +2496,7 @@ protected: * information thus makes only sense if a shape function is non-zero in only * one component. */ - std::vector> + std::vector>> face_system_to_component_table; /** @@ -2521,7 +2521,8 @@ protected: /** * Likewise for the indices on faces. */ - std::vector, unsigned int>> + std::vector< + std::vector, unsigned int>>> face_system_to_base_table; /** @@ -3122,7 +3123,7 @@ inline std::pair FiniteElement::face_system_to_component_index( const unsigned int index) const { - AssertIndexRange(index, face_system_to_component_table.size()); + AssertIndexRange(index, face_system_to_component_table[0].size()); // in debug mode, check whether the // function is primitive, since @@ -3141,7 +3142,7 @@ FiniteElement::face_system_to_component_index( (typename FiniteElement::ExcShapeFunctionNotPrimitive( index))); - return face_system_to_component_table[index]; + return face_system_to_component_table[0][index]; } @@ -3162,8 +3163,8 @@ inline std::pair, unsigned int> FiniteElement::face_system_to_base_index( const unsigned int index) const { - AssertIndexRange(index, face_system_to_base_table.size()); - return face_system_to_base_table[index]; + AssertIndexRange(index, face_system_to_base_table[0].size()); + return face_system_to_base_table[0][index]; } diff --git a/source/fe/fe.cc b/source/fe/fe.cc index bb90aaf84c..e46a6ef4c7 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -59,14 +59,9 @@ FiniteElement::FiniteElement( const std::vector & r_i_a_f, const std::vector &nonzero_c) : FiniteElementData(fe_data) - , adjust_quad_dof_index_for_face_orientation_table(dim == 3 ? - this->n_dofs_per_quad() : - 0, - dim == 3 ? 8 : 0) , adjust_line_dof_index_for_line_orientation_table( dim == 3 ? this->n_dofs_per_line() : 0) , system_to_base_table(this->n_dofs_per_cell()) - , face_system_to_base_table(this->n_dofs_per_face()) , component_to_base_table(this->components, std::make_pair(std::make_pair(0U, 0U), 0U)) , @@ -109,17 +104,23 @@ FiniteElement::FiniteElement( if (this->is_primitive()) { system_to_component_table.resize(this->n_dofs_per_cell()); - face_system_to_component_table.resize(this->n_dofs_per_face()); 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[j] = std::pair(0, j); + face_system_to_component_table[0][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[j] = std::make_pair(std::make_pair(0U, 0U), j); + face_system_to_base_table[0][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); @@ -140,7 +141,17 @@ FiniteElement::FiniteElement( FullMatrix()); } - adjust_quad_dof_index_for_face_orientation_table.fill(0); + if (dim == 3) + { + adjust_quad_dof_index_for_face_orientation_table.resize(1); + + 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); + } + + unit_face_support_points.resize(1); + generalized_face_support_points.resize(1); } @@ -653,10 +664,10 @@ FiniteElement::adjust_quad_dof_index_for_face_orientation( // the function should also not have been // called AssertIndexRange(index, this->n_dofs_per_quad()); - Assert(adjust_quad_dof_index_for_face_orientation_table.n_elements() == + Assert(adjust_quad_dof_index_for_face_orientation_table[0].n_elements() == 8 * this->n_dofs_per_quad(), ExcInternalError()); - return index + adjust_quad_dof_index_for_face_orientation_table( + return index + adjust_quad_dof_index_for_face_orientation_table[0]( index, 4 * face_orientation + 2 * face_flip + face_rotation); } @@ -1057,10 +1068,10 @@ FiniteElement::get_unit_face_support_points() const // support points, but only if // there are as many as there are // degrees of freedom on a face - Assert((unit_face_support_points.size() == 0) || - (unit_face_support_points.size() == this->n_dofs_per_face()), + Assert((unit_face_support_points[0].size() == 0) || + (unit_face_support_points[0].size() == this->n_dofs_per_face()), ExcInternalError()); - return unit_face_support_points; + return unit_face_support_points[0]; } @@ -1069,7 +1080,7 @@ template bool FiniteElement::has_face_support_points() const { - return (unit_face_support_points.size() != 0); + return (unit_face_support_points[0].size() != 0); } @@ -1080,9 +1091,9 @@ FiniteElement::unit_face_support_point( const unsigned int index) const { AssertIndexRange(index, this->n_dofs_per_face()); - Assert(unit_face_support_points.size() == this->n_dofs_per_face(), + Assert(unit_face_support_points[0].size() == this->n_dofs_per_face(), ExcFEHasNoSupportPoints()); - return unit_face_support_points[index]; + return unit_face_support_points[0][index]; } diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc index 83732630e5..c55552234b 100644 --- a/source/fe/fe_abf.cc +++ b/source/fe/fe_abf.cc @@ -160,7 +160,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.resize(n_face_points); + this->generalized_face_support_points[0].resize(n_face_points); // These might be required when the faces contribution is computed @@ -194,7 +194,7 @@ FE_ABF::initialize_support_points(const unsigned int deg) for (unsigned int k = 0; k < n_face_points; ++k) { - this->generalized_face_support_points[k] = face_points.point(k); + this->generalized_face_support_points[0][k] = face_points.point(k); // Compute its quadrature // contribution for each // moment. diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc index a85ede4162..dd1e4b469e 100644 --- a/source/fe/fe_bdm.cc +++ b/source/fe/fe_bdm.cc @@ -333,9 +333,9 @@ FE_BDM::initialize_support_points(const unsigned int deg) QGauss face_points(deg + 1); // Copy the quadrature formula to the face points. - this->generalized_face_support_points.resize(face_points.size()); + this->generalized_face_support_points[0].resize(face_points.size()); for (unsigned int k = 0; k < face_points.size(); ++k) - this->generalized_face_support_points[k] = face_points.point(k); + this->generalized_face_support_points[0][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 diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc index 3a37f8fd2b..012f4af037 100644 --- a/source/fe/fe_enriched.cc +++ b/source/fe/fe_enriched.cc @@ -440,7 +440,7 @@ 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.resize(this->n_dofs_per_face()); + 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, @@ -449,8 +449,8 @@ FE_Enriched::initialize( false); FETools::Compositing::build_face_tables( - this->face_system_to_base_table, - this->face_system_to_component_table, + this->face_system_to_base_table[0], + this->face_system_to_component_table[0], *this, false); } @@ -468,8 +468,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 = fe_system->unit_face_support_points; + this->unit_support_points = fe_system->unit_support_points; + this->unit_face_support_points[0] = fe_system->unit_face_support_points[0]; } // 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 8eefdcbbb6..fb4343c0f4 100644 --- a/source/fe/fe_face.cc +++ b/source/fe/fe_face.cc @@ -61,12 +61,12 @@ FE_FaceQ::FE_FaceQ(const unsigned int degree) { // initialize unit face support points const unsigned int codim = dim - 1; - this->unit_face_support_points.resize( + this->unit_face_support_points[0].resize( Utilities::fixed_power(this->degree + 1)); if (this->degree == 0) for (unsigned int d = 0; d < codim; ++d) - this->unit_face_support_points[0][d] = 0.5; + this->unit_face_support_points[0][0][d] = 0.5; else { std::vector> points = @@ -85,16 +85,16 @@ FE_FaceQ::FE_FaceQ(const unsigned int degree) if (codim > 2) p(2) = points[iz][0]; - this->unit_face_support_points[k++] = p; + this->unit_face_support_points[0][k++] = p; } - AssertDimension(k, this->unit_face_support_points.size()); + AssertDimension(k, this->unit_face_support_points[0].size()); } // initialize unit support points (this makes it possible to assign initial // values to FE_FaceQ) this->unit_support_points.resize(GeometryInfo::faces_per_cell * - this->unit_face_support_points.size()); - const unsigned int n_face_dofs = this->unit_face_support_points.size(); + this->unit_face_support_points[0].size()); + const unsigned int n_face_dofs = this->unit_face_support_points[0].size(); for (unsigned int i = 0; i < n_face_dofs; ++i) for (unsigned int d = 0; d < dim; ++d) { @@ -106,9 +106,9 @@ FE_FaceQ::FE_FaceQ(const unsigned int degree) if (dim == 3 && d == 1) renumber = i / (degree + 1) + (degree + 1) * (i % (degree + 1)); this->unit_support_points[n_face_dofs * 2 * d + i][e] = - this->unit_face_support_points[renumber][c]; + this->unit_face_support_points[0][renumber][c]; this->unit_support_points[n_face_dofs * (2 * d + 1) + i][e] = - this->unit_face_support_points[renumber][c]; + this->unit_face_support_points[0][renumber][c]; this->unit_support_points[n_face_dofs * (2 * d + 1) + i][d] = 1; ++c; } @@ -514,7 +514,7 @@ FE_FaceQ<1, spacedim>::FE_FaceQ(const unsigned int degree) std::vector(1, true), std::vector(1, ComponentMask(1, true))) { - this->unit_face_support_points.resize(1); + this->unit_face_support_points[0].resize(1); // initialize unit support points (this makes it possible to assign initial // values to FE_FaceQ) diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc index 9645e6a0c9..f5da3d679c 100644 --- a/source/fe/fe_nedelec.cc +++ b/source/fe/fe_nedelec.cc @@ -273,11 +273,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.resize(n_edge_points); + this->generalized_face_support_points[0].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[q_point] = + this->generalized_face_support_points[0][q_point] = reference_edge_quadrature.point(q_point); if (order > 0) @@ -310,7 +310,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[q_point](0)); + this->generalized_face_support_points[0][q_point](0)); } for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point) @@ -380,8 +380,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.resize(4 * n_edge_points + - n_face_points); + this->generalized_face_support_points[0].resize(4 * n_edge_points + + n_face_points); this->generalized_support_points.resize(n_boundary_points + n_interior_points); @@ -391,8 +391,8 @@ 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[line * n_edge_points + - q_point] = + this->generalized_face_support_points[0][line * n_edge_points + + q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( ReferenceCell::get_hypercube(dim - 1), @@ -421,13 +421,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[q_point](1)); + this->generalized_face_support_points[0][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[q_point + 4 * n_edge_points] = + this + ->generalized_face_support_points[0][q_point + 4 * n_edge_points] = reference_face_quadrature.point(q_point); for (unsigned int i = 0; i <= order; ++i) @@ -436,23 +437,23 @@ 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[q_point + - 4 * n_edge_points]( + this->generalized_face_support_points[0][q_point + + 4 * n_edge_points]( 0)) * lobatto_polynomials[j + 2].value( - this->generalized_face_support_points[q_point + - 4 * n_edge_points]( + this->generalized_face_support_points[0][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[q_point + - 4 * n_edge_points]( + this->generalized_face_support_points[0][q_point + + 4 * n_edge_points]( 1)) * lobatto_polynomials[j + 2].value( - this->generalized_face_support_points[q_point + - 4 * n_edge_points]( + this->generalized_face_support_points[0][q_point + + 4 * n_edge_points]( 0)); } } @@ -485,7 +486,7 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order) else { - this->generalized_face_support_points.resize(4 * n_edge_points); + this->generalized_face_support_points[0].resize(4 * n_edge_points); this->generalized_support_points.resize( GeometryInfo::lines_per_cell * n_edge_points); @@ -494,8 +495,8 @@ 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[line * n_edge_points + - q_point] = + this->generalized_face_support_points[0][line * n_edge_points + + q_point] = edge_quadrature.point( QProjector::DataSetDescriptor::face( ReferenceCell::get_hypercube(dim - 1), @@ -3182,7 +3183,7 @@ 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[q_point](0)); + this->generalized_face_support_points[0][q_point](0)); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3461,7 +3462,7 @@ 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[q_point](1)); + this->generalized_face_support_points[0][q_point](1)); FullMatrix system_matrix_inv(this->degree - 1, this->degree - 1); @@ -3540,10 +3541,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 - [q_point + 4 * n_edge_points](0)) * + [0][q_point + 4 * n_edge_points](0)) * lobatto_polynomials[j + 2].value( this->generalized_face_support_points - [q_point + 4 * n_edge_points](1)); + [0][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_p1nc.cc b/source/fe/fe_p1nc.cc index e86950d1c7..8d40cbce7c 100644 --- a/source/fe/fe_p1nc.cc +++ b/source/fe/fe_p1nc.cc @@ -28,9 +28,9 @@ FE_P1NC::FE_P1NC() std::vector(4, ComponentMask(1, true))) { // face support points: 2 end vertices - unit_face_support_points.resize(2); - unit_face_support_points[0][0] = 0.0; - unit_face_support_points[1][0] = 1.0; + unit_face_support_points[0].resize(2); + unit_face_support_points[0][0][0] = 0.0; + unit_face_support_points[0][1][0] = 1.0; // initialize constraints matrix initialize_constraints(); diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index d6664d6ea9..9528d3958e 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -942,7 +942,7 @@ FE_Q_Base::initialize_unit_face_support_points( if (dim == 1) return; - this->unit_face_support_points.resize( + this->unit_face_support_points[0].resize( Utilities::fixed_power(q_degree + 1)); // find renumbering of faces and assign from values of quadrature @@ -958,9 +958,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.resize(support_quadrature.size()); + this->unit_face_support_points[0].resize(support_quadrature.size()); for (unsigned int k = 0; k < support_quadrature.size(); ++k) - this->unit_face_support_points[face_index_map[k]] = + this->unit_face_support_points[0][face_index_map[k]] = support_quadrature.point(k); } @@ -975,8 +975,8 @@ FE_Q_Base:: if (dim < 3) return; - Assert(this->adjust_quad_dof_index_for_face_orientation_table.n_elements() == - 8 * this->n_dofs_per_quad(), + Assert(this->adjust_quad_dof_index_for_face_orientation_table[0] + .n_elements() == 8 * this->n_dofs_per_quad(), ExcInternalError()); const unsigned int n = q_degree - 1; @@ -1007,27 +1007,27 @@ FE_Q_Base:: 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(local, 0) = + this->adjust_quad_dof_index_for_face_orientation_table[0](local, 0) = j + i * n - local; // face_orientation=false, face_flip=false, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table(local, 1) = + this->adjust_quad_dof_index_for_face_orientation_table[0](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(local, 2) = + this->adjust_quad_dof_index_for_face_orientation_table[0](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(local, 3) = + this->adjust_quad_dof_index_for_face_orientation_table[0](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(local, 4) = 0; + this->adjust_quad_dof_index_for_face_orientation_table[0](local, 4) = 0; // face_orientation=true, face_flip=false, face_rotation=true - this->adjust_quad_dof_index_for_face_orientation_table(local, 5) = + this->adjust_quad_dof_index_for_face_orientation_table[0](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(local, 6) = + this->adjust_quad_dof_index_for_face_orientation_table[0](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(local, 7) = + this->adjust_quad_dof_index_for_face_orientation_table[0](local, 7) = (n - 1 - j) + i * n - local; } diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc index 964aae3ba3..099c58900e 100644 --- a/source/fe/fe_q_hierarchical.cc +++ b/source/fe/fe_q_hierarchical.cc @@ -1903,7 +1903,7 @@ FE_Q_Hierarchical::initialize_generalized_face_support_points() for (unsigned int i = 1; i < codim; ++i) n *= this->degree + 1; - this->generalized_face_support_points.resize(n); + this->generalized_face_support_points[0].resize(n); Point p; @@ -1936,7 +1936,7 @@ FE_Q_Hierarchical::initialize_generalized_face_support_points() else p(2) = .5; } - this->generalized_face_support_points[face_renumber[k++]] = p; + this->generalized_face_support_points[0][face_renumber[k++]] = p; } } diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc index 310806e5bc..5b4334ae6c 100644 --- a/source/fe/fe_raviart_thomas.cc +++ b/source/fe/fe_raviart_thomas.cc @@ -160,7 +160,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.resize(n_face_points); + this->generalized_face_support_points[0].resize(n_face_points); // Number of the point being entered unsigned int current = 0; @@ -178,7 +178,7 @@ FE_RaviartThomas::initialize_support_points(const unsigned int deg) for (unsigned int k = 0; k < n_face_points; ++k) { - this->generalized_face_support_points[k] = face_points.point(k); + this->generalized_face_support_points[0][k] = face_points.point(k); // Compute its quadrature // contribution for each // moment. diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 0ddeb83502..f633cc23b2 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -148,7 +148,7 @@ void FE_RaviartThomasNodal::initialize_support_points(const unsigned int deg) { this->generalized_support_points.resize(this->n_dofs_per_cell()); - this->generalized_face_support_points.resize(this->n_dofs_per_face()); + this->generalized_face_support_points[0].resize(this->n_dofs_per_face()); // Number of the point being entered unsigned int current = 0; @@ -164,7 +164,7 @@ FE_RaviartThomasNodal::initialize_support_points(const unsigned int deg) 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[k] = face_points.point(k); + this->generalized_face_support_points[0][k] = face_points.point(k); Quadrature faces = QProjector::project_to_all_faces(this->reference_cell_type(), face_points); @@ -614,7 +614,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); + source_fe.generalized_face_support_points[0]); // Rule of thumb for FP accuracy, // that can be expected for a @@ -721,7 +721,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); + source_fe.generalized_face_support_points[0]); // Rule of thumb for FP accuracy, // that can be expected for a diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc index b7c66bdf63..9f0b2025db 100644 --- a/source/fe/fe_rt_bubbles.cc +++ b/source/fe/fe_rt_bubbles.cc @@ -132,7 +132,7 @@ void FE_RT_Bubbles::initialize_support_points(const unsigned int deg) { this->generalized_support_points.resize(this->n_dofs_per_cell()); - this->generalized_face_support_points.resize(this->n_dofs_per_face()); + this->generalized_face_support_points[0].resize(this->n_dofs_per_face()); // Index of the point being entered unsigned int current = 0; @@ -145,7 +145,7 @@ FE_RT_Bubbles::initialize_support_points(const unsigned int deg) 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[k] = face_points.point(k); + this->generalized_face_support_points[0][k] = face_points.point(k); Quadrature faces = QProjector::project_to_all_faces(this->reference_cell_type(), face_points); diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index ef668c0171..f0392cccaa 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -1422,7 +1422,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[n]; + n_index = this->face_system_to_base_table[0][n]; // likewise for the m index. this is more complicated due to the // strange ordering we have for the dofs on the refined faces. @@ -1460,7 +1460,8 @@ FESystem::build_interface_constraints() // then use face_system_to_base_table const unsigned int tmp1 = 2 * this->n_dofs_per_vertex() + index_in_line; - m_index.first = this->face_system_to_base_table[tmp1].first; + m_index.first = + this->face_system_to_base_table[0][tmp1].first; // what we are still missing is the index of m within the // base elements interface_constraints table @@ -1471,12 +1472,12 @@ 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[tmp1].second >= + Assert(this->face_system_to_base_table[0][tmp1].second >= 2 * base_element(m_index.first.first) .n_dofs_per_vertex(), ExcInternalError()); const unsigned int tmp2 = - this->face_system_to_base_table[tmp1].second - + this->face_system_to_base_table[0][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(), @@ -1515,14 +1516,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[tmp1].first; + m_index.first = + this->face_system_to_base_table[0][tmp1].first; - Assert(this->face_system_to_base_table[tmp1].second >= + Assert(this->face_system_to_base_table[0][tmp1].second >= 4 * base_element(m_index.first.first) .n_dofs_per_vertex(), ExcInternalError()); const unsigned int tmp2 = - this->face_system_to_base_table[tmp1].second - + this->face_system_to_base_table[0][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(), @@ -1553,18 +1555,19 @@ FESystem::build_interface_constraints() 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.size(), + Assert(tmp1 < this->face_system_to_base_table[0].size(), ExcInternalError()); - m_index.first = this->face_system_to_base_table[tmp1].first; + m_index.first = + this->face_system_to_base_table[0][tmp1].first; - Assert(this->face_system_to_base_table[tmp1].second >= + 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()); const unsigned int tmp2 = - this->face_system_to_base_table[tmp1].second - + 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(); @@ -1644,7 +1647,9 @@ FESystem::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.resize(this->n_dofs_per_face()); + + 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, @@ -1652,8 +1657,8 @@ FESystem::initialize( *this); FETools::Compositing::build_face_tables( - this->face_system_to_base_table, - this->face_system_to_component_table, + this->face_system_to_base_table[0], + this->face_system_to_component_table[0], *this); } @@ -1712,28 +1717,28 @@ FESystem::initialize( if (!base_element(base_el).has_support_points() && (base_element(base_el).n_dofs_per_face() > 0)) { - this->unit_face_support_points.resize(0); + this->unit_face_support_points[0].resize(0); return; } // generate unit face support points from unit support points of sub // elements - this->unit_face_support_points.resize(this->n_dofs_per_face()); + 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[i].first.first; + this->face_system_to_base_table[0][i].first.first; const unsigned int index_in_base = - this->face_system_to_base_table[i].second; + this->face_system_to_base_table[0][i].second; Assert(index_in_base < - base_element(base_i).unit_face_support_points.size(), + base_element(base_i).unit_face_support_points[0].size(), ExcInternalError()); - this->unit_face_support_points[i] = - base_element(base_i).unit_face_support_points[index_in_base]; + this->unit_face_support_points[0][i] = + base_element(base_i).unit_face_support_points[0][index_in_base]; } }); @@ -1810,7 +1815,7 @@ FESystem::initialize( 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 + Assert(this->adjust_quad_dof_index_for_face_orientation_table[0] .n_elements() == 8 * this->n_dofs_per_quad(), ExcInternalError()); @@ -1821,12 +1826,12 @@ FESystem::initialize( { const Table<2, int> &temp = this->base_element(b) - .adjust_quad_dof_index_for_face_orientation_table; + .adjust_quad_dof_index_for_face_orientation_table[0]; 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( + this->adjust_quad_dof_index_for_face_orientation_table[0]( index + i, j) = temp(i, j); index += temp.size(0); } @@ -2343,13 +2348,14 @@ Point FESystem::unit_face_support_point(const unsigned int index) const { AssertIndexRange(index, this->n_dofs_per_face()); - Assert((this->unit_face_support_points.size() == this->n_dofs_per_face()) || - (this->unit_face_support_points.size() == 0), + Assert((this->unit_face_support_points[0].size() == + this->n_dofs_per_face()) || + (this->unit_face_support_points[0].size() == 0), (typename FiniteElement::ExcFEHasNoSupportPoints())); // let's see whether we have the information pre-computed - if (this->unit_face_support_points.size() != 0) - return this->unit_face_support_points[index]; + if (this->unit_face_support_points[0].size() != 0) + return this->unit_face_support_points[0][index]; else // no. ask the base element whether it would like to provide this // information diff --git a/source/fe/fe_trace.cc b/source/fe/fe_trace.cc index 577c3ae778..0f990292a5 100644 --- a/source/fe/fe_trace.cc +++ b/source/fe/fe_trace.cc @@ -53,7 +53,7 @@ FE_TraceQ::FE_TraceQ(const unsigned int degree) FETools::hierarchic_to_lexicographic_numbering(degree)); // Initialize face support points - this->unit_face_support_points = fe_q.get_unit_face_support_points(); + this->unit_face_support_points[0] = fe_q.get_unit_face_support_points(); // 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 diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 75af427581..a3545d4e8d 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -120,8 +120,8 @@ namespace Simplex this->unit_support_points.emplace_back(0.0, 0.0); // TODO - this->unit_face_support_points.emplace_back(0.0); - this->unit_face_support_points.emplace_back(1.0); + this->unit_face_support_points[0].emplace_back(0.0); + this->unit_face_support_points[0].emplace_back(1.0); } else if (degree == 2) { @@ -133,9 +133,9 @@ namespace Simplex this->unit_support_points.emplace_back(0.5, 0.0); // TODO - this->unit_face_support_points.emplace_back(0.0); - this->unit_face_support_points.emplace_back(1.0); - this->unit_face_support_points.emplace_back(0.5); + this->unit_face_support_points[0].emplace_back(0.0); + this->unit_face_support_points[0].emplace_back(1.0); + this->unit_face_support_points[0].emplace_back(0.5); } else { @@ -152,9 +152,9 @@ namespace Simplex this->unit_support_points.emplace_back(0.0, 0.0, 1.0); // TODO - this->unit_face_support_points.emplace_back(1.0, 0.0); - this->unit_face_support_points.emplace_back(0.0, 1.0); - this->unit_face_support_points.emplace_back(0.0, 0.0); + this->unit_face_support_points[0].emplace_back(1.0, 0.0); + this->unit_face_support_points[0].emplace_back(0.0, 1.0); + this->unit_face_support_points[0].emplace_back(0.0, 0.0); } else if (degree == 2) { @@ -170,12 +170,12 @@ namespace Simplex this->unit_support_points.emplace_back(0.0, 0.5, 0.5); // TODO - this->unit_face_support_points.emplace_back(1.0, 0.0); - this->unit_face_support_points.emplace_back(0.0, 1.0); - this->unit_face_support_points.emplace_back(0.0, 0.0); - this->unit_face_support_points.emplace_back(0.5, 0.5); - this->unit_face_support_points.emplace_back(0.0, 0.5); - this->unit_face_support_points.emplace_back(0.5, 0.0); + this->unit_face_support_points[0].emplace_back(1.0, 0.0); + this->unit_face_support_points[0].emplace_back(0.0, 1.0); + this->unit_face_support_points[0].emplace_back(0.0, 0.0); + this->unit_face_support_points[0].emplace_back(0.5, 0.5); + this->unit_face_support_points[0].emplace_back(0.0, 0.5); + this->unit_face_support_points[0].emplace_back(0.5, 0.0); } else { -- 2.39.5