diff --git a/CMakeLists.txt b/CMakeLists.txt index b385901..3e16c15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -220,6 +220,13 @@ 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(TriangleIntegrationTest test/testTriangleIntegration.cpp) +meshfields_add_exe(WedgeIntegrationTest test/testWedgeIntegration.cpp) +meshfields_add_exe(ExtrudedWedgeTest test/testExtrudedWedge.cpp) meshfields_add_exe(Interpolate test/testInterpolate.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -246,6 +253,13 @@ 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(TriangleIntegrationTest ./TriangleIntegrationTest) +test_func(WedgeIntegrationTest ./WedgeIntegrationTest) +test_func(ExtrudedWedgeTest ./ExtrudedWedgeTest) test_func(Interpolate ./Interpolate) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error diff --git a/src/MeshField.hpp b/src/MeshField.hpp index a967837..252acd1 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -284,6 +284,193 @@ template auto getTetrahedronElement(Omega_h::Mesh &mesh) { } } +/** + * @brief Maps the nodes of LinearWedgeShape wedges, extruded from a 2D + * triangle mesh, to vertices + * @details + * The 2D mesh is stacked into numLayers layers. Every plane has the vertex + * numbering and triangle connectivity of the 2D mesh; only the coordinates, + * which the caller sets, differ, so the triangles may change shape from plane + * to plane. Vertex v on plane p is vertex p * nVtx2d + v. Wedge w is in layer + * k = w / nTri2d and comes from triangle t = w % nTri2d. Layer k runs from + * plane k to plane k + 1, or to plane (k + 1) % numLayers when periodic. + * Wedge corner a is local vertex a of triangle t, or local vertex {0, 2, 1}[a] + * when reverseOrientation is true. Node a is corner a on plane k and node + * a + 3 is corner a on the top plane of the layer. These vertex ids, and the + * vertical edge ids of ExtrudedLinearCubicWedgeToField, number the extruded + * mesh; they are not entity ids of an Omega_h mesh. + * + * With coordinates (R, phi, Z), phi increasing with the layer and triangles + * counterclockwise in (R, Z), det J < 0. Setting reverseOrientation = true + * makes it positive. + */ +struct ExtrudedLinearWedgeToVertexField { + Omega_h::LOs triVerts; + int nVtx2d; + int nTri2d; + int numLayers; + bool periodic; + bool reverseOrientation; + ExtrudedLinearWedgeToVertexField(Omega_h::Mesh &mesh2d, int numLayersIn, + bool periodicIn, bool reverseOrientationIn) + : triVerts(mesh2d.ask_elem_verts()), nVtx2d(mesh2d.nverts()), + nTri2d(mesh2d.nfaces()), numLayers(numLayersIn), periodic(periodicIn), + reverseOrientation(reverseOrientationIn) { + if (mesh2d.dim() != 2 || mesh2d.family() != OMEGA_H_SIMPLEX) { + MeshField::fail( + "The mesh passed to %s must be 2D and simplex (triangles)\n", + __func__); + } + if (numLayers < 1 || (periodic && numLayers < 2)) { + MeshField::fail("%s needs numLayers >= 1, and >= 2 when periodic; got " + "%d\n", + __func__, numLayers); + } + } + + 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 layer = wedge / nTri2d; + const MeshField::LO tri = wedge % nTri2d; + const MeshField::LO corner = wedgeNodeIdx % 3; + // reverseOrientation swaps corners 1 and 2 + const MeshField::LO triNodeIdx = + (reverseOrientation && corner != 0) ? 3 - corner : corner; + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored); + const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); + const MeshField::LO vtx2d = triVerts[(tri * triToVtxDegree) + localVtxIdx]; + const MeshField::LO top = periodic ? (layer + 1) % numLayers : layer + 1; + const MeshField::LO plane = (wedgeNodeIdx < 3) ? layer : top; + return {0, wedgeCompIdx, plane * nVtx2d + vtx2d, MeshField::Vertex}; + } +}; + +/** + * @brief Maps the nodes of LinearCubicWedgeShape wedges, extruded from a 2D + * triangle mesh, to vertices and vertical edges + * @details + * Nodes 0-5 map as in ExtrudedLinearWedgeToVertexField. The vertical edge + * above 2D vertex v in layer k is edge k * nVtx2d + v. Nodes 6 + 2a and + * 7 + 2a are DOFs 0 (zeta = -1/3) and 1 (zeta = +1/3) of the edge above + * corner a. The orientation note of ExtrudedLinearWedgeToVertexField applies. + */ +struct ExtrudedLinearCubicWedgeToField { + Omega_h::LOs triVerts; + int nVtx2d; + int nTri2d; + int numLayers; + bool periodic; + bool reverseOrientation; + ExtrudedLinearCubicWedgeToField(Omega_h::Mesh &mesh2d, int numLayersIn, + bool periodicIn, bool reverseOrientationIn) + : triVerts(mesh2d.ask_elem_verts()), nVtx2d(mesh2d.nverts()), + nTri2d(mesh2d.nfaces()), numLayers(numLayersIn), periodic(periodicIn), + reverseOrientation(reverseOrientationIn) { + if (mesh2d.dim() != 2 || mesh2d.family() != OMEGA_H_SIMPLEX) { + MeshField::fail( + "The mesh passed to %s must be 2D and simplex (triangles)\n", + __func__); + } + if (numLayers < 1 || (periodic && numLayers < 2)) { + MeshField::fail("%s needs numLayers >= 1, and >= 2 when periodic; got " + "%d\n", + __func__, numLayers); + } + } + + 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 layer = wedge / nTri2d; + const MeshField::LO tri = wedge % nTri2d; + const bool isVtx = wedgeNodeIdx < 6; + const MeshField::LO corner = + isVtx ? wedgeNodeIdx % 3 : (wedgeNodeIdx - 6) / 2; + // reverseOrientation swaps corners 1 and 2 + const MeshField::LO triNodeIdx = + (reverseOrientation && corner != 0) ? 3 - corner : corner; + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored); + const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); + const MeshField::LO vtx2d = triVerts[(tri * triToVtxDegree) + localVtxIdx]; + if (isVtx) { + const MeshField::LO top = periodic ? (layer + 1) % numLayers : layer + 1; + const MeshField::LO plane = (wedgeNodeIdx < 3) ? layer : top; + return {0, wedgeCompIdx, plane * nVtx2d + vtx2d, MeshField::Vertex}; + } + return {(wedgeNodeIdx - 6) % 2, wedgeCompIdx, layer * nVtx2d + vtx2d, + MeshField::Edge}; + } +}; + +/** + * @brief MeshInfo of the wedges from getExtrudedWedgeElement + * @details numEdge counts only the vertical edges, which hold the edge DOFs + * of LinearCubicWedgeShape + */ +inline MeshField::MeshInfo getExtrudedWedgeMeshInfo(Omega_h::Mesh &mesh2d, + int numLayers, + bool periodic = false) { + const int numPlanes = periodic ? numLayers : numLayers + 1; + MeshField::MeshInfo meshInfo; + meshInfo.dim = 3; + meshInfo.numVtx = numPlanes * mesh2d.nverts(); + meshInfo.numEdge = numLayers * mesh2d.nverts(); + meshInfo.numWedge = numLayers * mesh2d.nfaces(); + return meshInfo; +} + +/** + * @brief Shape and element map for wedges extruded from a 2D triangle mesh + * @tparam ShapeOrder 1 for LinearWedgeShape, 3 for LinearCubicWedgeShape + * @details See ExtrudedLinearWedgeToVertexField for the numbering and + * orientation. The caller fills the coordinate field, sized by + * getExtrudedWedgeMeshInfo. + */ +template +auto getExtrudedWedgeElement(Omega_h::Mesh &mesh2d, int numLayers, + bool periodic = false, + bool reverseOrientation = false) { + static_assert(ShapeOrder == 1 || ShapeOrder == 3); + if constexpr (ShapeOrder == 1) { + struct result { + MeshField::LinearWedgeShape shp; + ExtrudedLinearWedgeToVertexField map; + }; + return result{MeshField::LinearWedgeShape(), + ExtrudedLinearWedgeToVertexField(mesh2d, numLayers, periodic, + reverseOrientation)}; + } else if constexpr (ShapeOrder == 3) { + struct result { + MeshField::LinearCubicWedgeShape shp; + ExtrudedLinearCubicWedgeToField map; + }; + return result{MeshField::LinearCubicWedgeShape(), + ExtrudedLinearCubicWedgeToField(mesh2d, numLayers, periodic, + reverseOrientation)}; + } +} + } // namespace Omegah /** 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_Integrate.hpp b/src/MeshField_Integrate.hpp index 6587644..cfc4db2 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -48,6 +48,102 @@ template class EntityIntegration { virtual Integration const *getIntegration(int i) const = 0; }; +// Gauss-Legendre rules on [-1, 1]. N1 to N6 are copied from SCOREC/core +// @ 7cd76473 apf/apfIntegrate.cc (EdgeIntegration, lines 27-129). The n-point +// rule is exact for degree 2n - 1, and the weights sum to 2, the length of +// the reference edge. +//! [EdgeIntegration] +class EdgeIntegration : public EntityIntegration<1> { +public: + class N1 : public Integration<1> { + public: + virtual int countPoints() const { return 1; } + virtual std::vector> getPoints() const { + return {IntegrationPoint(Vector1{0}, 2, 1, 0)}; + } + virtual int getAccuracy() const { return 1; } + }; // end N1 + class N2 : public Integration<1> { + public: + virtual int countPoints() const { return 2; } + virtual std::vector> getPoints() const { + return {IntegrationPoint(Vector1{0.577350269189626}, 1, 1, 0), + IntegrationPoint(Vector1{-0.577350269189626}, 1, 1, 0)}; + } + virtual int getAccuracy() const { return 3; } + }; // end N2 + class N3 : public Integration<1> { + public: + virtual int countPoints() const { return 3; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.000000000000000}, 0.888888888888889, 1, 0), + IntegrationPoint(Vector1{0.774596669241483}, 0.555555555555556, 1, 0), + IntegrationPoint(Vector1{-0.774596669241483}, 0.555555555555556, 1, + 0)}; + } + virtual int getAccuracy() const { return 5; } + }; // end N3 + class N4 : public Integration<1> { + public: + virtual int countPoints() const { return 4; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.339981043584856}, 0.652145154862546, 1, 0), + IntegrationPoint(Vector1{-0.339981043584856}, 0.652145154862546, 1, + 0), + IntegrationPoint(Vector1{0.861136311594053}, 0.347854845137454, 1, 0), + IntegrationPoint(Vector1{-0.861136311594053}, 0.347854845137454, 1, + 0)}; + } + virtual int getAccuracy() const { return 7; } + }; // end N4 + class N5 : public Integration<1> { + public: + virtual int countPoints() const { return 5; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.000000000000000}, 0.568888888888889, 1, 0), + IntegrationPoint(Vector1{0.538469310105683}, 0.478628670499366, 1, 0), + IntegrationPoint(Vector1{-0.538469310105683}, 0.478628670499366, 1, + 0), + IntegrationPoint(Vector1{0.906179845938664}, 0.236926885056189, 1, 0), + IntegrationPoint(Vector1{-0.906179845938664}, 0.236926885056189, 1, + 0)}; + } + virtual int getAccuracy() const { return 9; } + }; // end N5 + class N6 : public Integration<1> { + public: + virtual int countPoints() const { return 6; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.238619186083197}, 0.467913934572691, 1, 0), + IntegrationPoint(Vector1{-0.238619186083197}, 0.467913934572691, 1, + 0), + IntegrationPoint(Vector1{0.661209386466265}, 0.360761573048139, 1, 0), + IntegrationPoint(Vector1{-0.661209386466265}, 0.360761573048139, 1, + 0), + IntegrationPoint(Vector1{0.932469514203152}, 0.171324492379170, 1, 0), + IntegrationPoint(Vector1{-0.932469514203152}, 0.171324492379170, 1, + 0)}; + } + virtual int getAccuracy() const { return 11; } + }; // end N6 + virtual int countIntegrations() const { return 6; } + virtual Integration<1> const *getIntegration(int i) const { + static N1 i1; + static N2 i2; + static N3 i3; + static N4 i4; + static N5 i5; + static N6 i6; + static Integration<1> *integrations[6] = {&i1, &i2, &i3, &i4, &i5, &i6}; + return integrations[i]; + } +}; +//! [EdgeIntegration] + // points are given in the reduced parametric coordinates (xi0=L1, xi1=L2, // ...); the redundant barycentric coordinate L0 = 1-sum(xi) is omitted //! [TriangleIntegration] @@ -74,11 +170,40 @@ class TriangleIntegration : public EntityIntegration<2> { } virtual int getAccuracy() const { return 2; } }; // end N2 - virtual int countIntegrations() const { return 2; } + + // Symmetric Gaussian quadrature rules for the triangle from + // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature + // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). + // N4 is the p = 4 rule of Appendix II (p. 1140), copied from + // SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). The weights + // sum to one, so they are divided by 2, the area of the reference triangle. + // The p = 3 rule has a negative weight, so order 3 uses N4. + class N4 : public Integration<2> { + public: + virtual int countPoints() const { return 6; } + virtual std::vector> getPoints() const { + return {IntegrationPoint(Vector2{0.816847572980459, 0.091576213509771}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.091576213509771, 0.816847572980459}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.091576213509771, 0.091576213509771}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.108103018168070, 0.445948490915965}, + 0.223381589678011 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.445948490915965, 0.108103018168070}, + 0.223381589678011 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.445948490915965, 0.445948490915965}, + 0.223381589678011 / 2.0, 2, 0)}; + } + virtual int getAccuracy() const { return 4; } + }; // end N4 + + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; static N2 i2; - static Integration<2> *integrations[2] = {&i1, &i2}; + static N4 i4; + static Integration<2> *integrations[3] = {&i1, &i2, &i4}; return integrations[i]; } }; @@ -122,9 +247,36 @@ class TetrahedronIntegration : public EntityIntegration<3> { return integrations[i]; } }; + +// Wedge rules are the tensor product of a triangle rule in (xi0, xi1) and an +// edge rule in zeta. The rule for (triOrder, zetaOrder) integrates +// xi0^a xi1^b zeta^c exactly for a + b <= triOrder and c <= zetaOrder. The +// weights sum to 1, the volume of the reference wedge (area 1/2 times length +// 2). +inline std::vector> +getWedgeIntegrationPoints(int triOrder, int zetaOrder) { + const auto *tri = TriangleIntegration().getAccurate(triOrder); + if (tri == nullptr) { + fail("getWedgeIntegrationPoints: no triangle rule of order %d\n", triOrder); + } + const auto *edge = EdgeIntegration().getAccurate(zetaOrder); + if (edge == nullptr) { + fail("getWedgeIntegrationPoints: no edge rule of order %d\n", zetaOrder); + } + std::vector> points; + for (auto const &t : tri->getPoints()) + for (auto const &z : edge->getPoints()) + points.push_back( + IntegrationPoint(Vector3{t.param[0], t.param[1], z.param[0]}, + t.weight * z.weight, 3, 0)); + return points; +} + //! [getIntegration] template auto const getIntegration() { - if constexpr (topo == Triangle) { + if constexpr (topo == Edge) { + return std::make_shared(); + } else if constexpr (topo == Triangle) { return std::make_shared(); } else if constexpr (topo == Tetrahedron) { return std::make_shared(); @@ -132,9 +284,16 @@ template auto const getIntegration() { fail("getIntegration does not support given topology\n"); } //! [getIntegration] -template auto getIntegrationPoints(int order) { - auto ip = getIntegration()->getAccurate(order)->getPoints(); - return ip; + +// zetaOrder is used only for wedges and defaults to order +template +auto getIntegrationPoints(int order, int zetaOrder = -1) { + if constexpr (topo == Wedge) { + return getWedgeIntegrationPoints(order, zetaOrder < 0 ? order : zetaOrder); + } else { + auto ip = getIntegration()->getAccurate(order)->getPoints(); + return ip; + } } template @@ -208,7 +367,12 @@ auto getJacobianDeterminants(FieldElement &fes, class Integrator { public: /** \brief Construct an Integrator given an order of accuracy. */ - Integrator(int o) : order(o) {} + Integrator(int o) : order(o), zetaOrder(o) {} + /** \brief Construct an Integrator for wedges, with separate orders of + * accuracy for the triangle (xi0, xi1) and zeta directions. */ + Integrator(int triOrderIn, int zetaOrderIn) + : order(triOrderIn), zetaOrder(zetaOrderIn) {} + virtual ~Integrator(){}; /** \brief User callback: process entry * @@ -248,7 +412,7 @@ class Integrator { template void process(FieldElement &fes) { constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); pre(); - auto ip = getIntegrationPoints(order); + auto ip = getIntegrationPoints(order, zetaOrder); auto localCoords = getIntegrationPointLocalCoords(fes, ip); auto weights = getIntegrationPointWeights(fes, ip); auto dV = getJacobianDeterminants(fes, localCoords); @@ -265,6 +429,7 @@ class Integrator { protected: int order; + int zetaOrder; }; } // namespace MeshField 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/testExtrudedWedge.cpp b/test/testExtrudedWedge.cpp new file mode 100644 index 0000000..d577415 --- /dev/null +++ b/test/testExtrudedWedge.cpp @@ -0,0 +1,413 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_For.hpp" +#include "MeshField_Integrate.hpp" +#include "MeshField_ShapeField.hpp" +#include "Omega_h_build.hpp" +#include "Omega_h_library.hpp" +#include +#include +#include +#include +#include + +// Checks the extruded wedge element maps, getExtrudedWedgeElement and +// getExtrudedWedgeMeshInfo, on a 2D triangle mesh stacked into numLayers +// layers: +// 1. the DOF holder of every wedge node, for open and periodic stacks, linear +// and cubic, with and without reverseOrientation, and the MeshInfo counts, +// 2. the field positions at the nodes and interior points of every wedge, +// 3. the volume from the Integrator, +// 4. the sign of det J with (R, phi, Z) = (x, s, y). +// The coordinates come from analytic maps X(x, y, s) that the elements +// reproduce exactly. Checks 2-4 use open stacks: in a periodic stack the last +// layer closes onto plane 0, which a single map X cannot describe. Failures +// are counted explicitly so the test also fails in builds where assert() is +// disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::LO; +using MeshField::Real; +using Vec3 = std::array; + +namespace { + +constexpr int numLayers = 3; +constexpr Real H = 1.5; // height of the stack +constexpr Real h = H / numLayers; // height of one layer +constexpr Real A = 1.0; // area of the 2D mesh + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +// Analytic maps X(x, y, s) of the point (x, y) of the 2D mesh at height s +using Map = Vec3 (*)(Real, Real, Real); +Vec3 straightMap(Real x, Real y, Real s) { return {x, y, s}; } +Vec3 taperedMap(Real x, Real y, Real s) { return {x * (1 + 0.1 * s), y, s}; } +Vec3 curvedMap(Real x, Real y, Real s) { + return {x * (1 + 0.1 * s) + 0.02 * s * s * s, y, s + 0.05 * s * s * s}; +} +// R = x, phi = s, Z = y +Vec3 cylindricalMap(Real x, Real y, Real s) { return {x, s, y}; } + +// Host copies of the 2D mesh +struct HostMesh2d { + HostMesh2d(Omega_h::Mesh &mesh2d) + : triVerts(mesh2d.ask_elem_verts()), coords(mesh2d.coords()), + nVtx2d(mesh2d.nverts()), nTri2d(mesh2d.nfaces()) {} + Omega_h::HostRead triVerts; + Omega_h::HostRead coords; + int nVtx2d; + int nTri2d; + + // the 2D vertex of wedge corner a above triangle t + LO cornerVtx(LO t, int a, bool reverse) const { + const int corner[2][3] = {{0, 1, 2}, {0, 2, 1}}; + return triVerts[t * 3 + corner[reverse][a]]; + } + + Real signedArea(LO t) const { + Real x[3], y[3]; + for (int a = 0; a < 3; ++a) { + const LO v = triVerts[t * 3 + a]; + x[a] = coords[2 * v]; + y[a] = coords[2 * v + 1]; + } + return 0.5 * + ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0])); + } + + // X at the local point xi of wedge w in an open stack + Vec3 position(Map X, LO w, Vec3 const &xi, bool reverse) const { + const LO k = w / nTri2d, t = w % nTri2d; + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + Real x = 0, y = 0; + for (int a = 0; a < 3; ++a) { + const LO v = cornerVtx(t, a, reverse); + x += L[a] * coords[2 * v]; + y += L[a] * coords[2 * v + 1]; + } + return X(x, y, k * h + h * (1 + xi[2]) / 2); + } +}; + +template +auto createCoordField(const MeshField::MeshInfo &info) { + if constexpr (ShapeOrder == 1) { + return MeshField::CreateCoordinateField( + info); + } else { + return MeshField::CreateLagrangeField< + ExecutionSpace, MeshField::KokkosController, Real, 3, 3, 3, + MeshField::LinearCubicWedgeShape>(info); + } +} + +// Writes X into the coordinate field of an open stack: vertex (p, v) at +// X(x_v, y_v, s_p) and, for the cubic field, DOFs 0 and 1 of edge (k, v) at +// X(x_v, y_v, s_k + h/3) and X(x_v, y_v, s_k + 2h/3). +template +void setExtrudedCoords(CoordField field, const HostMesh2d &m, Map X) { + const int numVtx = (numLayers + 1) * m.nVtx2d; + const int numEdges = numLayers * m.nVtx2d; + Kokkos::View vtx_h("vtx_h", numVtx); + Kokkos::View edge_h("edge_h", numEdges); + for (int v = 0; v < m.nVtx2d; ++v) { + const Real x = m.coords[2 * v], y = m.coords[2 * v + 1]; + for (int p = 0; p <= numLayers; ++p) { + const Vec3 pos = X(x, y, p * h); + for (int d = 0; d < 3; ++d) + vtx_h(p * m.nVtx2d + v, d) = pos[d]; + } + for (int k = 0; k < numLayers; ++k) + for (int dof = 0; dof < 2; ++dof) { + const Vec3 pos = X(x, y, k * h + (dof + 1) * h / 3); + for (int d = 0; d < 3; ++d) + edge_h(k * m.nVtx2d + v, dof, d) = pos[d]; + } + } + auto vtx = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), vtx_h); + auto setVtxCoords = KOKKOS_LAMBDA(const int &i) { + for (int d = 0; d < 3; ++d) + field(i, 0, d, MeshField::Vertex) = vtx(i, d); + }; + MeshField::parallel_for( + {0}, {numVtx}, setVtxCoords, "setVtxCoords"); + if constexpr (ShapeOrder == 3) { + auto edge = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), edge_h); + auto setEdgeCoords = 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); + }; + MeshField::parallel_for( + {0}, {numEdges}, setEdgeCoords, "setEdgeCoords"); + } +} + +const std::vector interiorPts = { + {1. / 3, 1. / 3, 0}, {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; + +// The nodes of the shape, then interior points +template std::vector nodesAndInteriorPoints() { + const auto xi = Shape().getNodeParametricCoords(); + std::vector pts; + for (size_t n = 0; n < Shape::numNodes; ++n) + pts.push_back({xi[3 * n], xi[3 * n + 1], xi[3 * n + 2]}); + pts.insert(pts.end(), interiorPts.begin(), interiorPts.end()); + return pts; +} + +// the same local points in every wedge, stored element-major +Kokkos::View makeLocalCoords(int numWedges, + const std::vector &pts) { + const int numPts = static_cast(pts.size()); + 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); + return lc; +} + +const char *shapeName(int shapeOrder) { + return shapeOrder == 1 ? "linear" : "cubic"; +} + +// Gets the DOF holder of every (wedge, node) from the map on the device and +// compares it with the layout computed on the host. Also checks that the ids +// are in range, that the top of layer k is the bottom of layer k + 1 (and the +// top of the last layer is the bottom of layer 0 when periodic), and the +// MeshInfo counts. +template +int checkIds(Omega_h::Mesh &mesh2d, const HostMesh2d &m, bool periodic, + bool reverse) { + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, periodic, reverse); + const auto map = elem.map; + const int numNodes = decltype(elem.shp)::numNodes; + const int numWedges = numLayers * m.nTri2d; + // node, entity, topology, and 1 if the component is passed through + Kokkos::View holders("holders", numWedges * numNodes, 4); + Kokkos::parallel_for( + "getHolders", + Kokkos::RangePolicy(0, numWedges * numNodes), + KOKKOS_LAMBDA(const int i) { + const LO w = i / numNodes, n = i % numNodes; + const auto holder = map(n, 0, w, MeshField::Wedge); + LO compOk = holder.component == 0; + for (LO c = 1; c < 3; ++c) + compOk = compOk && map(n, c, w, MeshField::Wedge).component == c; + holders(i, 0) = holder.node; + holders(i, 1) = holder.entity; + holders(i, 2) = holder.topo; + holders(i, 3) = compOk; + }); + const auto holders_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), holders); + + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers, periodic); + const int numPlanes = periodic ? numLayers : numLayers + 1; + int idFailed = 0, rangeFailed = 0, stackFailed = 0, countFailed = 0; + for (LO w = 0; w < numWedges; ++w) { + const LO k = w / m.nTri2d, t = w % m.nTri2d; + const LO top = periodic ? (k + 1) % numLayers : k + 1; + for (int n = 0; n < numNodes; ++n) { + const int i = w * numNodes + n; + const LO dof = holders_h(i, 0), ent = holders_h(i, 1), + topo = holders_h(i, 2); + LO expDof = 0, expEnt = 0, expTopo = MeshField::Vertex; + if (n < 6) { + const LO plane = n < 3 ? k : top; + expEnt = plane * m.nVtx2d + m.cornerVtx(t, n % 3, reverse); + } else { + expDof = (n - 6) % 2; + expEnt = k * m.nVtx2d + m.cornerVtx(t, (n - 6) / 2, reverse); + expTopo = MeshField::Edge; + } + idFailed += !(dof == expDof && ent == expEnt && topo == expTopo && + holders_h(i, 3) == 1); + if (topo == MeshField::Vertex) + rangeFailed += !(dof == 0 && ent >= 0 && ent < info.numVtx); + else + rangeFailed += !(topo == MeshField::Edge && (dof == 0 || dof == 1) && + ent >= 0 && ent < info.numEdge); + } + } + for (LO k = 0; k < numLayers; ++k) { + if (k + 1 == numLayers && !periodic) + continue; + const LO next = (k + 1) % numLayers; + for (LO t = 0; t < m.nTri2d; ++t) + for (int a = 0; a < 3; ++a) + stackFailed += holders_h((k * m.nTri2d + t) * numNodes + a + 3, 1) != + holders_h((next * m.nTri2d + t) * numNodes + a, 1); + } + countFailed += info.dim != 3; + countFailed += info.numVtx != numPlanes * m.nVtx2d; + countFailed += info.numEdge != numLayers * m.nVtx2d; + countFailed += info.numWedge != numWedges; + countFailed += (info.numTri || info.numQuad || info.numTet || info.numHex || + info.numPrism || info.numPyramid); + std::cout << "ids, " << shapeName(ShapeOrder) << ", " + << (periodic ? "periodic" : "open") << ", reverseOrientation " + << reverse << ": id failures " << idFailed << ", range failures " + << rangeFailed << ", stacking failures " << stackFailed + << ", MeshInfo failures " << countFailed << "\n"; + return idFailed + rangeFailed + stackFailed + countFailed; +} + +// Evaluates the field at the nodes and interior points of every wedge of an +// open stack and compares with X. +template +int checkPositions(Omega_h::Mesh &mesh2d, const HostMesh2d &m, Map X, + bool reverse, const char *mapName) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, X); + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, false, reverse); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + const auto pts = nodesAndInteriorPoints(); + const int numPts = static_cast(pts.size()); + const auto lc = makeLocalCoords(info.numWedge, pts); + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + int failed = 0; + for (int w = 0; w < info.numWedge; ++w) + for (int p = 0; p < numPts; ++p) { + const Vec3 expected = m.position(X, w, pts[p], reverse); + for (int d = 0; d < 3; ++d) + failed += !isClose(x_h(w * numPts + p, d), expected[d], 1e-12); + } + std::cout << "positions, " << shapeName(ShapeOrder) << ", " << mapName + << ", reverseOrientation " << reverse << ": failures " << failed + << "\n"; + return failed; +} + +// sums w * dV over every point of every wedge +class VolumeIntegrator : public MeshField::Integrator { +public: + VolumeIntegrator(int triOrder, int zetaOrder) + : MeshField::Integrator(triOrder, zetaOrder) {} + void atPoints(Kokkos::View, Kokkos::View w, + Kokkos::View dV) override { + // w holds one weight per point of a single element; dV holds one entry + // per point of every element + const int numPtsPerElem = w.extent(0); + Real sum = 0; + Kokkos::parallel_reduce( + "volume", dV.extent(0), + KOKKOS_LAMBDA(const int i, Real &s) { + s += w(i % numPtsPerElem) * dV(i); + }, + sum); + volume = sum; + } + Real volume = 0; +}; + +// dV is signed, so X must give positively oriented wedges +template +int checkVolume(Omega_h::Mesh &mesh2d, const HostMesh2d &m, Map X, + Real expected, const char *mapName) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, X); + const auto elem = + MeshField::Omegah::getExtrudedWedgeElement(mesh2d, numLayers); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + VolumeIntegrator vi(3, 11); + vi.process(fe); + const bool ok = isClose(vi.volume, expected, 1e-12); + std::cout << "volume, " << shapeName(ShapeOrder) << ", " << mapName << ": " + << vi.volume << ", expected " << expected << (ok ? "" : " FAIL") + << "\n"; + return !ok; +} + +// With (R, phi, Z) = (x, s, y), det J has the sign -sign(area) of the 2D +// triangle, and +sign(area) with reverseOrientation. Uses getJacobians, not +// the Integrator, since Debug builds assert dV >= 0 there. +template +int checkOrientation(Omega_h::Mesh &mesh2d, const HostMesh2d &m, bool reverse) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, cylindricalMap); + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, false, reverse); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + const int numPts = static_cast(interiorPts.size()); + const auto lc = makeLocalCoords(info.numWedge, interiorPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + int failed = 0; + for (int w = 0; w < info.numWedge; ++w) { + const Real area = m.signedArea(w % m.nTri2d); + const Real areaSign = (area > 0) - (area < 0); + const Real expectedSign = reverse ? areaSign : -areaSign; + for (int p = 0; p < numPts; ++p) + failed += !(detJ_h(w * numPts + p) * expectedSign > 0); + } + std::cout << "orientation, " << shapeName(ShapeOrder) + << ", reverseOrientation " << reverse << ": wrong signs " << failed + << " of " << info.numWedge * numPts << "\n"; + return failed; +} + +int runExtrudedWedgeChecks(Omega_h::Mesh &mesh2d) { + const HostMesh2d m(mesh2d); + int failed = 0; + for (bool periodic : {false, true}) + for (bool reverse : {false, true}) { + failed += checkIds<1>(mesh2d, m, periodic, reverse); + failed += checkIds<3>(mesh2d, m, periodic, reverse); + } + for (bool reverse : {false, true}) { + failed += checkPositions<1>(mesh2d, m, taperedMap, reverse, "tapered"); + failed += checkPositions<3>(mesh2d, m, curvedMap, reverse, "curved"); + } + std::cout.precision(16); + failed += checkVolume<1>(mesh2d, m, straightMap, A * H, "straight"); + failed += + checkVolume<1>(mesh2d, m, taperedMap, A * (H + 0.05 * H * H), "tapered"); + failed += checkVolume<3>( + mesh2d, m, curvedMap, + A * (H + 0.05 * H * H + 0.05 * H * H * H + 0.00375 * H * H * H * H), + "sheared and curved"); + for (bool reverse : {false, true}) { + failed += checkOrientation<1>(mesh2d, m, reverse); + failed += checkOrientation<3>(mesh2d, m, reverse); + } + return failed; +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int numFailed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + auto mesh2d = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, + 0.0, 2, 2, 0); + numFailed = runExtrudedWedgeChecks(mesh2d); + } + Kokkos::finalize(); + return numFailed ? 1 : 0; +} 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 diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp new file mode 100644 index 0000000..8aacc8d --- /dev/null +++ b/test/testTriangleIntegration.cpp @@ -0,0 +1,139 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Integrate.hpp" +#include "Omega_h_build.hpp" +#include +#include +#include + +// Checks MeshField::TriangleIntegration: +// 1. every rule has positive weights and points inside the reference +// triangle, and integrates xi0^a xi1^b exactly for a + b up to its +// accuracy, +// 2. getAccurate(order) picks the expected rule for orders 1 to 4 and none +// for order 5, +// 3. Integrator::process with orders 3 and 4 integrates monomials of that +// degree exactly over a triangulated unit square. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +Real factorial(int n) { return n <= 1 ? 1.0 : n * factorial(n - 1); } + +// exact integral of xi0^a xi1^b over the reference triangle +Real triangleMonomial(int a, int b) { + return factorial(a) * factorial(b) / factorial(a + b + 2); +} + +// A rule of accuracy p must integrate xi0^a xi1^b exactly for a + b <= p. +// The a = b = 0 case checks that the weights sum to 1/2, the area of the +// reference triangle. +int checkRule(MeshField::Integration<2> const &rule) { + const auto pts = rule.getPoints(); + const int p = rule.getAccuracy(); + int failed = (static_cast(pts.size()) != rule.countPoints()); + for (auto const &pt : pts) { + const Real x0 = pt.param[0], x1 = pt.param[1]; + failed += !(pt.weight > 0 && x0 >= 0 && x1 >= 0 && x0 + x1 <= 1); + } + for (int a = 0; a <= p; ++a) { + for (int b = 0; a + b <= p; ++b) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], a) * std::pow(pt.param[1], b); + failed += !isClose(sum, triangleMonomial(a, b), 1e-14); + } + } + std::cout << "triangle rule, accuracy " << p << ", " << pts.size() + << " points: failures " << failed << "\n"; + return failed; +} + +// getAccurate(order) returns the first rule with accuracy >= order +int checkRuleSelection() { + const auto tri = MeshField::getIntegration(); + const int numPts[] = {1, 3, 6, 6}; // orders 1 to 4 + int failed = 0; + for (int order = 1; order <= 4; ++order) + failed += tri->getAccurate(order)->countPoints() != numPts[order - 1]; + failed += tri->getAccurate(5) != nullptr; + std::cout << "rule selection: failures " << failed << "\n"; + return failed; +} + +// Integrates x^px y^py over the mesh. The coordinate map is affine on each +// triangle, so the integrand has degree px + py in the parametric +// coordinates and a rule of that order is exact. +template +class MonomialIntegrator : public MeshField::Integrator { +public: + MonomialIntegrator(FieldElement &feIn, int order, int pxIn, int pyIn) + : MeshField::Integrator(order), fe(feIn), px(pxIn), py(pyIn) {} + void atPoints(Kokkos::View p, Kokkos::View w, + Kokkos::View dV) override { + // p and w hold one entry per point of a single element; x and dV hold + // one entry per point of every element + const auto x = MeshField::evaluateStaticLength(fe, p); + const int numPtsPerElem = w.extent(0); + const Real a = px, b = py; + Real sum = 0; + Kokkos::parallel_reduce( + "integrateMonomial", dV.extent(0), + KOKKOS_LAMBDA(const int i, Real &lsum) { + lsum += w(i % numPtsPerElem) * dV(i) * Kokkos::pow(x(i, 0), a) * + Kokkos::pow(x(i, 1), b); + }, + sum); + integral += sum; + } + FieldElement &fe; + int px, py; + Real integral = 0; +}; + +// integral of x^px y^py over the unit square +Real squareMonomial(int px, int py) { return 1.0 / ((px + 1) * (py + 1)); } + +int checkIntegrator(Omega_h::Mesh &mesh) { + MeshField::OmegahMeshField omf(mesh); + auto coordField = omf.getCoordField(); + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement fe(mesh.nelems(), coordField.field, shp, map); + // {order, px, py}: each monomial has degree equal to the order + const int cases[][3] = {{3, 2, 1}, {3, 0, 3}, {4, 2, 2}, {4, 3, 1}}; + int failed = 0; + std::cout.precision(16); + for (auto const &c : cases) { + MonomialIntegrator integrator(fe, c[0], c[1], c[2]); + integrator.process(fe); + const Real expected = squareMonomial(c[1], c[2]); + const bool ok = isClose(integrator.integral, expected, 1e-13); + std::cout << "order " << c[0] << ": integral of x^" << c[1] << " y^" << c[2] + << " = " << integrator.integral << ", expected " << expected + << (ok ? "" : " FAIL") << "\n"; + failed += !ok; + } + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int failed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + failed += checkRule(MeshField::TriangleIntegration::N1()); + failed += checkRule(MeshField::TriangleIntegration::N2()); + failed += checkRule(MeshField::TriangleIntegration::N4()); + failed += checkRuleSelection(); + auto tris = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, + 3, 3, 0); + failed += checkIntegrator(tris); + } + Kokkos::finalize(); + return failed ? 1 : 0; +} diff --git a/test/testWedgeIntegration.cpp b/test/testWedgeIntegration.cpp new file mode 100644 index 0000000..063175d --- /dev/null +++ b/test/testWedgeIntegration.cpp @@ -0,0 +1,202 @@ +#include "KokkosController.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_For.hpp" +#include "MeshField_Integrate.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include +#include + +// Checks the wedge quadrature: +// 1. every EdgeIntegration rule has positive weights and points inside +// (-1, 1), and integrates zeta^c exactly for c up to its accuracy, +// 2. getWedgeIntegrationPoints(triOrder, zetaOrder) has the expected number +// of points, positive weights and points inside the reference wedge, and +// integrates xi0^a xi1^b zeta^c exactly for a + b <= triOrder and +// c <= zetaOrder, +// 3. Integrator::process with LinearWedgeShape gives the exact volume of a +// right prism and a tapered wedge, and a zeta order that is too low gives +// the wrong volume for the tapered wedge. +// Failures are counted explicitly so the test also fails in builds where +// assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +Real factorial(int n) { return n <= 1 ? 1.0 : n * factorial(n - 1); } + +// exact integral of xi0^a xi1^b over the reference triangle +Real triangleMonomial(int a, int b) { + return factorial(a) * factorial(b) / factorial(a + b + 2); +} + +// exact integral of zeta^c over [-1, 1] +Real lineMonomial(int c) { return (c % 2 == 1) ? 0.0 : 2.0 / (c + 1); } + +// A rule of accuracy p must integrate zeta^c exactly for c <= p. The c = 0 +// case checks that the weights sum to 2, the length of the reference edge. +int checkEdgeRule(MeshField::Integration<1> const &rule) { + const auto pts = rule.getPoints(); + const int p = rule.getAccuracy(); + int failed = (static_cast(pts.size()) != rule.countPoints()); + for (auto const &pt : pts) + failed += !(pt.weight > 0 && pt.param[0] > -1 && pt.param[0] < 1); + for (int c = 0; c <= p; ++c) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], c); + failed += !isClose(sum, lineMonomial(c), 1e-14); + } + std::cout << "edge rule, accuracy " << p << ", " << pts.size() + << " points: failures " << failed << "\n"; + return failed; +} + +int checkWedgeRule(int triOrder, int zetaOrder, int numTriPts, + int numGaussPts) { + const auto pts = MeshField::getWedgeIntegrationPoints(triOrder, zetaOrder); + int failed = (static_cast(pts.size()) != numTriPts * numGaussPts); + for (auto const &pt : pts) { + const Real x0 = pt.param[0], x1 = pt.param[1], z = pt.param[2]; + failed += !(pt.weight > 0 && x0 >= 0 && x1 >= 0 && x0 + x1 <= 1 && z > -1 && + z < 1); + } + // a = b = c = 0 checks that the weights sum to 1 + for (int a = 0; a <= triOrder; ++a) + for (int b = 0; a + b <= triOrder; ++b) + for (int c = 0; c <= zetaOrder; ++c) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], a) * + std::pow(pt.param[1], b) * std::pow(pt.param[2], c); + failed += + !isClose(sum, triangleMonomial(a, b) * lineMonomial(c), 1e-14); + } + std::cout << "wedge rule (" << triOrder << ", " << zetaOrder << "), " + << pts.size() << " points: failures " << failed << "\n"; + return failed; +} + +// node n of wedge w is vertex 6*w + n +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); + return {0, wedgeCompIdx, 6 * wedge + wedgeNodeIdx, MeshField::Vertex}; + } +}; + +// sums w * |det J| over every point of every wedge +class VolumeIntegrator : public MeshField::Integrator { +public: + VolumeIntegrator(int triOrder, int zetaOrder) + : MeshField::Integrator(triOrder, zetaOrder) {} + void atPoints(Kokkos::View, Kokkos::View w, + Kokkos::View dV) override { + // w holds one weight per point of a single element; dV holds one entry + // per point of every element + const int numPtsPerElem = w.extent(0); + Real sum = 0; + Kokkos::parallel_reduce( + "volume", dV.extent(0), + KOKKOS_LAMBDA(const int i, Real &s) { + s += w(i % numPtsPerElem) * dV(i); + }, + sum); + volume = sum; + } + Real volume = 0; +}; + +// clang-format off +// wedge 0: right prism, triangle (0,0),(2,0),(0,3) times height 4, volume 12 +// wedge 1: tapered wedge, bottom triangle (0,0),(2,0),(0,2) at z = 0 and top +// triangle (0,0),(1,0),(0,1) at z = 1, a frustum with volume +// h/3 (A0 + A1 + sqrt(A0 A1)) = (2 + 1/2 + 1)/3 = 7/6 +constexpr int numWedges = 2; +const Real wedgeCoords[numWedges][6][3] = { + {{0, 0, 0}, {2, 0, 0}, {0, 3, 0}, {0, 0, 4}, {2, 0, 4}, {0, 3, 4}}, + {{0, 0, 0}, {2, 0, 0}, {0, 2, 0}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}}}; +// clang-format on +const Real exactVolume = 12.0 + 7.0 / 6.0; + +template void setVtxCoords(CoordField field) { + const int numVtx = 6 * numWedges; + Kokkos::View coords_h("coords_h", numVtx); + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < 6; ++n) + for (int d = 0; d < 3; ++d) + coords_h(6 * 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"); +} + +int checkIntegrator() { + const MeshField::MeshInfo meshInfo{ + .numVtx = 6 * numWedges, .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()); + // det J of a linear wedge is degree 1 in (xi0, xi1) and degree 2 in zeta + const int cases[][2] = {{1, 2}, {2, 2}, {3, 11}}; + int failed = 0; + std::cout.precision(16); + for (auto const &c : cases) { + VolumeIntegrator vi(c[0], c[1]); + vi.process(fe); + const bool ok = isClose(vi.volume, exactVolume, 1e-13); + failed += !ok; + std::cout << "Integrator (" << c[0] << ", " << c[1] << "): volume " + << vi.volume << ", expected " << exactVolume + << (ok ? "" : " FAIL") << "\n"; + } + // one Gauss point in zeta is not enough for the tapered wedge + VolumeIntegrator low(1, 1); + low.process(fe); + const bool lowDiffers = !isClose(low.volume, exactVolume, 1e-6); + failed += !lowDiffers; + std::cout << "Integrator (1, 1): volume " << low.volume + << (lowDiffers ? " (not exact, as expected)" : " FAIL") << "\n"; + return failed; +} + +int runWedgeIntegrationChecks() { + int failed = checkEdgeRule(MeshField::EdgeIntegration::N1()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N2()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N3()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N4()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N5()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N6()); + failed += checkWedgeRule(1, 1, 1, 1); + failed += checkWedgeRule(2, 2, 3, 2); + failed += checkWedgeRule(2, 5, 3, 3); + failed += checkWedgeRule(4, 4, 6, 3); + failed += checkWedgeRule(3, 11, 6, 6); + failed += checkIntegrator(); + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runWedgeIntegrationChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +}