From 0dbbb18bd9903a514ba3d21c1c5949c6bd0986cb Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 1 Aug 2020 14:43:56 +0200 Subject: [PATCH] Pass face_no to FE functions --- include/deal.II/dofs/dof_accessor.templates.h | 16 +- .../matrix_free/cuda_hanging_nodes_internal.h | 24 +- .../matrix_free/matrix_free.templates.h | 3 +- .../deal.II/multigrid/mg_constrained_dofs.h | 2 +- .../numerics/matrix_creator.templates.h | 9 +- .../vector_tools_boundary.templates.h | 100 +++++---- .../vector_tools_constraints.templates.h | 58 +++-- .../numerics/vector_tools_project.templates.h | 4 +- source/dofs/dof_handler.cc | 52 +++-- source/dofs/dof_handler_policy.cc | 32 ++- source/dofs/dof_tools.cc | 31 ++- source/dofs/dof_tools_constraints.cc | 206 ++++++++++++------ source/dofs/dof_tools_sparsity.cc | 6 +- source/grid/grid_tools_dof_handlers.cc | 6 +- source/multigrid/mg_tools.cc | 89 +++++--- 15 files changed, 392 insertions(+), 246 deletions(-) diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index c7c9cc567f..e84ffdbda9 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -351,7 +351,7 @@ namespace internal dof_handler.object_dof_ptr[obj_level][d][obj_index]; const unsigned int ptr_1 = ptr_0 + - dof_handler.get_fe(fe_index).template n_dofs_per_object(); + dof_handler.get_fe(fe_index).template n_dofs_per_object(0); return {ptr_0, ptr_1}; } @@ -819,11 +819,11 @@ namespace internal { const auto &fe = accessor.get_fe(fe_index_); - const unsigned int // - dofs_per_vertex = fe.n_dofs_per_vertex(), // - dofs_per_line = fe.n_dofs_per_line(), // - dofs_per_quad = fe.n_dofs_per_quad(), // - dofs_per_hex = fe.n_dofs_per_hex(); // + const unsigned int // + dofs_per_vertex = fe.n_dofs_per_vertex(), // + dofs_per_line = fe.n_dofs_per_line(), // + dofs_per_quad = fe.n_dofs_per_quad(0 /*dummy*/), // + dofs_per_hex = fe.n_dofs_per_hex(); // const unsigned int inner_dofs = structdim == 1 ? dofs_per_line : @@ -1211,7 +1211,7 @@ namespace internal (void)fe_index; for (unsigned int d = 0; - d < fe.template n_dofs_per_object(); + d < fe.template n_dofs_per_object(0); ++d, ++index) index_value[index] = accessor.mg_dof_index(level, mapping(d)); } @@ -1291,7 +1291,7 @@ namespace internal (void)fe_index; for (unsigned int d = 0; - d < fe.template n_dofs_per_object(); + d < fe.template n_dofs_per_object(0); ++d, ++index) accessor.set_mg_dof_index(level, mapping(d), index_value[index]); } diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h index 7e993a6d11..ab14ffdbb2 100644 --- a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h @@ -140,25 +140,31 @@ namespace CUDAWrappers void setup_constraint_weigths(unsigned int fe_degree) { + const unsigned int face_no = + 0; // we assume that all faces have the same number of dofs + FE_Q<2> fe_q(fe_degree); - FullMatrix interpolation_matrix(fe_q.n_dofs_per_face(), - fe_q.n_dofs_per_face()); - fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix); + FullMatrix interpolation_matrix(fe_q.n_dofs_per_face(face_no), + fe_q.n_dofs_per_face(face_no)); + fe_q.get_subface_interpolation_matrix(fe_q, + 0, + interpolation_matrix, + face_no); std::vector mapping = FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree); - FullMatrix mapped_matrix(fe_q.n_dofs_per_face(), - fe_q.n_dofs_per_face()); - for (unsigned int i = 0; i < fe_q.n_dofs_per_face(); ++i) - for (unsigned int j = 0; j < fe_q.n_dofs_per_face(); ++j) + FullMatrix mapped_matrix(fe_q.n_dofs_per_face(face_no), + fe_q.n_dofs_per_face(face_no)); + for (unsigned int i = 0; i < fe_q.n_dofs_per_face(face_no); ++i) + for (unsigned int j = 0; j < fe_q.n_dofs_per_face(face_no); ++j) mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]); cudaError_t error_code = cudaMemcpyToSymbol(internal::constraint_weights, &mapped_matrix[0][0], - sizeof(double) * fe_q.n_dofs_per_face() * - fe_q.n_dofs_per_face()); + sizeof(double) * fe_q.n_dofs_per_face(face_no) * + fe_q.n_dofs_per_face(face_no)); AssertCuda(error_code); } } // namespace internal diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 19e6710dac..8519cb123c 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -760,7 +760,8 @@ namespace internal const FiniteElement &fe = fes[fe_index]; // cache number of finite elements and dofs_per_cell dof_info[no].dofs_per_cell.push_back(fe.n_dofs_per_cell()); - dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face()); + dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face( + 0)); // we assume that all faces have the same number of dofs dof_info[no].dimension = dim; dof_info[no].n_base_elements = fe.n_base_elements(); dof_info[no].n_components.resize(dof_info[no].n_base_elements); diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 7baf0b58b6..53e9d88dc7 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -273,7 +273,7 @@ MGConstrainedDoFs::initialize(const DoFHandler &dof) } const unsigned int dofs_per_face = - cell->face(f)->get_fe(0).n_dofs_per_face(); + dof.get_fe(0).n_dofs_per_face(f); std::vector dofs_1(dofs_per_face); std::vector dofs_2(dofs_per_face); diff --git a/include/deal.II/numerics/matrix_creator.templates.h b/include/deal.II/numerics/matrix_creator.templates.h index 3b37e5b399..57f067889e 100644 --- a/include/deal.II/numerics/matrix_creator.templates.h +++ b/include/deal.II/numerics/matrix_creator.templates.h @@ -948,8 +948,6 @@ namespace MatrixCreator const bool fe_is_system = (n_components != 1); const bool fe_is_primitive = fe.is_primitive(); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); - copy_data.cell = cell; copy_data.dofs_per_cell = fe.n_dofs_per_cell(); @@ -973,7 +971,7 @@ namespace MatrixCreator copy_data.dofs.resize(copy_data.dofs_per_cell); cell->get_dof_indices(copy_data.dofs); - std::vector dofs_on_face_vector(dofs_per_face); + std::vector dofs_on_face_vector; // Because CopyData objects are reused and emplace_back is // used, dof_is_on_face, cell_matrix, and cell_vector must be @@ -1147,6 +1145,7 @@ namespace MatrixCreator } + dofs_on_face_vector.resize(fe.n_dofs_per_face(face)); cell->face(face)->get_dof_indices(dofs_on_face_vector); // for each dof on the cell, have a flag whether it is on // the face @@ -1379,7 +1378,6 @@ namespace MatrixCreator const FiniteElement &fe = cell->get_fe(); const bool fe_is_system = (n_components != 1); const bool fe_is_primitive = fe.is_primitive(); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); copy_data.cell = cell; copy_data.dofs_per_cell = fe.n_dofs_per_cell(); @@ -1407,7 +1405,7 @@ namespace MatrixCreator std::vector rhs_values_scalar; std::vector> rhs_values_system; - std::vector dofs_on_face_vector(dofs_per_face); + std::vector dofs_on_face_vector; copy_data.dofs.resize(copy_data.dofs_per_cell); cell->get_dof_indices(copy_data.dofs); @@ -1599,6 +1597,7 @@ namespace MatrixCreator } } + dofs_on_face_vector.resize(fe.n_dofs_per_face(face)); cell->face(face)->get_dof_indices(dofs_on_face_vector, cell->active_fe_index()); // for each dof on the cell, have a diff --git a/include/deal.II/numerics/vector_tools_boundary.templates.h b/include/deal.II/numerics/vector_tools_boundary.templates.h index 0e8c79faf0..7801f0d42e 100644 --- a/include/deal.II/numerics/vector_tools_boundary.templates.h +++ b/include/deal.II/numerics/vector_tools_boundary.templates.h @@ -131,12 +131,14 @@ namespace VectorTools function_values); for (unsigned int i = 0; i < fe.n_dofs_per_vertex(); ++i) - if (component_mask[fe.face_system_to_component_index(i) + if (component_mask[fe.face_system_to_component_index( + i, direction) .first]) boundary_values[cell->vertex_dof_index( direction, i, cell->active_fe_index())] = function_values( - fe.face_system_to_component_index(i).first); + fe.face_system_to_component_index(i, direction) + .first); } } else // dim > 1 @@ -157,6 +159,9 @@ namespace VectorTools dof_values_system.reserve( dof.get_fe_collection().max_dofs_per_face()); + // TODO: get support for each face -> PR #10764 + const unsigned int face_no = 0; + // before we start with the loop over all cells create an hp::FEValues // object that holds the interpolation points of all finite elements // that may ever be in use @@ -173,9 +178,9 @@ namespace VectorTools // // to do this, we check whether the FE has support points on the // face at all: - if (fe.has_face_support_points()) - q_collection.push_back( - Quadrature(fe.get_unit_face_support_points())); + if (fe.has_face_support_points(face_no)) + q_collection.push_back(Quadrature( + fe.get_unit_face_support_points(face_no))); else { // if not, then we should try a more clever way. the idea is @@ -194,13 +199,15 @@ namespace VectorTools // components we are interested in are primitive (by the above // check), we can safely put such a check in front std::vector> unit_support_points( - fe.n_dofs_per_face()); + fe.n_dofs_per_face(face_no)); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) - if (fe.is_primitive(fe.face_to_cell_index(i, 0))) - if (component_mask[fe.face_system_to_component_index(i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i) + if (fe.is_primitive(fe.face_to_cell_index(i, face_no))) + if (component_mask[fe.face_system_to_component_index( + i, face_no) .first] == true) - unit_support_points[i] = fe.unit_face_support_point(i); + unit_support_points[i] = + fe.unit_face_support_point(i, face_no); q_collection.push_back( Quadrature(unit_support_points)); @@ -260,7 +267,7 @@ namespace VectorTools // element in use here has DoFs on the face at all if ((function_map.find(boundary_component) != function_map.end()) && - (cell->get_fe().n_dofs_per_face() > 0)) + (cell->get_fe().n_dofs_per_face(face_no) > 0)) { // face is of the right component x_fe_values.reinit(cell, face_no); @@ -269,7 +276,7 @@ namespace VectorTools // get indices, physical location and boundary values of // dofs on this face - face_dofs.resize(fe.n_dofs_per_face()); + face_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_dof_indices(face_dofs, cell->active_fe_index()); const std::vector> &dof_locations = fe_values.get_quadrature_points(); @@ -278,12 +285,14 @@ namespace VectorTools { // resize array. avoid construction of a memory // allocating temporary if possible - if (dof_values_system.size() < fe.n_dofs_per_face()) - dof_values_system.resize(fe.n_dofs_per_face(), - Vector( - fe.n_components())); + if (dof_values_system.size() < + fe.n_dofs_per_face(face_no)) + dof_values_system.resize( + fe.n_dofs_per_face(face_no), + Vector(fe.n_components())); else - dof_values_system.resize(fe.n_dofs_per_face()); + dof_values_system.resize( + fe.n_dofs_per_face(face_no)); function_map.find(boundary_component) ->second->vector_value_list(dof_locations, @@ -298,7 +307,8 @@ namespace VectorTools unsigned int component; if (fe.is_primitive()) component = - fe.face_system_to_component_index(i).first; + fe.face_system_to_component_index(i, face_no) + .first; else { // non-primitive case. make sure that this @@ -360,7 +370,7 @@ namespace VectorTools // fe has only one component, so save some computations { // get only the one component that this function has - dof_values_scalar.resize(fe.n_dofs_per_face()); + dof_values_scalar.resize(fe.n_dofs_per_face(face_no)); function_map.find(boundary_component) ->second->value_list(dof_locations, dof_values_scalar, @@ -1026,7 +1036,7 @@ namespace VectorTools // Compute the degrees of // freedom. - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index(fe.face_to_cell_index(i, face)) .first == base_indices) && @@ -1201,7 +1211,7 @@ namespace VectorTools // Compute the degrees // of freedom. - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index( @@ -1271,7 +1281,7 @@ namespace VectorTools for (unsigned int d = 0; d < dim; ++d) tmp[d] = values[q_point](first_vector_component + d); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index( @@ -1337,7 +1347,7 @@ namespace VectorTools unsigned int index = 0; - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index( @@ -1391,7 +1401,7 @@ namespace VectorTools { unsigned int index = 0; - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index(fe.face_to_cell_index(i, face)) .first == base_indices) && @@ -1431,7 +1441,7 @@ namespace VectorTools for (unsigned int d = 0; d < dim; ++d) tmp[d] = values[q_point](first_vector_component + d); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index( @@ -1487,7 +1497,7 @@ namespace VectorTools unsigned int index = 0; - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index( @@ -1527,7 +1537,7 @@ namespace VectorTools unsigned int index = 0; - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (((dynamic_cast *>(&fe) != nullptr) && (fe.system_to_base_index(fe.face_to_cell_index(i, face)) .first == base_indices) && @@ -1583,7 +1593,12 @@ namespace VectorTools // corresponding degrees of freedom. const unsigned int superdegree = dof_handler.get_fe().degree; const QGauss reference_face_quadrature(2 * superdegree); - const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face(); + + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(dof_handler.get_fe().n_unique_faces(), 1); + const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face(0); + const hp::FECollection &fe_collection(dof_handler.get_fe_collection()); const hp::MappingCollection mapping_collection(mapping); hp::QCollection face_quadrature_collection; @@ -1860,7 +1875,7 @@ namespace VectorTools } const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face); dofs_processed.resize(dofs_per_face); dof_values.resize(dofs_per_face); @@ -1953,7 +1968,7 @@ namespace VectorTools const unsigned int superdegree = cell->get_fe().degree; const unsigned int degree = superdegree - 1; const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face); dofs_processed.resize(dofs_per_face); dof_values.resize(dofs_per_face); @@ -2428,7 +2443,8 @@ namespace VectorTools 1); unsigned int associated_edge_dof_index = 0; - for (unsigned int face_idx = 0; face_idx < fe.n_dofs_per_face(); + for (unsigned int face_idx = 0; + face_idx < fe.n_dofs_per_face(face); ++face_idx) { const unsigned int cell_idx = @@ -2663,7 +2679,7 @@ namespace VectorTools // Loop over these quad-interior dofs. unsigned int associated_face_dof_index = 0; for (unsigned int quad_dof_idx = 0; - quad_dof_idx < fe.n_dofs_per_quad(); + quad_dof_idx < fe.n_dofs_per_quad(face); ++quad_dof_idx) { const unsigned int face_idx = @@ -2925,7 +2941,7 @@ namespace VectorTools } const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face); dofs_processed.resize(dofs_per_face); dof_values.resize(dofs_per_face); @@ -3051,7 +3067,7 @@ namespace VectorTools cell->get_fe().degree; const unsigned int degree = superdegree - 1; const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face); dofs_processed.resize(dofs_per_face); dof_values.resize(dofs_per_face); @@ -3206,7 +3222,7 @@ namespace VectorTools 0}; std::vector> values(fe_values.n_quadrature_points, Vector(2)); - Vector dof_values(fe.n_dofs_per_face()); + Vector dof_values(fe.n_dofs_per_face(face)); // Get the values of the boundary function at the quadrature points. { @@ -3231,7 +3247,7 @@ namespace VectorTools jacobians[q_point][1][face_coordinate_direction[face]] * jacobians[q_point][1][face_coordinate_direction[face]]); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) dof_values(i) += tmp * (normals[q_point] * fe_values[vec].value( @@ -3244,14 +3260,14 @@ namespace VectorTools } std::vector face_dof_indices( - fe.n_dofs_per_face()); + fe.n_dofs_per_face(face)); cell->face(face)->get_dof_indices(face_dof_indices, cell->active_fe_index()); // Copy the computed values in the AffineConstraints only, if the degree // of freedom is not already constrained. - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (!(constraints.is_constrained(face_dof_indices[i])) && fe.get_nonzero_components(fe.face_to_cell_index( i, @@ -3307,7 +3323,7 @@ namespace VectorTools {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}}; std::vector> values(fe_values.n_quadrature_points, Vector(3)); - Vector dof_values_local(fe.n_dofs_per_face()); + Vector dof_values_local(fe.n_dofs_per_face(face)); { const std::vector> &quadrature_points = @@ -3340,7 +3356,7 @@ namespace VectorTools jacobians[q_point][2][face_coordinate_directions[face][1]] * jacobians[q_point][2][face_coordinate_directions[face][1]])); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) dof_values_local(i) += tmp * (normals[q_point] * fe_values[vec].value( @@ -3353,12 +3369,12 @@ namespace VectorTools } std::vector face_dof_indices( - fe.n_dofs_per_face()); + fe.n_dofs_per_face(face)); cell->face(face)->get_dof_indices(face_dof_indices, cell->active_fe_index()); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (projected_dofs[face_dof_indices[i]] < fe.degree && fe.get_nonzero_components(fe.face_to_cell_index( i, diff --git a/include/deal.II/numerics/vector_tools_constraints.templates.h b/include/deal.II/numerics/vector_tools_constraints.templates.h index 410f034c41..24477aac63 100644 --- a/include/deal.II/numerics/vector_tools_constraints.templates.h +++ b/include/deal.II/numerics/vector_tools_constraints.templates.h @@ -490,6 +490,11 @@ namespace VectorTools for (unsigned int i = 0; i < fe_collection.size(); ++i) mapping_collection.push_back(mapping); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(dof_handler.get_fe().n_unique_faces(), 1); + const unsigned int face_no = 0; + // now also create a quadrature collection for the faces of a cell. fill // it with a quadrature formula with the support points on faces for each // FE @@ -497,9 +502,10 @@ namespace VectorTools for (unsigned int i = 0; i < fe_collection.size(); ++i) { const std::vector> &unit_support_points = - fe_collection[i].get_unit_face_support_points(); + fe_collection[i].get_unit_face_support_points(face_no); - Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(), + Assert(unit_support_points.size() == + fe_collection[i].n_dofs_per_face(face_no), ExcInternalError()); face_quadrature_collection.push_back( @@ -549,7 +555,7 @@ namespace VectorTools cell->face(face_no); // get the indices of the dofs on this cell... - face_dofs.resize(fe.n_dofs_per_face()); + face_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_dof_indices(face_dofs, cell->active_fe_index()); x_fe_face_values.reinit(cell, face_no); @@ -559,7 +565,7 @@ namespace VectorTools // then identify which of them correspond to the selected set of // vector components for (unsigned int i = 0; i < face_dofs.size(); ++i) - if (fe.face_system_to_component_index(i).first == + if (fe.face_system_to_component_index(i, face_no).first == first_vector_component) { // find corresponding other components of vector @@ -572,18 +578,19 @@ namespace VectorTools "Error: the finite element does not have enough components " "to define a normal direction.")); - for (unsigned int k = 0; k < fe.n_dofs_per_face(); ++k) + for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no); + ++k) if ((k != i) && (face_quadrature_collection[cell->active_fe_index()] .point(k) == face_quadrature_collection[cell->active_fe_index()] .point(i)) && - (fe.face_system_to_component_index(k).first >= - first_vector_component) && - (fe.face_system_to_component_index(k).first < + (fe.face_system_to_component_index(k, face_no) + .first >= first_vector_component) && + (fe.face_system_to_component_index(k, face_no).first < first_vector_component + dim)) vector_dofs.dof_indices - [fe.face_system_to_component_index(k).first - + [fe.face_system_to_component_index(k, face_no).first - first_vector_component] = face_dofs[k]; for (unsigned int d = 0; d < dim; ++d) @@ -1076,6 +1083,11 @@ namespace VectorTools for (unsigned int i = 0; i < fe_collection.size(); ++i) mapping_collection.push_back(mapping); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(dof_handler.get_fe().n_unique_faces(), 1); + const unsigned int face_no = 0; + // now also create a quadrature collection for the faces of a cell. fill // it with a quadrature formula with the support points on faces for each // FE @@ -1083,9 +1095,10 @@ namespace VectorTools for (unsigned int i = 0; i < fe_collection.size(); ++i) { const std::vector> &unit_support_points = - fe_collection[i].get_unit_face_support_points(); + fe_collection[i].get_unit_face_support_points(face_no); - Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(), + Assert(unit_support_points.size() == + fe_collection[i].n_dofs_per_face(face_no), ExcInternalError()); face_quadrature_collection.push_back( @@ -1123,7 +1136,7 @@ namespace VectorTools cell->face(face_no); // get the indices of the dofs on this cell... - face_dofs.resize(fe.n_dofs_per_face()); + face_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_dof_indices(face_dofs, cell->active_fe_index()); x_fe_face_values.reinit(cell, face_no); @@ -1133,22 +1146,23 @@ namespace VectorTools std::map dof_to_b_value; unsigned int n_scalar_indices = 0; - cell_vector_dofs.resize(fe.n_dofs_per_face()); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + cell_vector_dofs.resize(fe.n_dofs_per_face(face_no)); + for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i) { - if (fe.face_system_to_component_index(i).first >= + if (fe.face_system_to_component_index(i, face_no).first >= first_vector_component && - fe.face_system_to_component_index(i).first < + fe.face_system_to_component_index(i, face_no).first < first_vector_component + dim) { const unsigned int component = - fe.face_system_to_component_index(i).first - + fe.face_system_to_component_index(i, face_no).first - first_vector_component; - n_scalar_indices = - std::max(n_scalar_indices, - fe.face_system_to_component_index(i).second + - 1); - cell_vector_dofs[fe.face_system_to_component_index(i) + n_scalar_indices = std::max( + n_scalar_indices, + fe.face_system_to_component_index(i, face_no).second + + 1); + cell_vector_dofs[fe.face_system_to_component_index( + i, face_no) .second][component] = face_dofs[i]; const Point point = fe_values.quadrature_point(i); diff --git a/include/deal.II/numerics/vector_tools_project.templates.h b/include/deal.II/numerics/vector_tools_project.templates.h index 6db682e4e1..56c522bcba 100644 --- a/include/deal.II/numerics/vector_tools_project.templates.h +++ b/include/deal.II/numerics/vector_tools_project.templates.h @@ -83,10 +83,10 @@ namespace VectorTools for (auto f : GeometryInfo::face_indices()) if (cell->at_boundary(f)) { - face_dof_indices.resize(cell->get_fe().n_dofs_per_face()); + face_dof_indices.resize(cell->get_fe().n_dofs_per_face(f)); cell->face(f)->get_dof_indices(face_dof_indices, cell->active_fe_index()); - for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face(); + for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face(f); ++i) // enter zero boundary values // for all boundary nodes diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index ea798e48a1..1335bd6f4c 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -309,7 +309,8 @@ namespace internal { dof_handler.object_dof_indices[i][2].resize( dof_handler.tria->n_raw_cells(i) * - dof_handler.get_fe().n_dofs_per_quad(), + dof_handler.get_fe().n_dofs_per_quad( + 0 /*note: in 2D there is only one quad*/), numbers::invalid_dof_index); dof_handler.object_dof_ptr[i][2].reserve( @@ -317,7 +318,8 @@ namespace internal for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1; j++) dof_handler.object_dof_ptr[i][2].push_back( - j * dof_handler.get_fe().n_dofs_per_quad()); + j * dof_handler.get_fe().n_dofs_per_quad( + 0 /*note: in 2D there is only one quad*/)); dof_handler.cell_dof_cache_indices[i].resize( dof_handler.tria->n_raw_cells(i) * @@ -421,7 +423,7 @@ namespace internal { const auto &face = cell->face(face_index); const auto n_dofs_per_quad = - dof_handler.get_fe().n_dofs_per_quad(/*face_index*/); + dof_handler.get_fe().n_dofs_per_quad(face_index); auto &n_dofs_per_quad_target = dof_handler.object_dof_ptr[0][2][face->index() + 1]; @@ -549,9 +551,10 @@ namespace internal std::make_unique< internal::DoFHandlerImplementation::DoFLevel<2>>()); dof_handler.mg_levels.back()->dof_object.dofs = - std::vector(tria.n_raw_quads(i) * - fe.n_dofs_per_quad(), - numbers::invalid_dof_index); + std::vector( + tria.n_raw_quads(i) * + fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/), + numbers::invalid_dof_index); } dof_handler.mg_faces = @@ -635,10 +638,13 @@ namespace internal std::vector(tria.n_raw_lines() * fe.n_dofs_per_line(), numbers::invalid_dof_index); - dof_handler.mg_faces->quads.dofs = - std::vector(tria.n_raw_quads() * - fe.n_dofs_per_quad(), - numbers::invalid_dof_index); + + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(fe.n_unique_faces(), 1); + dof_handler.mg_faces->quads.dofs = std::vector( + tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/), + numbers::invalid_dof_index); const unsigned int n_vertices = tria.n_vertices(); @@ -1331,17 +1337,18 @@ namespace internal fe_slots_needed = 1; n_face_slots += dof_handler.get_fe(cell->active_fe_index()) - .template n_dofs_per_object(); + .template n_dofs_per_object(face); } else { fe_slots_needed = 2; n_face_slots += dof_handler.get_fe(cell->active_fe_index()) - .template n_dofs_per_object() + + .template n_dofs_per_object(face) + dof_handler .get_fe(cell->neighbor(face)->active_fe_index()) - .template n_dofs_per_object(); + .template n_dofs_per_object( + cell->neighbor_face_no(face)); } // mark this face as visited @@ -1412,7 +1419,7 @@ namespace internal const unsigned int fe = cell->active_fe_index(); const unsigned int n_dofs = dof_handler.get_fe(fe) - .template n_dofs_per_object(); + .template n_dofs_per_object(face); const unsigned int offset = dof_handler .hp_object_fe_ptr[d][cell->face(face)->index()]; @@ -1426,20 +1433,25 @@ namespace internal } else { - unsigned int fe_1 = cell->active_fe_index(); + unsigned int fe_1 = cell->active_fe_index(); + unsigned int face_no_1 = face; unsigned int fe_2 = cell->neighbor(face)->active_fe_index(); + unsigned int face_no_2 = cell->neighbor_face_no(face); if (fe_2 < fe_1) - std::swap(fe_1, fe_2); + { + std::swap(fe_1, fe_2); + std::swap(face_no_1, face_no_2); + } const unsigned int n_dofs_1 = dof_handler.get_fe(fe_1) - .template n_dofs_per_object(); + .template n_dofs_per_object(face_no_1); const unsigned int n_dofs_2 = dof_handler.get_fe(fe_2) - .template n_dofs_per_object(); + .template n_dofs_per_object(face_no_2); const unsigned int offset = dof_handler @@ -2316,7 +2328,7 @@ DoFHandler::n_boundary_dofs() const if (face->at_boundary()) { const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(iface); dofs_on_face.resize(dofs_per_face); face->get_dof_indices(dofs_on_face, cell->active_fe_index()); @@ -2368,7 +2380,7 @@ DoFHandler::n_boundary_dofs( (boundary_ids.find(boundary_id) != boundary_ids.end())) { const unsigned int dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(iface); dofs_on_face.resize(dofs_per_face); face->get_dof_indices(dofs_on_face, cell->active_fe_index()); diff --git a/source/dofs/dof_handler_policy.cc b/source/dofs/dof_handler_policy.cc index e230c03396..c3fae71238 100644 --- a/source/dofs/dof_handler_policy.cc +++ b/source/dofs/dof_handler_policy.cc @@ -145,8 +145,12 @@ namespace internal ensure_existence_and_return_dof_identities( const FiniteElement &fe1, const FiniteElement &fe2, - std::unique_ptr & identities) + std::unique_ptr & identities, + const unsigned int face_no = numbers::invalid_unsigned_int) { + Assert(structdim == 2 || face_no == numbers::invalid_unsigned_int, + ExcInternalError()); + // see if we need to fill this entry, or whether it already // exists if (identities.get() == nullptr) @@ -170,7 +174,7 @@ namespace internal case 2: { identities = std::make_unique( - fe1.hp_quad_dof_identities(fe2)); + fe1.hp_quad_dof_identities(fe2, face_no)); break; } @@ -183,10 +187,10 @@ namespace internal for (unsigned int i = 0; i < identities->size(); ++i) { Assert((*identities)[i].first < - fe1.template n_dofs_per_object(), + fe1.template n_dofs_per_object(face_no), ExcInternalError()); Assert((*identities)[i].second < - fe2.template n_dofs_per_object(), + fe2.template n_dofs_per_object(face_no), ExcInternalError()); } } @@ -807,6 +811,11 @@ namespace internal fe_indices, /*codim=*/dim - 2); + const unsigned int most_dominating_fe_index_face_no = + cell->active_fe_index() == most_dominating_fe_index ? + q : + cell->neighbor_face_no(q); + // if we found the most dominating element, then use // this to eliminate some of the degrees of freedom // by identification. otherwise, the code that @@ -829,7 +838,8 @@ namespace internal quad_dof_identities [most_dominating_fe_index][other_fe_index] [cell->quad(q)->reference_cell_type() == - ReferenceCell::Type::Quad]); + ReferenceCell::Type::Quad], + most_dominating_fe_index_face_no); for (const auto &identity : identities) { @@ -1558,6 +1568,11 @@ namespace internal fe_indices, /*codim=*/dim - 2); + const unsigned int most_dominating_fe_index_face_no = + cell->active_fe_index() == most_dominating_fe_index ? + q : + cell->neighbor_face_no(q); + // if we found the most dominating element, then use // this to eliminate some of the degrees of freedom // by identification. otherwise, the code that @@ -1580,7 +1595,8 @@ namespace internal quad_dof_identities [most_dominating_fe_index][other_fe_index] [cell->quad(q)->reference_cell_type() == - ReferenceCell::Type::Quad]); + ReferenceCell::Type::Quad], + most_dominating_fe_index_face_no); for (const auto &identity : identities) { @@ -2312,7 +2328,7 @@ namespace internal for (unsigned int d = 0; d < - dof_handler.get_fe(fe_index).n_dofs_per_quad(); + dof_handler.get_fe(fe_index).n_dofs_per_quad(q); ++d) { const types::global_dof_index old_dof_index = @@ -2674,7 +2690,7 @@ namespace internal if (cell->quad(l)->user_flag_set()) { for (unsigned int d = 0; - d < dof_handler.get_fe().n_dofs_per_quad(); + d < dof_handler.get_fe().n_dofs_per_quad(l); ++d) { const dealii::types::global_dof_index idx = diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index 341edb4396..fce2cbdfa8 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -712,12 +712,12 @@ namespace DoFTools { const FiniteElement &fe = cell->get_fe(); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); + const unsigned int dofs_per_face = fe.n_dofs_per_face(face); face_dof_indices.resize(dofs_per_face); cell->face(face)->get_dof_indices(face_dof_indices, cell->active_fe_index()); - for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) if (!check_vector_component) selected_dofs.add_index(face_dof_indices[i]); else @@ -744,9 +744,9 @@ namespace DoFTools numbers::invalid_unsigned_int))); if (fe.is_primitive(cell_index)) { - if (component_mask - [fe.face_system_to_component_index(i).first] == - true) + if (component_mask[fe.face_system_to_component_index( + i, face) + .first] == true) selected_dofs.add_index(face_dof_indices[i]); } else // not primitive @@ -1041,7 +1041,7 @@ namespace DoFTools { // simply take all DoFs that live on this subface std::vector ldi( - fe.n_dofs_per_face()); + fe.n_dofs_per_face(f, child)); face->child(child)->get_dof_indices(ldi); selected_dofs.add_indices(ldi.begin(), ldi.end()); } @@ -2111,7 +2111,8 @@ namespace DoFTools for (const unsigned int f : cell->face_indices()) if (cell->at_boundary(f)) { - const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); + const unsigned int dofs_per_face = + cell->get_fe().n_dofs_per_face(f); dofs_on_face.resize(dofs_per_face); cell->face(f)->get_dof_indices(dofs_on_face, cell->active_fe_index()); @@ -2156,7 +2157,8 @@ namespace DoFTools if (boundary_ids.find(cell->face(f)->boundary_id()) != boundary_ids.end()) { - const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); + const unsigned int dofs_per_face = + cell->get_fe().n_dofs_per_face(f); dofs_on_face.resize(dofs_per_face); cell->face(f)->get_dof_indices(dofs_on_face, cell->active_fe_index()); @@ -2533,12 +2535,11 @@ namespace DoFTools // Exclude degrees of freedom on faces opposite to the vertex exclude.resize(fe.n_dofs_per_cell()); std::fill(exclude.begin(), exclude.end(), false); - const unsigned int dpf = fe.n_dofs_per_face(); for (const unsigned int face : cell->face_indices()) if (cell->at_boundary(face) || cell->neighbor(face)->level() != cell->level()) - for (unsigned int i = 0; i < dpf; ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) exclude[fe.face_to_cell_index(i, face)] = true; for (types::global_dof_index j = 0; j < indices.size(); ++j) if (!exclude[j]) @@ -2594,13 +2595,11 @@ namespace DoFTools { // Eliminate dofs on faces of the child which are on faces // of the parent - const unsigned int dpf = fe.n_dofs_per_face(); - for (unsigned int d = 0; d < dim; ++d) { const unsigned int face = GeometryInfo::vertex_to_face[child][d]; - for (unsigned int i = 0; i < dpf; ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) exclude[fe.face_to_cell_index(i, face)] = true; } @@ -2610,7 +2609,8 @@ namespace DoFTools for (const unsigned int face : GeometryInfo::face_indices()) if (cell->at_boundary(face)) - for (unsigned int i = 0; i < dpf; ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); + ++i) exclude[fe.face_to_cell_index(i, face)] = false; } @@ -2750,7 +2750,6 @@ namespace DoFTools // vertex exclude.resize(fe.n_dofs_per_cell()); std::fill(exclude.begin(), exclude.end(), false); - const unsigned int dpf = fe.n_dofs_per_face(); for (unsigned int d = 0; d < dim; ++d) { @@ -2758,7 +2757,7 @@ namespace DoFTools GeometryInfo::vertex_to_face[v][d]; const unsigned int face = GeometryInfo::opposite_face[a_face]; - for (unsigned int i = 0; i < dpf; ++i) + for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i) { // For each dof, get the block it is in and decide to // exclude it or not diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 3fd168895c..ae4f919a88 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -175,15 +175,23 @@ namespace DoFTools const FullMatrix & face_interpolation_matrix, std::vector & primary_dof_mask) { - Assert(fe1.n_dofs_per_face() >= fe2.n_dofs_per_face(), + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(fe1.n_unique_faces(), 1); + AssertDimension(fe2.n_unique_faces(), 1); + const unsigned int face_no = 0; + (void)face_no; + + Assert(fe1.n_dofs_per_face(face_no) >= fe2.n_dofs_per_face(face_no), ExcInternalError()); - AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face()); + AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face(face_no)); Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(), ExcInternalError()); Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(), ExcInternalError()); - Assert((dim < 3) || (fe2.n_dofs_per_quad() <= fe1.n_dofs_per_quad()), + Assert((dim < 3) || + (fe2.n_dofs_per_quad(face_no) <= fe1.n_dofs_per_quad(face_no)), ExcInternalError()); // the idea here is to designate as many DoFs in fe1 per object (vertex, @@ -268,9 +276,9 @@ namespace DoFTools // same algorithm as above unsigned int dofs_added = 0; unsigned int i = 0; - while (dofs_added < fe2.n_dofs_per_quad()) + while (dofs_added < fe2.n_dofs_per_quad(q)) { - Assert(i < fe1.n_dofs_per_quad(), ExcInternalError()); + Assert(i < fe1.n_dofs_per_quad(q), ExcInternalError()); primary_dof_list.push_back(index + i); if (check_primary_dof_list(face_interpolation_matrix, @@ -281,11 +289,11 @@ namespace DoFTools ++i; } - index += fe1.n_dofs_per_quad(); + index += fe1.n_dofs_per_quad(q); } - AssertDimension(index, fe1.n_dofs_per_face()); - AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face()); + AssertDimension(index, fe1.n_dofs_per_face(face_no)); + AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face(face_no)); // finally copy the list into the mask std::fill(primary_dof_mask.begin(), primary_dof_mask.end(), false); @@ -307,10 +315,16 @@ namespace DoFTools const FullMatrix & face_interpolation_matrix, std::unique_ptr> &primary_dof_mask) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(fe1.n_unique_faces(), 1); + AssertDimension(fe2.n_unique_faces(), 1); + const unsigned int face_no = 0; + if (primary_dof_mask == nullptr) { primary_dof_mask = - std::make_unique>(fe1.n_dofs_per_face()); + std::make_unique>(fe1.n_dofs_per_face(face_no)); select_primary_dofs_for_face_restriction(fe1, fe2, face_interpolation_matrix, @@ -332,12 +346,17 @@ namespace DoFTools const FiniteElement & fe2, std::unique_ptr> &matrix) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(fe1.n_unique_faces(), 1); + AssertDimension(fe2.n_unique_faces(), 1); + const unsigned int face_no = 0; + if (matrix == nullptr) { - matrix = - std::make_unique>(fe2.n_dofs_per_face(), - fe1.n_dofs_per_face()); - fe1.get_face_interpolation_matrix(fe2, *matrix); + matrix = std::make_unique>( + fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no)); + fe1.get_face_interpolation_matrix(fe2, *matrix, face_no); } } @@ -354,12 +373,20 @@ namespace DoFTools const unsigned int subface, std::unique_ptr> &matrix) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension(fe1.n_unique_faces(), 1); + AssertDimension(fe2.n_unique_faces(), 1); + const unsigned int face_no = 0; + if (matrix == nullptr) { - matrix = - std::make_unique>(fe2.n_dofs_per_face(), - fe1.n_dofs_per_face()); - fe1.get_subface_interpolation_matrix(fe2, subface, *matrix); + matrix = std::make_unique>( + fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no)); + fe1.get_subface_interpolation_matrix(fe2, + subface, + *matrix, + face_no); } } @@ -774,7 +801,7 @@ namespace DoFTools // first of all, make sure that we treat a case which is // possible, i.e. either no dofs on the face at all or no // anisotropic refinement - if (cell->get_fe().n_dofs_per_face() == 0) + if (cell->get_fe().n_dofs_per_face(face) == 0) continue; Assert(cell->face(face)->refinement_case() == @@ -836,10 +863,10 @@ namespace DoFTools const FiniteElement &fe = cell->get_fe(); const unsigned int fe_index = cell->active_fe_index(); - const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(); + const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face); const unsigned int n_dofs_on_children = (5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() + - 4 * fe.n_dofs_per_quad()); + 4 * fe.n_dofs_per_quad(face)); // TODO[TL]: think about this and the following in case of // anisotropic refinement @@ -872,7 +899,8 @@ namespace DoFTools for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof) dofs_on_mother[next_index++] = this_face->line(line)->dof_index(dof, fe_index); - for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); ++dof) + for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face); + ++dof) dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index); AssertDimension(next_index, dofs_on_mother.size()); @@ -937,7 +965,7 @@ namespace DoFTools if (cell->neighbor_child_on_subface(face, child) ->is_artificial()) continue; - for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); + for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face); ++dof) dofs_on_children.push_back( this_face->child(child)->dof_index(dof, fe_index)); @@ -1048,7 +1076,7 @@ namespace DoFTools // first of all, make sure that we treat a case which is // possible, i.e. either no dofs on the face at all or no // anisotropic refinement - if (cell->get_fe().n_dofs_per_face() == 0) + if (cell->get_fe().n_dofs_per_face(face) == 0) continue; Assert(cell->face(face)->refinement_case() == @@ -1120,7 +1148,8 @@ namespace DoFTools // // so we are going to constrain the DoFs on the face // children against the DoFs on the face itself - primary_dofs.resize(cell->get_fe().n_dofs_per_face()); + primary_dofs.resize( + cell->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( primary_dofs, cell->active_fe_index()); @@ -1166,7 +1195,7 @@ namespace DoFTools // the face DoFs from the cell DoFs. dependent_dofs.resize( subface->get_fe(subface_fe_index) - .n_dofs_per_face()); + .n_dofs_per_face(face, c)); subface->get_dof_indices(dependent_dofs, subface_fe_index); @@ -1275,8 +1304,8 @@ namespace DoFTools // first get the interpolation matrix from the mother to // the virtual dofs - Assert(dominating_fe.n_dofs_per_face() <= - cell->get_fe().n_dofs_per_face(), + Assert(dominating_fe.n_dofs_per_face(face) <= + cell->get_fe().n_dofs_per_face(face), ExcInternalError()); ensure_existence_of_face_matrix( @@ -1318,16 +1347,17 @@ namespace DoFTools // now compute the constraint matrix as the product // between the inverse matrix and the dependent part constraint_matrix.reinit( - cell->get_fe().n_dofs_per_face() - - dominating_fe.n_dofs_per_face(), - dominating_fe.n_dofs_per_face()); + cell->get_fe().n_dofs_per_face(face) - + dominating_fe.n_dofs_per_face(face), + dominating_fe.n_dofs_per_face(face)); restrict_mother_to_virtual_dependent.mmult( constraint_matrix, restrict_mother_to_virtual_primary_inv); // then figure out the global numbers of primary and // dependent dofs and apply constraints - scratch_dofs.resize(cell->get_fe().n_dofs_per_face()); + scratch_dofs.resize( + cell->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( scratch_dofs, cell->active_fe_index()); @@ -1335,7 +1365,7 @@ namespace DoFTools primary_dofs.clear(); dependent_dofs.clear(); for (unsigned int i = 0; - i < cell->get_fe().n_dofs_per_face(); + i < cell->get_fe().n_dofs_per_face(face); ++i) if ((*primary_dof_masks[dominating_fe_index] [cell @@ -1346,10 +1376,10 @@ namespace DoFTools dependent_dofs.push_back(scratch_dofs[i]); AssertDimension(primary_dofs.size(), - dominating_fe.n_dofs_per_face()); + dominating_fe.n_dofs_per_face(face)); AssertDimension(dependent_dofs.size(), - cell->get_fe().n_dofs_per_face() - - dominating_fe.n_dofs_per_face()); + cell->get_fe().n_dofs_per_face(face) - + dominating_fe.n_dofs_per_face(face)); filter_constraints(primary_dofs, dependent_dofs, @@ -1386,8 +1416,8 @@ namespace DoFTools // first get the interpolation matrix from the // subface to the virtual dofs - Assert(dominating_fe.n_dofs_per_face() <= - subface_fe.n_dofs_per_face(), + Assert(dominating_fe.n_dofs_per_face(face) <= + subface_fe.n_dofs_per_face(face), ExcInternalError()); ensure_existence_of_subface_matrix( dominating_fe, @@ -1402,14 +1432,15 @@ namespace DoFTools [dominating_fe_index][subface_fe_index][sf]); constraint_matrix.reinit( - subface_fe.n_dofs_per_face(), - dominating_fe.n_dofs_per_face()); + subface_fe.n_dofs_per_face(face), + dominating_fe.n_dofs_per_face(face)); restrict_subface_to_virtual.mmult( constraint_matrix, restrict_mother_to_virtual_primary_inv); - dependent_dofs.resize(subface_fe.n_dofs_per_face()); + dependent_dofs.resize( + subface_fe.n_dofs_per_face(face)); cell->face(face)->child(sf)->get_dof_indices( dependent_dofs, subface_fe_index); @@ -1471,7 +1502,7 @@ namespace DoFTools // Get DoFs on dominating and dominated side of the // face primary_dofs.resize( - cell->get_fe().n_dofs_per_face()); + cell->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( primary_dofs, cell->active_fe_index()); @@ -1482,7 +1513,7 @@ namespace DoFTools break; dependent_dofs.resize( - neighbor->get_fe().n_dofs_per_face()); + neighbor->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( dependent_dofs, neighbor->active_fe_index()); @@ -1589,8 +1620,8 @@ namespace DoFTools // first get the interpolation matrix from main FE // to the virtual dofs - Assert(dominating_fe.n_dofs_per_face() <= - cell->get_fe().n_dofs_per_face(), + Assert(dominating_fe.n_dofs_per_face(face) <= + cell->get_fe().n_dofs_per_face(face), ExcInternalError()); ensure_existence_of_face_matrix( @@ -1633,9 +1664,9 @@ namespace DoFTools // now compute the constraint matrix as the product // between the inverse matrix and the dependent part constraint_matrix.reinit( - cell->get_fe().n_dofs_per_face() - - dominating_fe.n_dofs_per_face(), - dominating_fe.n_dofs_per_face()); + cell->get_fe().n_dofs_per_face(face) - + dominating_fe.n_dofs_per_face(face), + dominating_fe.n_dofs_per_face(face)); restrict_mother_to_virtual_dependent.mmult( constraint_matrix, restrict_mother_to_virtual_primary_inv); @@ -1643,7 +1674,7 @@ namespace DoFTools // then figure out the global numbers of primary and // dependent dofs and apply constraints scratch_dofs.resize( - cell->get_fe().n_dofs_per_face()); + cell->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( scratch_dofs, cell->active_fe_index()); @@ -1651,7 +1682,7 @@ namespace DoFTools primary_dofs.clear(); dependent_dofs.clear(); for (unsigned int i = 0; - i < cell->get_fe().n_dofs_per_face(); + i < cell->get_fe().n_dofs_per_face(face); ++i) if ((*primary_dof_masks[dominating_fe_index] [cell->active_fe_index()]) @@ -1661,10 +1692,12 @@ namespace DoFTools dependent_dofs.push_back(scratch_dofs[i]); AssertDimension(primary_dofs.size(), - dominating_fe.n_dofs_per_face()); - AssertDimension(dependent_dofs.size(), - cell->get_fe().n_dofs_per_face() - - dominating_fe.n_dofs_per_face()); + dominating_fe.n_dofs_per_face( + face)); + AssertDimension( + dependent_dofs.size(), + cell->get_fe().n_dofs_per_face(face) - + dominating_fe.n_dofs_per_face(face)); filter_constraints(primary_dofs, dependent_dofs, @@ -1674,8 +1707,8 @@ namespace DoFTools // now do the same for another FE this is pretty // much the same we do above to resolve h-refinement // constraints - Assert(dominating_fe.n_dofs_per_face() <= - neighbor->get_fe().n_dofs_per_face(), + Assert(dominating_fe.n_dofs_per_face(face) <= + neighbor->get_fe().n_dofs_per_face(face), ExcInternalError()); ensure_existence_of_face_matrix( @@ -1692,15 +1725,15 @@ namespace DoFTools [neighbor->active_fe_index()]); constraint_matrix.reinit( - neighbor->get_fe().n_dofs_per_face(), - dominating_fe.n_dofs_per_face()); + neighbor->get_fe().n_dofs_per_face(face), + dominating_fe.n_dofs_per_face(face)); restrict_secondface_to_virtual.mmult( constraint_matrix, restrict_mother_to_virtual_primary_inv); dependent_dofs.resize( - neighbor->get_fe().n_dofs_per_face()); + neighbor->get_fe().n_dofs_per_face(face)); cell->face(face)->get_dof_indices( dependent_dofs, neighbor->active_fe_index()); @@ -1796,6 +1829,14 @@ namespace DoFTools Assert(face_1->n_active_fe_indices() == 1, ExcInternalError()); + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension( + face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1); + AssertDimension( + face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1); + const unsigned int face_no = 0; + // If face_2 does have children, then we need to iterate over these // children and set periodic constraints in the inverse direction: @@ -1806,7 +1847,8 @@ namespace DoFTools ExcNotImplemented()); const unsigned int dofs_per_face = - face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face(); + face_1->get_fe(face_1->nth_active_fe_index(0)) + .n_dofs_per_face(face_no); FullMatrix child_transformation(dofs_per_face, dofs_per_face); FullMatrix subface_interpolation(dofs_per_face, dofs_per_face); @@ -1817,7 +1859,10 @@ namespace DoFTools // interpolated from face_1 to face_2 by multiplying from the left // with the one that interpolates from face_2 to its child const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0)); - fe.get_subface_interpolation_matrix(fe, c, subface_interpolation); + fe.get_subface_interpolation_matrix(fe, + c, + subface_interpolation, + face_no); subface_interpolation.mmult(child_transformation, transformation); set_periodicity_constraints(face_1, @@ -1852,7 +1897,7 @@ namespace DoFTools "zero or equal to the number of components in the finite " "element.")); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); + const unsigned int dofs_per_face = fe.n_dofs_per_face(face_no); std::vector dofs_1(dofs_per_face); std::vector dofs_2(dofs_per_face); @@ -1925,7 +1970,8 @@ namespace DoFTools // Obey the component mask if ((component_mask.n_selected_components(fe.n_components()) != fe.n_components()) && - !component_mask[fe.face_system_to_component_index(i).first]) + !component_mask[fe.face_system_to_component_index(i, face_no) + .first]) continue; // We have to be careful to treat so called "identity @@ -2120,9 +2166,14 @@ namespace DoFTools const FullMatrix & matrix, const std::vector & first_vector_components) { + // 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; + Assert(matrix.m() == matrix.n(), ExcInternalError()); - const unsigned int n_dofs_per_face = fe.n_dofs_per_face(); + const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no); if (matrix.m() == n_dofs_per_face) { @@ -2144,7 +2195,7 @@ namespace DoFTools Assert(matrix.m() == spacedim, ExcInternalError()) Quadrature - quadrature(fe.get_unit_face_support_points()); + quadrature(fe.get_unit_face_support_points(face_no)); // have an array that stores the location of each vector-dof tuple we want // to rotate. @@ -2158,7 +2209,7 @@ namespace DoFTools std::vector::const_iterator comp_it = std::find(first_vector_components.begin(), first_vector_components.end(), - fe.face_system_to_component_index(i).first); + fe.face_system_to_component_index(i, face_no).first); if (comp_it != first_vector_components.end()) { const unsigned int first_vector_component = *comp_it; @@ -2176,12 +2227,13 @@ namespace DoFTools for (unsigned int k = 0; k < n_dofs_per_face; ++k) if ((k != i) && (quadrature.point(k) == quadrature.point(i)) && - (fe.face_system_to_component_index(k).first >= + (fe.face_system_to_component_index(k, face_no).first >= first_vector_component) && - (fe.face_system_to_component_index(k).first < + (fe.face_system_to_component_index(k, face_no).first < first_vector_component + spacedim)) { - vector_dofs[fe.face_system_to_component_index(k).first - + vector_dofs[fe.face_system_to_component_index(k, face_no) + .first - first_vector_component] = k; n_found++; if (n_found == dim) @@ -2221,6 +2273,14 @@ namespace DoFTools const std::vector & first_vector_components, const number periodicity_factor) { + // TODO: the implementation makes the assumption that all faces have the + // same number of dofs + AssertDimension( + face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1); + AssertDimension( + face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1); + const unsigned int face_no = 0; + static const int dim = FaceIterator::AccessorType::dimension; static const int spacedim = FaceIterator::AccessorType::space_dimension; @@ -2256,7 +2316,8 @@ namespace DoFTools { Assert(face_1->n_active_fe_indices() == 1, ExcInternalError()); const unsigned int n_dofs_per_face = - face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face(); + face_1->get_fe(face_1->nth_active_fe_index(0)) + .n_dofs_per_face(face_no); Assert(matrix.m() == 0 || (first_vector_components.empty() && @@ -2273,7 +2334,8 @@ namespace DoFTools { Assert(face_2->n_active_fe_indices() == 1, ExcInternalError()); const unsigned int n_dofs_per_face = - face_2->get_fe(face_2->nth_active_fe_index(0)).n_dofs_per_face(); + face_2->get_fe(face_2->nth_active_fe_index(0)) + .n_dofs_per_face(face_no); Assert(matrix.m() == 0 || (first_vector_components.empty() && @@ -2372,7 +2434,7 @@ namespace DoFTools face_2->get_fe(face_2->nth_active_fe_index(0)) : face_1->get_fe(face_1->nth_active_fe_index(0)); - const unsigned int n_dofs_per_face = fe.n_dofs_per_face(); + const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no); // Sometimes we just have nothing to do (for all finite elements, or // systems which accidentally don't have any dofs on the boundary). @@ -3408,7 +3470,7 @@ namespace DoFTools (face->boundary_id() == boundary_id))) { // get indices and physical location on this face - face_dofs.resize(fe.n_dofs_per_face()); + face_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_dof_indices(face_dofs, cell->active_fe_index()); // enter those dofs into the list that match the component diff --git a/source/dofs/dof_tools_sparsity.cc b/source/dofs/dof_tools_sparsity.cc index 165a15ca66..c61801583d 100644 --- a/source/dofs/dof_tools_sparsity.cc +++ b/source/dofs/dof_tools_sparsity.cc @@ -407,7 +407,8 @@ namespace DoFTools for (const unsigned int f : cell->face_indices()) if (cell->at_boundary(f)) { - const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); + const unsigned int dofs_per_face = + cell->get_fe().n_dofs_per_face(f); dofs_on_this_face.resize(dofs_per_face); cell->face(f)->get_dof_indices(dofs_on_this_face, cell->active_fe_index()); @@ -504,7 +505,8 @@ namespace DoFTools if (boundary_ids.find(cell->face(f)->boundary_id()) != boundary_ids.end()) { - const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face(); + const unsigned int dofs_per_face = + cell->get_fe().n_dofs_per_face(f); dofs_on_this_face.resize(dofs_per_face); cell->face(f)->get_dof_indices(dofs_on_this_face, cell->active_fe_index()); diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 88598a00cf..a4d99be7b8 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -1830,7 +1830,7 @@ namespace GridTools ExcInternalError()); const unsigned int n_dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(f, c); local_face_dof_indices.resize(n_dofs_per_face); cell->face(f)->child(c)->get_dof_indices( @@ -1865,7 +1865,7 @@ namespace GridTools unsigned int subface = neighbor_face_no_subface_no.second; const unsigned int n_dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face_no); local_face_dof_indices.resize(n_dofs_per_face); cell->neighbor(f)->face(face_no)->get_dof_indices( @@ -1884,7 +1884,7 @@ namespace GridTools // original cell { const unsigned int n_dofs_per_face = - cell->get_fe().n_dofs_per_face(); + cell->get_fe().n_dofs_per_face(face_no, c); local_face_dof_indices.resize(n_dofs_per_face); Assert(cell->neighbor(f) diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc index a1f6556517..09d532ae6f 100644 --- a/source/multigrid/mg_tools.cc +++ b/source/multigrid/mg_tools.cc @@ -162,10 +162,14 @@ namespace MGTools // n_dofs_per_vertex() is // arbitrary, not the other way // round. - // TODO: This assumes that even in hp context, the dofs per face - // coincide! + // TODO: This assumes that the dofs per face on all faces coincide! + const unsigned int face_no = 0; + + Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + unsigned int increment = - fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face(); + fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face(face_no); while (i < fe.get_first_line_index()) row_lengths[cell_indices[i++]] += increment; // From now on, if an object is @@ -178,23 +182,25 @@ namespace MGTools // In all other cases we // subtract adjacent faces to be // added in the loop below. - increment = (dim > 1) ? - fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face() : - fe.n_dofs_per_cell() - GeometryInfo::faces_per_cell * - fe.n_dofs_per_face(); - while (i < fe.get_first_quad_index()) + increment = + (dim > 1) ? + fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face(face_no) : + fe.n_dofs_per_cell() - + GeometryInfo::faces_per_cell * fe.n_dofs_per_face(face_no); + while (i < fe.get_first_quad_index(face_no)) row_lengths[cell_indices[i++]] += increment; // Now quads in 2D and 3D - increment = (dim > 2) ? - fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face() : - fe.n_dofs_per_cell() - GeometryInfo::faces_per_cell * - fe.n_dofs_per_face(); + increment = + (dim > 2) ? + fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face(face_no) : + fe.n_dofs_per_cell() - + GeometryInfo::faces_per_cell * fe.n_dofs_per_face(face_no); while (i < fe.get_first_hex_index()) row_lengths[cell_indices[i++]] += increment; // Finally, cells in 3D - increment = fe.n_dofs_per_cell() - - GeometryInfo::faces_per_cell * fe.n_dofs_per_face(); + increment = fe.n_dofs_per_cell() - GeometryInfo::faces_per_cell * + fe.n_dofs_per_face(face_no); while (i < fe.n_dofs_per_cell()) row_lengths[cell_indices[i++]] += increment; @@ -225,7 +231,8 @@ namespace MGTools for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell(); ++local_dof) - row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face(); + row_lengths[cell_indices[local_dof]] += + fe.n_dofs_per_face(face_no); continue; } @@ -244,7 +251,7 @@ namespace MGTools if (flux_coupling != DoFTools::none) { const unsigned int dof_increment = - nfe.n_dofs_per_cell() - nfe.n_dofs_per_face(); + nfe.n_dofs_per_cell() - nfe.n_dofs_per_face(face_no); for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell(); ++local_dof) @@ -272,10 +279,12 @@ namespace MGTools neighbor->get_mg_dof_indices(neighbor_indices); for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell(); ++local_dof) - row_lengths[cell_indices[local_dof]] += nfe.n_dofs_per_face(); + row_lengths[cell_indices[local_dof]] += + nfe.n_dofs_per_face(face_no); for (unsigned int local_dof = 0; local_dof < nfe.n_dofs_per_cell(); ++local_dof) - row_lengths[neighbor_indices[local_dof]] += fe.n_dofs_per_face(); + row_lengths[neighbor_indices[local_dof]] += + fe.n_dofs_per_face(face_no); } } user_flags_triangulation.load_user_flags(old_flags); @@ -333,6 +342,12 @@ namespace MGTools const FiniteElement &fe = cell->get_fe(); const unsigned int fe_index = cell->active_fe_index(); + + // TODO: This assumes that the dofs per face on all faces coincide! + const unsigned int face_no = 0; + Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim), + ExcNotImplemented()); + Assert(couplings.n_rows() == fe.n_components(), ExcDimensionMismatch(couplings.n_rows(), fe.n_components())); Assert(couplings.n_cols() == fe.n_components(), @@ -378,8 +393,9 @@ namespace MGTools fe.first_block_of_base(base) + mult) != DoFTools::none) { - increment = fe.base_element(base).n_dofs_per_cell() - - dim * fe.base_element(base).n_dofs_per_face(); + increment = + fe.base_element(base).n_dofs_per_cell() - + dim * fe.base_element(base).n_dofs_per_face(face_no); row_lengths[cell_indices[i]] += increment; } ++i; @@ -394,7 +410,7 @@ namespace MGTools // In all other cases we // subtract adjacent faces to be // added in the loop below. - while (i < fe.get_first_quad_index()) + while (i < fe.get_first_quad_index(face_no)) { for (unsigned int base = 0; base < fe.n_base_elements(); ++base) for (unsigned int mult = 0; mult < fe.element_multiplicity(base); @@ -407,7 +423,7 @@ namespace MGTools fe.base_element(base).n_dofs_per_cell() - ((dim > 1) ? (dim - 1) : GeometryInfo::faces_per_cell) * - fe.base_element(base).n_dofs_per_face(); + fe.base_element(base).n_dofs_per_face(face_no); row_lengths[cell_indices[i]] += increment; } ++i; @@ -427,7 +443,7 @@ namespace MGTools fe.base_element(base).n_dofs_per_cell() - ((dim > 2) ? (dim - 2) : GeometryInfo::faces_per_cell) * - fe.base_element(base).n_dofs_per_face(); + fe.base_element(base).n_dofs_per_face(face_no); row_lengths[cell_indices[i]] += increment; } ++i; @@ -443,9 +459,10 @@ namespace MGTools fe.first_block_of_base(base) + mult) != DoFTools::none) { - increment = fe.base_element(base).n_dofs_per_cell() - - GeometryInfo::faces_per_cell * - fe.base_element(base).n_dofs_per_face(); + increment = + fe.base_element(base).n_dofs_per_cell() - + GeometryInfo::faces_per_cell * + fe.base_element(base).n_dofs_per_face(face_no); row_lengths[cell_indices[i]] += increment; } ++i; @@ -478,7 +495,8 @@ namespace MGTools for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell(); ++local_dof) - row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face(); + row_lengths[cell_indices[local_dof]] += + fe.n_dofs_per_face(face_no); continue; } @@ -506,7 +524,7 @@ namespace MGTools { const unsigned int dof_increment = nfe.base_element(base).n_dofs_per_cell() - - nfe.base_element(base).n_dofs_per_face(); + nfe.base_element(base).n_dofs_per_face(face_no); row_lengths[cell_indices[local_dof]] += dof_increment; } @@ -548,7 +566,7 @@ namespace MGTools fe.system_to_component_index(local_dof).first, nfe.first_block_of_base(base) + mult) != DoFTools::none) row_lengths[cell_indices[local_dof]] += - nfe.base_element(base).n_dofs_per_face(); + nfe.base_element(base).n_dofs_per_face(face_no); for (unsigned int base = 0; base < fe.n_base_elements(); ++base) for (unsigned int mult = 0; mult < fe.element_multiplicity(base); ++mult) @@ -559,7 +577,7 @@ namespace MGTools nfe.system_to_component_index(local_dof).first, fe.first_block_of_base(base) + mult) != DoFTools::none) row_lengths[neighbor_indices[local_dof]] += - fe.base_element(base).n_dofs_per_face(); + fe.base_element(base).n_dofs_per_face(face_no); } } user_flags_triangulation.load_user_flags(old_flags); @@ -1332,7 +1350,6 @@ namespace MGTools continue; const FiniteElement &fe = cell->get_fe(); const unsigned int level = cell->level(); - local_dofs.resize(fe.n_dofs_per_face()); for (const unsigned int face_no : GeometryInfo::face_indices()) if (cell->at_boundary(face_no) == true) @@ -1343,6 +1360,7 @@ namespace MGTools // Face is listed in boundary map if (boundary_ids.find(bi) != boundary_ids.end()) { + local_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_mg_dof_indices(level, local_dofs); dofs_by_level[level].insert(dofs_by_level[level].end(), local_dofs.begin(), @@ -1405,7 +1423,7 @@ namespace MGTools // get indices, physical location and boundary values of // dofs on this face - local_dofs.resize(fe.n_dofs_per_face()); + local_dofs.resize(fe.n_dofs_per_face(face_no)); face->get_mg_dof_indices(level, local_dofs); if (fe_is_system) { @@ -1415,7 +1433,8 @@ namespace MGTools numbers::invalid_unsigned_int; if (fe.is_primitive()) component = - fe.face_system_to_component_index(i).first; + fe.face_system_to_component_index(i, face_no) + .first; else { // Just pick the first of the components @@ -1473,7 +1492,6 @@ namespace MGTools const FiniteElement &fe = mg_dof_handler.get_fe(); const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); - const unsigned int dofs_per_face = fe.n_dofs_per_face(); std::vector local_dof_indices(dofs_per_cell); @@ -1517,7 +1535,8 @@ namespace MGTools // Do refinement face from the coarse side if (neighbor->level() < cell->level()) { - for (unsigned int j = 0; j < dofs_per_face; ++j) + for (unsigned int j = 0; j < fe.n_dofs_per_face(face_nr); + ++j) cell_dofs[fe.face_to_cell_index(j, face_nr)] = true; has_coarser_neighbor = true; -- 2.39.5