From: Peter Munch Date: Sun, 8 May 2022 15:29:24 +0000 (+0200) Subject: Generalize do_interpolate_boundary_values() for pyramid and wedge X-Git-Tag: v9.4.0-rc1~236^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=059deb72e13adff0a13eb1ad5a96cb9da55a7cd5;p=dealii.git Generalize do_interpolate_boundary_values() for pyramid and wedge --- diff --git a/include/deal.II/numerics/vector_tools_boundary.templates.h b/include/deal.II/numerics/vector_tools_boundary.templates.h index ffde7f2f6b..bc6f01ba6c 100644 --- a/include/deal.II/numerics/vector_tools_boundary.templates.h +++ b/include/deal.II/numerics/vector_tools_boundary.templates.h @@ -160,60 +160,65 @@ 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 const dealii::hp::FECollection &finite_elements = dof.get_fe_collection(); - dealii::hp::QCollection q_collection; + std::vector> q_collection( + finite_elements.size()); for (unsigned int f = 0; f < finite_elements.size(); ++f) - { - const FiniteElement &fe = finite_elements[f]; - - // generate a quadrature rule on the face from the unit support - // points. this will be used to obtain the quadrature points on - // the real cell's face - // - // to do this, we check whether the FE has support points on the - // face at all: - 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 - // that a finite element may not offer support points for all - // its shape functions, but maybe only some. if it offers - // support points for the components we are interested in in - // this function, then that's fine. if not, the function we - // call in the finite element will raise an exception. the - // support points for the other shape functions are left - // uninitialized (well, initialized by the default - // constructor), since we don't need them anyway. - // - // As a detour, we must make sure we only query - // face_system_to_component_index if the index corresponds to - // a primitive shape function. since we know that all the - // 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(face_no)); - - 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, face_no); - - q_collection.push_back( - Quadrature(unit_support_points)); - } - } + for (unsigned int face_no = 0; + face_no < (finite_elements[f].n_unique_faces() == 1 ? + 1 : + finite_elements[f].reference_cell().n_faces()); + ++face_no) + { + const FiniteElement &fe = finite_elements[f]; + + // generate a quadrature rule on the face from the unit support + // points. this will be used to obtain the quadrature points on + // the real cell's face + // + // to do this, we check whether the FE has support points on the + // face at all: + if (fe.has_face_support_points(face_no)) + q_collection[f].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 + // that a finite element may not offer support points for + // all its shape functions, but maybe only some. if it + // offers support points for the components we are + // interested in in this function, then that's fine. if not, + // the function we call in the finite element will raise an + // exception. the support points for the other shape + // functions are left uninitialized (well, initialized by + // the default constructor), since we don't need them + // anyway. + // + // As a detour, we must make sure we only query + // face_system_to_component_index if the index corresponds + // to a primitive shape function. since we know that all the + // 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(face_no)); + + 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, face_no); + + q_collection[f].push_back( + Quadrature(unit_support_points)); + } + } // now that we have a q_collection object with all the right // quadrature points, create an hp::FEFaceValues object that we can // use to evaluate the boundary values at @@ -283,16 +288,9 @@ namespace VectorTools if (fe_is_system) { - // resize array. avoid construction of a memory - // allocating temporary if possible - 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(face_no)); + dof_values_system.resize(fe.n_dofs_per_face(face_no), + Vector( + fe.n_components())); function_map.find(boundary_component) ->second->vector_value_list(dof_locations,