diff --git a/fem/fe/face_map_utils.cpp b/fem/fe/face_map_utils.cpp index 17e33b93ea..0d72407dd3 100644 --- a/fem/fe/face_map_utils.cpp +++ b/fem/fe/face_map_utils.cpp @@ -12,6 +12,7 @@ // Finite Element Base classes #include "face_map_utils.hpp" +#include // std::pow namespace mfem { @@ -57,4 +58,53 @@ void FillFaceMap(const int n_face_dofs_per_component, } } +void GetNodalTensorFaceMap(const int dim, const int order, const int face_id, + Array &face_map) +{ + const int dof1d = order + 1; + int n_face_dofs = int(std::pow(dof1d, dim - 1)); + std::vector offsets, strides; + switch (dim) + { + case 1: + offsets = {(face_id == 0) ? 0 : dof1d - 1}; + break; + case 2: + strides = {(face_id == 0 || face_id == 2) ? 1 : dof1d}; + switch (face_id) + { + case 0: offsets = {0}; break; // y = 0 + case 1: offsets = {dof1d - 1}; break; // x = 1 + case 2: offsets = {(dof1d-1)*dof1d}; break; // y = 1 + case 3: offsets = {0}; break; // x = 0 + } + break; + case 3: + { + const auto f = GetFaceNormal3D(face_id); + const int face_normal = f.first, level = f.second; + if (face_normal == 0) // x-normal + { + offsets = {level ? dof1d-1 : 0}; + strides = {dof1d, dof1d*dof1d}; + } + else if (face_normal == 1) // y-normal + { + offsets = {level ? (dof1d-1)*dof1d : 0}; + strides = {1, dof1d*dof1d}; + } + else if (face_normal == 2) // z-normal + { + offsets = {level ? (dof1d-1)*dof1d*dof1d : 0}; + strides = {1, dof1d}; + } + break; + } + } + + // same number of DOFs in each dimension, repeat dof1d (dim - 1) times + std::vector n_dofs(dim - 1, dof1d); + FillFaceMap(n_face_dofs, offsets, strides, n_dofs, face_map); } + +} // namespace mfem diff --git a/fem/fe/face_map_utils.hpp b/fem/fe/face_map_utils.hpp index 3ecb8d1dc7..dbaabdff24 100644 --- a/fem/fe/face_map_utils.hpp +++ b/fem/fe/face_map_utils.hpp @@ -44,6 +44,10 @@ void FillFaceMap(const int n_face_dofs_per_component, const std::vector &n_dofs_per_dim, Array &face_map); +/// Return the face map for nodal tensor elements (H1, L2, and Bernstein basis). +void GetNodalTensorFaceMap(const int dim, const int order, const int face_id, + Array &face_map); + } // namespace mfem #endif diff --git a/fem/fe/fe_base.cpp b/fem/fe/fe_base.cpp index 74eb2a9e26..94a8acfd9a 100644 --- a/fem/fe/fe_base.cpp +++ b/fem/fe/fe_base.cpp @@ -2519,50 +2519,7 @@ void NodalTensorFiniteElement::SetMapType(const int map_type) void NodalTensorFiniteElement::GetFaceMap(const int face_id, Array &face_map) const { - const int dof1d = order + 1; - int n_face_dofs = pow(dof1d, dim - 1); - std::vector offsets, strides; - switch (dim) - { - case 1: - offsets = {(face_id == 0) ? 0 : dof1d - 1}; - break; - case 2: - strides = {(face_id == 0 || face_id == 2) ? 1 : dof1d}; - switch (face_id) - { - case 0: offsets = {0}; break; // y = 0 - case 1: offsets = {dof1d - 1}; break; // x = 1 - case 2: offsets = {(dof1d-1)*dof1d}; break; // y = 1 - case 3: offsets = {0}; break; // x = 0 - } - break; - case 3: - { - const auto f = GetFaceNormal3D(face_id); - const int face_normal = f.first, level = f.second; - if (face_normal == 0) // x-normal - { - offsets = {level ? dof1d-1 : 0}; - strides = {dof1d, dof1d*dof1d}; - } - else if (face_normal == 1) // y-normal - { - offsets = {level ? (dof1d-1)*dof1d : 0}; - strides = {1, dof1d*dof1d}; - } - else if (face_normal == 2) // z-normal - { - offsets = {level ? (dof1d-1)*dof1d*dof1d : 0}; - strides = {1, dof1d}; - } - break; - } - } - - // same number of DOFs in each dimension, repeat dof1d (dim - 1) times - std::vector n_dofs(dim - 1, dof1d); - FillFaceMap(n_face_dofs, offsets, strides, n_dofs, face_map); + GetNodalTensorFaceMap(dim, order, face_id, face_map); } VectorTensorFiniteElement::VectorTensorFiniteElement(const int dims, diff --git a/fem/fe/fe_pos.cpp b/fem/fe/fe_pos.cpp index ad5dac5d2f..44ea937635 100644 --- a/fem/fe/fe_pos.cpp +++ b/fem/fe/fe_pos.cpp @@ -12,6 +12,7 @@ // H1 Finite Element classes utilizing the Bernstein basis #include "fe_pos.hpp" +#include "face_map_utils.hpp" #include "../bilininteg.hpp" #include "../lininteg.hpp" #include "../coefficient.hpp" @@ -84,6 +85,12 @@ PositiveTensorFiniteElement::PositiveTensorFiniteElement( dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk), TensorBasisElement(dims, p, BasisType::Positive, dmtype) { } +void PositiveTensorFiniteElement::GetFaceMap(const int face_id, + Array &face_map) const +{ + GetNodalTensorFaceMap(dim, order, face_id, face_map); +} + BiQuadPos2DFiniteElement::BiQuadPos2DFiniteElement() : PositiveFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk) diff --git a/fem/fe/fe_pos.hpp b/fem/fe/fe_pos.hpp index 4af7e53b11..a1d102a567 100644 --- a/fem/fe/fe_pos.hpp +++ b/fem/fe/fe_pos.hpp @@ -70,12 +70,14 @@ public: const DofMapType dmtype); const DofToQuad &GetDofToQuad(const IntegrationRule &ir, - DofToQuad::Mode mode) const + DofToQuad::Mode mode) const override { return (mode == DofToQuad::FULL) ? FiniteElement::GetDofToQuad(ir, mode) : GetTensorDofToQuad(*this, ir, mode, basis1d, true, dof2quad_array); } + + virtual void GetFaceMap(const int face_id, Array &face_map) const override; };