diff --git a/CMakeLists.txt b/CMakeLists.txt index b385901..4635e68 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -220,6 +220,10 @@ meshfields_add_exe(IntegratorTest test/testIntegrator.cpp) meshfields_add_exe(BoxIntegratorTest test/testBoxIntegrator.cpp) meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) +meshfields_add_exe(LinearCubicWedgeShapeTest test/testLinearCubicWedgeShape.cpp) +meshfields_add_exe(LinearCubicWedgeElementTest test/testLinearCubicWedgeElement.cpp) +meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) +meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) meshfields_add_exe(Interpolate test/testInterpolate.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -246,6 +250,10 @@ test_func(BoxIntegratorTest ./BoxIntegratorTest) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) +test_func(LinearCubicWedgeShapeTest ./LinearCubicWedgeShapeTest) +test_func(LinearCubicWedgeElementTest ./LinearCubicWedgeElementTest) +test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest) +test_func(LinearWedgeElementTest ./LinearWedgeElementTest) test_func(Interpolate ./Interpolate) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error diff --git a/src/MeshField_Defines.hpp b/src/MeshField_Defines.hpp index 8467a0b..1c79020 100644 --- a/src/MeshField_Defines.hpp +++ b/src/MeshField_Defines.hpp @@ -13,7 +13,8 @@ enum Mesh_Topology { Tetrahedron, Hexahedron, Prism, - Pyramid + Pyramid, + Wedge }; static bool Debug = false; const Real MachinePrecision = 1e-15; diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 6fd4d47..419c70e 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -557,5 +557,276 @@ struct QuadraticTetrahedronShape { } }; +/** + * @brief Linear shape functions for 3D wedge (triangular prism) elements + * + * Defines linear basis functions for a wedge element with 6 nodes: the tensor product + * of the linear (P1) triangle in (xi0, xi1) and the linear edge in + * zeta = xi2. + * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with + * xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1] + * + * Node ordering (node a+3 is joined to node a by an edge along \zeta): + * - Node 0 at (0, 0 , -1) + * - Node 1 at (1, 0 , -1) + * - Node 2 at (0, 1 , -1) + * - Node 3 at (0, 0 , 1) + * - Node 4 at (1, 0 , 1) + * - Node 5 at (0, 1 , 1) + * + * Shape functions use barycentric coordinates L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 of the triangle and, for a = 0, 1, 2: + * - \f$N_a = L_a (1 - \zeta) / 2\f$ + * - \f$N_{a+3} = L_a (1 + \zeta) / 2\f$ + * + * Shape functions and ordering from: + * Zienkiewicz, Taylor, and Zhu + * 'The Finite Element Method: Its Basis and Fundamentals', 2013, Section 6.4, + * Fig. 6.16a: \f$N_a = \frac{1}{2} L_a (1 + \zeta_a \zeta)\f$, with + * \f$\zeta_a = -1\f$ for nodes 0-2 and \f$+1\f$ for nodes 3-5 + */ +struct LinearWedgeShape { + // clang-format off + static const size_t numNodes = 6; ///< Number of nodes (6 for linear wedge) + static const size_t meshEntDim = 3; ///< Mesh entity dimension (3D) + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) + // clang-format on + + /** + * @brief Get parametric coordinates of element nodes + * @return Array of node coordinates [node0_xi0, node0_xi1, node0_zeta, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getNodeParametricCoords() const { + // clang-format off + return {0,0,-1, //node 0 + 1,0,-1, //node 1 + 0,1,-1, //node 2 + 0,0, 1, //node 3 + 1,0, 1, //node 4 + 0,1, 1}; //node 5 + // clang-format on + } + + /** + * @brief Evaluate shape functions at parametric coordinate + * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1 + * \f$\in\f$ [0,1], xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1] + * @return Array of shape function values [N0, ..., N5] where + * N_a = L_a ( 1 - \zeta)/2 and N_{a+3} = L_a (1+\zeta)/2 for a = 0,1,2 + * with L0 = 1 - x0 - xi1, L1 = xi0 and L2 = xi1 + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getValues(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + const Real down = (1 - xi[2]) / 2; // linear edge function of the -1 face + const Real up = (1 + xi[2]) / 2; // linear edge function of the +1 face + // clang-format off + return {L0*down, //node 0 + xi[0]*down, //node 1 + xi[1]*down, //node 2 + L0*up, //node 3 + xi[0]*up, //node 4 + xi[1]*up}; //node 5 + // clang-format on + } + + /** + * @brief Get gradients of shape functions in parametric coordinates + * @param xi Parametric coordinates (xi0, xi1, zeta), see getValues + * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN0/dzeta, dN1/dxi0, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array + getLocalGradients(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + const Real L1 = xi[0]; + const Real L2 = xi[1]; + const Real down = (1 - xi[2]) / 2; + const Real up = (1 + xi[2]) / 2; + // clang-format off + return {-down, -down, -L0/2, //node 0 + down, 0, -L1/2, //node 1 + 0, down, -L2/2, //node 2 + -up, -up, L0/2, //node 3 + up, 0, L1/2, //node 4 + 0, up, L2/2}; //node 5 + // clang-format on + } +}; + + +/** + * @brief Linear in the triangle and cubic in the zeta direction shape functions for 3D tensor product wedge elements + * + * Defines shape functions for a wedge element with 12 nodes: the tensor product + * of the linear triangle basis in (xi0, xi1) and the cubic lagrange basis in + * zeta = xi2. + * + * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with + * xi0, xi1 \f$\ge\f$ 0, xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1] + * + * Node ordering (nodes 0-5 as in LinearWedgeShape): + * - Node 0-2: triangle vertices at zeta = -1 + * - Node 3-5: triangle vertices at zeta = +1 (node a + 3 is directly above node a) + * - Nodes 6 + 2a and 7+2a , a = 0, 1, 2 : on the vertical edge through vertex a, + * at zeta = -1/3 and zeta = +1/3 + * + * Barycentric coordinates of the triangle: L0 = 1 - xi0 - xi1, L1 = xi0, L2 = xi1. + * Cubic Lagrange functions in zeta on the stations -1, -1/3, 1/3, 1: + * \f[ + * \ell_0 = -\tfrac{9}{16}(\zeta+\tfrac{1}{3})(\zeta-\tfrac{1}{3})(\zeta-1), \qquad + * \ell_1 = \tfrac{27}{16}(\zeta+1)(\zeta-\tfrac{1}{3})(\zeta-1), + * \f] + * \f[ + * \ell_2 = -\tfrac{27}{16}(\zeta+1)(\zeta+\tfrac{1}{3})(\zeta-1), \qquad + * \ell_3 = \tfrac{9}{16}(\zeta+1)(\zeta+\tfrac{1}{3})(\zeta-\tfrac{1}{3}). + * \f] + * Shape functions, for a = 0, 1, 2: + * - \f$N_a = L_a \, \ell_0\f$ and \f$N_{a+3} = L_a \, \ell_3\f$ + * - \f$N_{6+2a} = L_a \, \ell_1\f$ and \f$N_{7+2a} = L_a \, \ell_2\f$ + * + */ + +struct LinearCubicWedgeShape { + // clang-format off + static const size_t numNodes = 12; ///< Number of nodes (12 for linear-cubic wedge) + static const size_t meshEntDim = 3; ///< Mesh entity dimension (3D) + constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge}; ///< DOFs at vertices and edges + constexpr static size_t NumDofHolders[2] = {6, 3}; ///< 6 vertices, 3 edges + constexpr static size_t DofsPerHolder[2] = {1, 2}; ///< 1 DOF per vertex, 2 DOF per edge + constexpr static size_t Order = 3; ///< Polynomial order (cubic) + // clang-format on + + /** + * @brief Get parametric coordinates of element nodes + * @return Array of node coordinates [node0_xi0, node0_xi1, node0_zeta, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getNodeParametricCoords() const { + // clang-format off + return { + // nodes at vertices + 0,0,-1, //node 0 + 1,0,-1, //node 1 + 0,1,-1, //node 2 + 0,0, 1, //node 3 + 1,0, 1, //node 4 + 0,1, 1, //node 5 + + // nodes at two stations of edges at -1/3 and +1/3 + 0, 0, -1.0/3, // node 6 + 0, 0, 1.0/3, // node 7 + 1, 0, -1.0/3, // node 8 + 1, 0, 1.0/3, // node 9 + 0, 1, -1.0/3, // node 10 + 0, 1, 1.0/3 // node 11 + }; + // clang-format on + } + + /** + * @brief Evaluate shape functions at parametric coordinate + * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1 + * \f$\ge\f$ 0, xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1] + * @return Array of shape function values [N0, ..., N11] where, for a = 0, 1, 2, + * \f$N_a = L_a \, \ell_0(\zeta)\f$, \f$N_{a+3} = L_a \, \ell_3(\zeta)\f$, + * \f$N_{6+2a} = L_a \, \ell_1(\zeta)\f$ and \f$N_{7+2a} = L_a \, \ell_2(\zeta)\f$, + * with L0 = 1 - xi0 - xi1, L1 = xi0, L2 = xi1 and the cubic Lagrange + * functions \f$\ell_0, \dots, \ell_3\f$ defined in the class description + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getValues(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + const Real L1 = xi[0]; + const Real L2 = xi[1]; + + const Real l0 = -1.0/16 * (3 * xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); //cubic edge function at -1 station + const Real l1 = 9.0/16 * (xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); // cubic edge function at -1/3 station + const Real l2 = -9.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (xi[2] - 1); // cubic edge function at +1/3 station + const Real l3 = 1.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (3 * xi[2] - 1); // cubic edge function at +1 station + + // clang-format off + return {L0 * l0, + L1 * l0, + L2 * l0, + L0 * l3, + L1 * l3, + L2 * l3, + L0 * l1, + L0 * l2, + L1 * l1, + L1 * l2, + L2 * l1, + L2 * l2 + }; + // clang-format on + } + + /** + * @brief Get gradients of shape functions in parametric coordinates + * @param xi Parametric coordinates (xi0, xi1, zeta), see getValues + * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN0/dzeta, dN1/dxi0, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array + getLocalGradients(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + const Real L1 = xi[0]; + const Real L2 = xi[1]; + + const Real l0 = -1.0/16 * (3 * xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); //cubic edge function at -1 station + const Real dl0 = 1.0/16 * (-27 * xi[2] * xi[2] + 18 * xi[2] + 1); //derivative of cubic edge function at -1 station + + const Real l1 = 9.0/16 * (xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); // cubic edge function at -1/3 station + const Real dl1 = 1.0/16 * (81 * xi[2] * xi[2] - 18 * xi[2] -27);// derivative of cubic edge function at -1/3 station + + const Real l2 = -9.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (xi[2] - 1); // cubic edge function at +1/3 station + const Real dl2 = 1.0/16 * (-81 * xi[2] * xi[2] - 18 * xi[2] + 27); // derivative of cubic edge function at +1/3 station + + const Real l3 = 1.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (3 * xi[2] - 1); // cubic edge function at +1 station + const Real dl3 = 1.0/16 * (27 * xi[2] * xi[2] + 18 * xi[2] - 1); // derivative of cubic edge function at +1 station + + + // clang-format off + return {-l0, -l0, L0 * dl0, //node 0 + l0, 0, L1 * dl0, //node 1 + 0, l0, L2 * dl0, //node 2 + -l3, -l3, L0 * dl3, //node 3 + l3, 0, L1 * dl3, //node 4 + 0, l3, L2 * dl3, //node 5 + // nodes at two stations of edges at -1/3 and +1/3 + -l1, -l1, L0 * dl1, // node 6 + -l2, -l2, L0 * dl2, // node 7 + l1, 0 , L1 * dl1, // node 8 + l2, 0 , L1 * dl2, // node 9 + 0, l1, L2 * dl1, // node 10 + 0, l2 , L2 * dl2 // node 11 + }; + // clang-format on + } +}; + } // namespace MeshField #endif diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 708e80a..c7b0b38 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -25,6 +25,7 @@ struct MeshInfo { int numHex = 0; // entDim = 3 int numPrism = 0; // entDim = 3 int numPyramid = 0; // entDim = 3 + int numWedge = 0; // entDim = 3 int dim = 0; /** @brief runtime map of number of entities based on topology */ @@ -46,6 +47,8 @@ struct MeshInfo { return numPrism; case Pyramid: return numPyramid; + case Wedge: + return numWedge; } return 0; } @@ -181,6 +184,18 @@ template struct LinearAccessor { }; //! [LinearAccessor] +/** + * @brief + * Shape functions used by CreateLagrangeField when the caller does not name + * one: triangles for 1d and 2d meshes and tetrahedra for 3d meshes. + */ +template +using DefaultLagrangeShape = std::conditional_t< + order == 2, + std::conditional_t, + std::conditional_t>; + /** * @brief * Create a field using linear or quadratic Lagrange shape functions @@ -197,6 +212,9 @@ template struct LinearAccessor { * @tparam order the order of the shape functions; linear(1) and quadratic(2) * are supported * @tparam dim the dimension of the mesh + * @tparam Shape the shape functions; the default picks the triangle or + * tetrahedron shape, and LinearCubicWedgeShape (order 3) stores one DOF per + * vertex and two per vertical edge * * @param meshInfo defines on-process mesh metadata * @return a linear or quadratic ShapeField @@ -204,19 +222,75 @@ template struct LinearAccessor { template typename Controller = MeshField::KokkosController, - typename DataType, size_t order, size_t dim, size_t numComp> + typename DataType, size_t order, size_t dim, size_t numComp, + typename Shape = DefaultLagrangeShape> auto CreateLagrangeField(const MeshInfo &meshInfo) { + constexpr bool isLinearCubicWedge = + std::is_same_v; static_assert((std::is_same_v == true || std::is_same_v == true), "CreateLagrangeField only supports single and double precision " "floating point fields\n"); - static_assert( - (order == 1 || order == 2), - "CreateLagrangeField only supports linear and quadratic fields\n"); + static_assert((order == 1 || order == 2 || isLinearCubicWedge), + "CreateLagrangeField only supports linear and quadratic " + "fields, and LinearCubicWedgeShape\n"); + static_assert(isLinearCubicWedge || + std::is_same_v>, + "CreateLagrangeField: Shape must be the default or " + "LinearCubicWedgeShape\n"); + static_assert(Shape::Order == order, + "CreateLagrangeField order must match Shape::Order\n"); static_assert((dim == 1 || dim == 2 || dim == 3), "CreateLagrangeField only supports 1d, 2d, and 3d meshes\n"); using MemorySpace = typename ExecutionSpace::memory_space; - if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) { + if constexpr (isLinearCubicWedge) { + static_assert(dim == 3, "LinearCubicWedgeShape fields need a 3d mesh\n"); + if (meshInfo.numVtx <= 0) { + fail("mesh has no vertices\n"); + } + if (meshInfo.numEdge <= 0) { + fail("mesh has no edges\n"); + } + if (meshInfo.numWedge <= 0) { + fail("mesh has no wedges\n"); + } + // 1 dof per vertex and Shape::DofsPerHolder[1] (2) dofs per vertical edge + constexpr int dofsPerEdge = Shape::DofsPerHolder[1]; +#ifdef MESHFIELDS_ENABLE_CABANA + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; + auto createController = [](auto numVtx, auto numEdge) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr({numVtx, numEdge}); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp, + /*field 1*/ numEdge, dofsPerEdge, numComp}); + } + }; + Ctrlr kk_ctrl = createController(meshInfo.numVtx, meshInfo.numEdge); +#else + using Ctrlr = + Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp, + /*field 1*/ meshInfo.numEdge, dofsPerEdge, numComp}); +#endif + auto vtxField = MeshField::makeField(kk_ctrl); + auto edgeField = MeshField::makeField(kk_ctrl); + // QuadraticAccessor forwards the dof index to the edge field, so it also + // serves edges with more than one dof + using QA = QuadraticAccessor; + using WedgeShapeField = ShapeField; + WedgeShapeField wsf(meshInfo, {vtxField, edgeField}); + return FieldWithController{kk_ctrl, wsf}; + } else if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } diff --git a/test/testLinearCubicWedgeElement.cpp b/test/testLinearCubicWedgeElement.cpp new file mode 100644 index 0000000..0ec2c17 --- /dev/null +++ b/test/testLinearCubicWedgeElement.cpp @@ -0,0 +1,220 @@ +#include "MeshField_Defines.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Shape.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include +#include +#include +#include +#include + +// Checks MeshField::LinearCubicWedgeShape inside a MeshField::FieldElement: +// field evaluation (MeshField::evaluate) and the Jacobian +// (MeshField::getJacobians, MeshField::getJacobianDeterminants) for two cubic +// wedges. The coordinates are stored in a field from +// MeshField::CreateLagrangeField with LinearCubicWedgeShape: one DOF per vertex +// and two per vertical edge. +// The nodes are placed by geometric maps that are linear in (xi0, xi1) +// and cubic in zeta, so the element map reproduces them exactly and the +// expected values are the analytic maps and their derivatives. Failures are +// counted explicitly so the test also fails in builds where assert() is +// disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Vec3 = Omega_h::Vector<3>; + +namespace { + +constexpr int numWedges = 2; +constexpr int numNodes = MeshField::LinearCubicWedgeShape::numNodes; +constexpr int numVtxPerWedge = 6; +constexpr int numEdgesPerWedge = 3; + +bool isClose(Real a, Real b) { return Omega_h::are_close(a, b, 1e-12, 1e-12); } + +// Local node n of wedge w: nodes 0-5 are vertex 6*w + n; nodes 6+2a and 7+2a +// are the two DOFs on vertical edge 3*w + a, at zeta = -1/3 and +1/3. +struct CubicWedgeToField { + static constexpr KOKKOS_FUNCTION Kokkos::Array + getTopology() { + return {MeshField::Wedge}; + } + + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap + operator()(MeshField::LO wedgeNodeIdx, MeshField::LO wedgeCompIdx, + MeshField::LO wedge, MeshField::Mesh_Topology topo) const { + assert(topo == MeshField::Wedge); + if (wedgeNodeIdx < numVtxPerWedge) { + const MeshField::LO vtx = numVtxPerWedge * wedge + wedgeNodeIdx; + return {0, wedgeCompIdx, vtx, MeshField::Vertex}; + } + const MeshField::LO a = (wedgeNodeIdx - numVtxPerWedge) / 2; + const MeshField::LO dof = (wedgeNodeIdx - numVtxPerWedge) % 2; + const MeshField::LO edge = numEdgesPerWedge * wedge + a; + return {dof, wedgeCompIdx, edge, MeshField::Edge}; + } +}; + +// Geometric maps that place the nodes of each wedge. Both are linear in +// (xi0, xi1) and cubic in zeta, so they lie in the element's function space. +Vec3 mapX(int w, Vec3 xi) { + const Real x0 = xi[0], x1 = xi[1], z = xi[2]; + if (w == 0) + return {2 * x0 + 0.3 * z * z * x1 + 0.2 * z * z * z, + 3 * x1 + 0.1 * z * z * z * x0 - 0.4 * z, + 2 + 2 * z + 0.25 * z * z * z}; + return {x0 + 0.5 * x1 + 0.3 * z * z * z * x0 + 0.1 * z, + x1 + 0.2 * z * z * x1 + 0.4 * z * z, + 1.5 + z + 0.2 * z * z * z * x1 + 0.1 * z * z * z}; +} + +// Analytic derivatives of mapX. Row i is dx/dxi_i, the same layout as the +// MeshField Jacobian: J(i, j) = dx_j / dxi_i. +std::array mapJ(int w, Vec3 xi) { + const Real x0 = xi[0], x1 = xi[1], z = xi[2]; + if (w == 0) + return {Vec3{2, 0.1 * z * z * z, 0}, Vec3{0.3 * z * z, 3, 0}, + Vec3{0.6 * z * x1 + 0.6 * z * z, 0.3 * z * z * x0 - 0.4, + 2 + 0.75 * z * z}}; + return {Vec3{1 + 0.3 * z * z * z, 0, 0}, + Vec3{0.5, 1 + 0.2 * z * z, 0.2 * z * z * z}, + Vec3{0.9 * z * z * x0 + 0.1, 0.4 * z * x1 + 0.8 * z, + 1 + 0.6 * z * z * x1 + 0.3 * z * z}}; +} + +// Reference coordinates of the 12 nodes: six vertices, then two edge nodes on +// each vertical edge at zeta = -1/3 and +1/3. +const Vec3 nodeRef[numNodes] = { + {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, + {1, 0, 1}, {0, 1, 1}, {0, 0, -1. / 3}, {0, 0, 1. / 3}, + {1, 0, -1. / 3}, {1, 0, 1. / 3}, {0, 1, -1. / 3}, {0, 1, 1. / 3}}; + +// The 12 nodes, then interior points. +std::vector makeLocalPoints() { + std::vector pts(nodeRef, nodeRef + numNodes); + pts.push_back({1. / 3, 1. / 3, 0}); + pts.push_back({0.2, 0.3, -0.4}); + pts.push_back({0.6, 0.1, 0.7}); + pts.push_back({0.05, 0.9, 0.95}); + return pts; +} + +// Writes the node positions from mapX into the coordinate field: vertex +// nodes to the vertex DOFs, and the zeta = -1/3 and +1/3 nodes to DOFs 0 and 1 +// of their vertical edge. +template void setCubicWedgeCoords(CoordField field) { + const int numVtx = numWedges * numVtxPerWedge; + const int numEdges = numWedges * numEdgesPerWedge; + Kokkos::View vtx_h("vtx_h", numVtx); + Kokkos::View edge_h("edge_h", numEdges); + + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < numNodes; ++n) { + const Vec3 X = mapX(w, nodeRef[n]); + for (int d = 0; d < 3; ++d) { + if (n < numVtxPerWedge) { + vtx_h(numVtxPerWedge * w + n, d) = X[d]; + } else { + const int a = (n - numVtxPerWedge) / 2; + const int dof = (n - numVtxPerWedge) % 2; + edge_h(numEdgesPerWedge * w + a, dof, d) = X[d]; + } + } + } + + auto vtx = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), vtx_h); + auto edge = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), edge_h); + Kokkos::parallel_for( + "setVtxCoords", Kokkos::RangePolicy(0, numVtx), + KOKKOS_LAMBDA(const int v) { + for (int d = 0; d < 3; ++d) + field(v, 0, d, MeshField::Vertex) = vtx(v, d); + }); + Kokkos::parallel_for( + "setEdgeCoords", Kokkos::RangePolicy(0, numEdges), + KOKKOS_LAMBDA(const int e) { + for (int dof = 0; dof < 2; ++dof) + for (int d = 0; d < 3; ++d) + field(e, dof, d, MeshField::Edge) = edge(e, dof, d); + }); +} + +// Compares the evaluated positions, Jacobians and determinants of wedge w +// with the analytic maps at every point in pts. Returns the failure count. +template +int checkWedge(int w, const std::vector &pts, const XView &x_h, + const JView &J_h, const DetView &detJ_h) { + const int numPts = static_cast(pts.size()); + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + const Vec3 xi = pts[p]; + const Vec3 xExpected = mapX(w, xi); + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + const auto jExpected = mapJ(w, xi); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); + const Real detExpected = + jExpected[0] * Omega_h::cross(jExpected[1], jExpected[2]); + detFailed += !isClose(detJ_h(row), detExpected); + } + std::cout << "LinearCubicWedgeShape FieldElement, cubic wedge " << w + << ": evaluate failures " << evalFailed << ", Jacobian failures " + << jacFailed << ", det J failures " << detFailed << "\n"; + return evalFailed + jacFailed + detFailed; +} + +// Builds the field and the FieldElement, runs the evaluation, Jacobian and +// determinant, and checks every wedge. +int runCubicWedgeElementChecks() { + const auto pts = makeLocalPoints(); + const int numPts = static_cast(pts.size()); + // numEdge counts only the vertical edges, which hold the zeta = +-1/3 DOFs + const MeshField::MeshInfo meshInfo{.numVtx = numWedges * numVtxPerWedge, + .numEdge = numWedges * numEdgesPerWedge, + .numWedge = numWedges, + .dim = 3}; + auto coordField = + MeshField::CreateLagrangeField( + meshInfo); + setCubicWedgeCoords(coordField.field); + MeshField::FieldElement fe(numWedges, coordField.field, + MeshField::LinearCubicWedgeShape(), + CubicWedgeToField()); + // the same local points in every wedge, stored element-major + Kokkos::View lc("localCoords", numWedges * numPts, 3); + auto lc_h = Kokkos::create_mirror_view(lc); + for (int w = 0; w < numWedges; ++w) + for (int p = 0; p < numPts; ++p) + for (int d = 0; d < 3; ++d) + lc_h(w * numPts + p, d) = pts[p][d]; + Kokkos::deep_copy(lc, lc_h); + + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + const auto J_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + + int numFailed = 0; + for (int w = 0; w < numWedges; ++w) + numFailed += checkWedge(w, pts, x_h, J_h, detJ_h); + return numFailed; +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runCubicWedgeElementChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} diff --git a/test/testLinearCubicWedgeShape.cpp b/test/testLinearCubicWedgeShape.cpp new file mode 100644 index 0000000..8b37471 --- /dev/null +++ b/test/testLinearCubicWedgeShape.cpp @@ -0,0 +1,147 @@ +#include +#include +#include +#include + +// Checks of MeshField::LinearCubicWedgeShape, linear in (xi0, xi1) and cubic +// in zeta. The checks run inside Kokkos kernels so the shape functions are +// exercised in the default execution space, and failures are counted +// explicitly so the test also fails in builds where assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Shape = MeshField::LinearCubicWedgeShape; + +KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { + return Kokkos::fabs(a - b) <= tol; +} + +// A function in the span of the cubic wedge basis, +// {1, xi0, xi1} x {1, zeta, zeta^2, zeta^3}, and its gradient +KOKKOS_INLINE_FUNCTION Real f(Real x0, Real x1, Real z) { + return 0.5 + 2 * x0 - 3 * x1 + 0.25 * z - 1.5 * z * z + 0.75 * z * z * z + + 4 * x0 * z - 1.5 * x1 * z * z - 2 * x0 * z * z * z + x1 * z * z * z; +} +KOKKOS_INLINE_FUNCTION MeshField::Vector3 gradF(Real x0, Real x1, Real z) { + return {2 + 4 * z - 2 * z * z * z, -3 - 1.5 * z * z + z * z * z, + 0.25 - 3 * z + 2.25 * z * z + 4 * x0 - 3 * x1 * z - 6 * x0 * z * z + + 3 * x1 * z * z}; +} + +// Each shape function is one at its own node and zero at the others. +int checkKroneckerDelta() { + int numFailed = 0; + Kokkos::parallel_reduce( + "cubicWedgeKronecker", + Kokkos::RangePolicy(0, Shape::numNodes), + KOKKOS_LAMBDA(const int node, int &failed) { + const Shape shp; + const auto coords = shp.getNodeParametricCoords(); + const MeshField::Vector3 xi{coords[3 * node], coords[3 * node + 1], + coords[3 * node + 2]}; + const auto N = shp.getValues(xi); + for (size_t m = 0; m < Shape::numNodes; ++m) { + const Real expected = (static_cast(node) == m) ? 1.0 : 0.0; + failed += !isClose(N[m], expected, 10 * MeshField::MachinePrecision); + } + }, + numFailed); + return numFailed; +} + +// Interior points, including points at the zeta stations and between them. +Kokkos::View makeInteriorPoints() { + // clang-format off + const Real pts[][3] = {{1. / 3, 1. / 3, 0.0}, + {0.2, 0.3, -0.4}, + {0.6, 0.1, 0.7}, + {0.05, 0.9, 0.95}, + {0.45, 0.45, -0.95}, + {0.2, 0.3, 1. / 3}, // at a station + {0.6, 0.1, -1. / 3}, // at a station + {0.3, 0.2, -0.7}, // between -1 and -1/3 + {0.1, 0.8, 0.5}}; // between 1/3 and 1 + // clang-format on + const int numPts = static_cast(sizeof(pts) / sizeof(pts[0])); + Kokkos::View pts_d("cubicWedgePts", numPts); + auto pts_h = Kokkos::create_mirror_view(pts_d); + for (int i = 0; i < numPts; ++i) + for (int d = 0; d < 3; ++d) + pts_h(i, d) = pts[i][d]; + Kokkos::deep_copy(pts_d, pts_h); + return pts_d; +} + +// At interior points: the values sum to one, the gradients sum to zero, the +// basis reproduces f and its gradient, and the analytic gradients match +// central finite differences of the values. +int checkInteriorPoints() { + const auto pts_d = makeInteriorPoints(); + const int numPts = static_cast(pts_d.extent(0)); + + int numFailed = 0; + Kokkos::parallel_reduce( + "cubicWedgeInterior", Kokkos::RangePolicy(0, numPts), + KOKKOS_LAMBDA(const int i, int &failed) { + const Real tol = 10 * MeshField::MachinePrecision; + const Shape shp; + const MeshField::Vector3 xi{pts_d(i, 0), pts_d(i, 1), pts_d(i, 2)}; + const auto N = shp.getValues(xi); + const auto G = shp.getLocalGradients(xi); + const auto nodes = shp.getNodeParametricCoords(); + + Real sumN = 0; + Real sumG[3] = {0, 0, 0}; + Real fInterp = 0; + Real gInterp[3] = {0, 0, 0}; + for (size_t n = 0; n < Shape::numNodes; ++n) { + const Real fn = f(nodes[3 * n], nodes[3 * n + 1], nodes[3 * n + 2]); + sumN += N[n]; + fInterp += N[n] * fn; + for (int d = 0; d < 3; ++d) { + sumG[d] += G[3 * n + d]; + gInterp[d] += G[3 * n + d] * fn; + } + } + failed += !isClose(sumN, 1.0, tol); + failed += !isClose(fInterp, f(xi[0], xi[1], xi[2]), tol); + const auto gExact = gradF(xi[0], xi[1], xi[2]); + for (int d = 0; d < 3; ++d) { + failed += !isClose(sumG[d], 0.0, tol); + failed += !isClose(gInterp[d], gExact[d], tol); + } + + // gradients against central finite differences of the values + const Real h = 1e-6; + for (int d = 0; d < 3; ++d) { + auto xp = xi; + auto xm = xi; + xp[d] += h; + xm[d] -= h; + const auto Np = shp.getValues(xp); + const auto Nm = shp.getValues(xm); + for (size_t n = 0; n < Shape::numNodes; ++n) + failed += !isClose((Np[n] - Nm[n]) / (2 * h), G[3 * n + d], 1e-8); + } + }, + numFailed); + return numFailed; +} + +// Runs both checks and prints their failure counts. +int runCubicWedgeShapeChecks() { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearCubicWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearCubicWedgeShape interior point failures: " << interior + << "\n"; + return kronecker + interior; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runCubicWedgeShapeChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} diff --git a/test/testLinearWedgeElement.cpp b/test/testLinearWedgeElement.cpp new file mode 100644 index 0000000..057a51a --- /dev/null +++ b/test/testLinearWedgeElement.cpp @@ -0,0 +1,193 @@ +#include "KokkosController.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Field.hpp" +#include "MeshField_For.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include +#include +#include +#include + +// Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field +// evaluation (MeshField::evaluate) and the Jacobian (MeshField::getJacobians, +// MeshField::getJacobianDeterminants) for two wedges, compared against the +// geometric map x(xi) differentiated by hand. Failures are counted explicitly +// so the test also fails in builds where assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Vec3 = Omega_h::Vector<3>; + +constexpr int numWedges = 2; +constexpr int numNodes = MeshField::LinearWedgeShape::numNodes; + +bool isClose(Real a, Real b) { + return Omega_h::are_close(a, b, 1e-12, 1e-12); +} + +// Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n. +// Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle, +// matching the LinearWedgeShape node ordering. +struct WedgeToVertexField { + static constexpr KOKKOS_FUNCTION Kokkos::Array + getTopology() { + return {MeshField::Wedge}; + } + + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap + operator()(MeshField::LO wedgeNodeIdx, MeshField::LO wedgeCompIdx, + MeshField::LO wedge, MeshField::Mesh_Topology topo) const { + assert(topo == MeshField::Wedge); + const MeshField::LO vtx = numNodes * wedge + wedgeNodeIdx; + return {0, wedgeCompIdx, vtx, MeshField::Vertex}; + } +}; + +// clang-format off +// wedge 0: right prism, triangle (0,0),(2,0),(0,3) swept by height 4 +// x(xi) = (2 xi0, 3 xi1, 2 (1 + zeta)), J = diag(2, 3, 2), det J = 12 +// wedge 1: the top triangle is shifted, scaled and tilted, so J varies with +// xi0, xi1 and zeta +const Vec3 wedgeCoords[numWedges][numNodes] = { + {{0, 0, 0}, {2, 0, 0}, {0, 3, 0}, + {0, 0, 4}, {2, 0, 4}, {0, 3, 4}}, + {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, + {0.2, 0.1, 1}, {1.7, 0.3, 1.2}, {0.1, 1.4, 0.9}}}; +const char *wedgeNames[numWedges] = {"right prism", "skewed wedge"}; + +// the six nodes, the centroid, and interior points +const Vec3 localPts[] = { + {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, + {1. / 3, 1. / 3, 0}, + {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; +// clang-format on +constexpr int numPts = static_cast(sizeof(localPts) / sizeof(localPts[0])); +constexpr int centroidPt = 6; + +// x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the +// bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates +Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + Vec3 x = Omega_h::zero_vector<3>(); + for (int a = 0; a < 3; ++a) { + x += L[a] * (1 - xi[2]) / 2 * nodes[a]; + x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; + } + return x; +} + +// rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating +// x(xi) directly: +// dx/dxi0 = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0) +// dx/dxi1 = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0) +// dx/dzeta = 1/2 sum_a L_a (T_a - B_a) +std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { + const Real down = (1 - xi[2]) / 2; + const Real up = (1 + xi[2]) / 2; + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + const Vec3 dXi0 = down * (nodes[1] - nodes[0]) + up * (nodes[4] - nodes[3]); + const Vec3 dXi1 = down * (nodes[2] - nodes[0]) + up * (nodes[5] - nodes[3]); + Vec3 dZeta = Omega_h::zero_vector<3>(); + for (int a = 0; a < 3; ++a) + dZeta += L[a] / 2 * (nodes[a + 3] - nodes[a]); + return {dXi0, dXi1, dZeta}; +} + +// Writes the wedge vertex coordinates into the coordinate field. +template +void setVtxCoords(CoordField field) { + const int numVtx = numWedges * numNodes; + Kokkos::View coords_h("coords_h", numVtx); + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < numNodes; ++n) + for (int d = 0; d < 3; ++d) + coords_h(numNodes * w + n, d) = wedgeCoords[w][n][d]; + auto coords = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), coords_h); + auto setCoordField = KOKKOS_LAMBDA(const int &vtx) { + for (int d = 0; d < 3; ++d) + field(vtx, 0, d, MeshField::Vertex) = coords(vtx, d); + }; + MeshField::parallel_for( + {0}, {numVtx}, setCoordField, "setCoordField"); +} + +// Compares the evaluated positions, Jacobians and determinants of wedge w with +// the hand-differentiated map at every local point. Returns the failure count. +template +int checkWedge(int w, const XView &x_h, const JView &J_h, + const DetView &detJ_h) { + const auto &nodes = wedgeCoords[w]; + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + // at a node, x is that node; at the centroid, x is the mean of the six + // vertices; everywhere else x follows the map + Vec3 xExpected = expectedX(nodes, localPts[p]); + if (p < numNodes) + xExpected = nodes[p]; + if (p == centroidPt) { + xExpected = Omega_h::zero_vector<3>(); + for (int n = 0; n < numNodes; ++n) + xExpected += nodes[n] / Real(numNodes); + } + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + + const auto jExpected = expectedJ(nodes, localPts[p]); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); + const Real detExpected = + jExpected[0] * Omega_h::cross(jExpected[1], jExpected[2]); + detFailed += !isClose(detJ_h(row), detExpected); + } + std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] + << ": evaluate failures " << evalFailed << ", Jacobian failures " + << jacFailed << ", det J failures " << detFailed << "\n"; + return evalFailed + jacFailed + detFailed; +} + +// Builds the coordinate field and the FieldElement, runs the evaluation, +// Jacobian and determinant, and checks every wedge. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeElementChecks() { + const MeshField::MeshInfo meshInfo{ + .numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; + auto coordField = MeshField::CreateCoordinateField< + ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); + setVtxCoords(coordField.field); + + MeshField::FieldElement fe(numWedges, coordField.field, + MeshField::LinearWedgeShape(), + WedgeToVertexField()); + + // the same local points in every wedge, stored element-major + Kokkos::View lc("localCoords", numWedges * numPts, 3); + auto lc_h = Kokkos::create_mirror_view(lc); + for (int w = 0; w < numWedges; ++w) + for (int p = 0; p < numPts; ++p) + for (int d = 0; d < 3; ++d) + lc_h(w * numPts + p, d) = localPts[p][d]; + Kokkos::deep_copy(lc, lc_h); + + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + const auto J_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + + int numFailed = 0; + for (int w = 0; w < numWedges; ++w) + numFailed += checkWedge(w, x_h, J_h, detJ_h); + return numFailed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runLinearWedgeElementChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} \ No newline at end of file diff --git a/test/testLinearWedgeShape.cpp b/test/testLinearWedgeShape.cpp new file mode 100644 index 0000000..267c966 --- /dev/null +++ b/test/testLinearWedgeShape.cpp @@ -0,0 +1,140 @@ +#include +#include +#include +#include + +// Checks of MeshField::LinearWedgeShape. The checks run inside Kokkos kernels +// so the shape functions are exercised in the default execution space, and +// failures are counted explicitly so the test also fails in builds where +// assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Shape = MeshField::LinearWedgeShape; + +KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { + return Kokkos::fabs(a - b) <= tol; +} + +// A function in the span of the wedge basis, +// {1, xi0, xi1, zeta, xi0*zeta, xi1*zeta}, and its gradient +KOKKOS_INLINE_FUNCTION Real f(Real x0, Real x1, Real z) { + return 0.5 + 2 * x0 - 3 * x1 + 0.25 * z + 4 * x0 * z - 1.5 * x1 * z; +} +KOKKOS_INLINE_FUNCTION MeshField::Vector3 gradF(Real x0, Real x1, Real z) { + return {2 + 4 * z, -3 - 1.5 * z, 0.25 + 4 * x0 - 1.5 * x1}; +} + +// Each shape function is one at its own node and zero at the others. +int checkKroneckerDelta() { + int numFailed = 0; + Kokkos::parallel_reduce( + "wedgeKronecker", Kokkos::RangePolicy(0, Shape::numNodes), + KOKKOS_LAMBDA(const int node, int &failed) { + const Shape shp; + const auto coords = shp.getNodeParametricCoords(); + const MeshField::Vector3 xi{coords[3 * node], coords[3 * node + 1], + coords[3 * node + 2]}; + const auto N = shp.getValues(xi); + for (size_t m = 0; m < Shape::numNodes; ++m) { + const Real expected = (static_cast(node) == m) ? 1.0 : 0.0; + failed += !isClose(N[m], expected, MeshField::MachinePrecision); + } + }, + numFailed); + return numFailed; +} + +// Interior points, spanning the zeta range. +Kokkos::View makeInteriorPoints() { + // clang-format off + const Real pts[][3] = {{1. / 3, 1. / 3, 0.0}, + {0.2, 0.3, -0.4}, + {0.6, 0.1, 0.7}, + {0.05, 0.9, 0.95}, + {0.45, 0.45, -0.95}}; + // clang-format on + const int numPts = static_cast(sizeof(pts) / sizeof(pts[0])); + Kokkos::View pts_d("wedgePts", numPts); + auto pts_h = Kokkos::create_mirror_view(pts_d); + for (int i = 0; i < numPts; ++i) + for (int d = 0; d < 3; ++d) + pts_h(i, d) = pts[i][d]; + Kokkos::deep_copy(pts_d, pts_h); + return pts_d; +} + +// At interior points: the values sum to one, the gradients sum to zero, the +// basis reproduces f and its gradient exactly, and the analytic gradients +// match central finite differences of the values. +int checkInteriorPoints() { + const auto pts_d = makeInteriorPoints(); + const int numPts = static_cast(pts_d.extent(0)); + + int numFailed = 0; + Kokkos::parallel_reduce( + "wedgeInterior", Kokkos::RangePolicy(0, numPts), + KOKKOS_LAMBDA(const int i, int &failed) { + const Shape shp; + const MeshField::Vector3 xi{pts_d(i, 0), pts_d(i, 1), pts_d(i, 2)}; + const auto N = shp.getValues(xi); + const auto G = shp.getLocalGradients(xi); + const auto nodes = shp.getNodeParametricCoords(); + + Real sumN = 0; + Real sumG[3] = {0, 0, 0}; + Real fInterp = 0; + Real gInterp[3] = {0, 0, 0}; + for (size_t n = 0; n < Shape::numNodes; ++n) { + const Real fn = f(nodes[3 * n], nodes[3 * n + 1], nodes[3 * n + 2]); + sumN += N[n]; + fInterp += N[n] * fn; + for (int d = 0; d < 3; ++d) { + sumG[d] += G[3 * n + d]; + gInterp[d] += G[3 * n + d] * fn; + } + } + const Real tol = 10 * MeshField::MachinePrecision; + failed += !isClose(sumN, 1.0, tol); + failed += !isClose(fInterp, f(xi[0], xi[1], xi[2]), tol); + const auto gExact = gradF(xi[0], xi[1], xi[2]); + for (int d = 0; d < 3; ++d) { + failed += !isClose(sumG[d], 0.0, tol); + failed += !isClose(gInterp[d], gExact[d], tol); + } + + // gradients against central finite differences of the values + const Real h = 1e-6; + for (int d = 0; d < 3; ++d) { + auto xp = xi; + auto xm = xi; + xp[d] += h; + xm[d] -= h; + const auto Np = shp.getValues(xp); + const auto Nm = shp.getValues(xm); + for (size_t n = 0; n < Shape::numNodes; ++n) + failed += !isClose((Np[n] - Nm[n]) / (2 * h), G[3 * n + d], 1e-8); + } + }, + numFailed); + return numFailed; +} + +// Runs both checks and prints their failure counts. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeShapeChecks() { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearWedgeShape interior point failures: " << interior + << "\n"; + return kronecker + interior; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runLinearWedgeShapeChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} \ No newline at end of file