diff --git a/CMakeLists.txt b/CMakeLists.txt index b385901..ea27c4c 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(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(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(LinearWedgeShapeTest ./LinearWedgeShapeTest) +test_func(LinearWedgeElementTest ./LinearWedgeElementTest) +test_func(TriangleIntegrationTest ./TriangleIntegrationTest) +test_func(WedgeIntegrationTest ./WedgeIntegrationTest) 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_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..8489814 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -557,5 +557,113 @@ 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 + } +}; + } // namespace MeshField #endif diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 708e80a..776fa19 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; } 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; +}