diff --git a/CMakeLists.txt b/CMakeLists.txt index b385901..32e97ee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -220,6 +220,8 @@ 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(Interpolate test/testInterpolate.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -246,6 +248,8 @@ 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(Interpolate ./Interpolate) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error diff --git a/src/MeshField_Defines.hpp b/src/MeshField_Defines.hpp index 8467a0b..1c79020 100644 --- a/src/MeshField_Defines.hpp +++ b/src/MeshField_Defines.hpp @@ -13,7 +13,8 @@ enum Mesh_Topology { Tetrahedron, Hexahedron, Prism, - Pyramid + Pyramid, + Wedge }; static bool Debug = false; const Real MachinePrecision = 1e-15; diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 6fd4d47..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