From 19546061099504f389f80aa1a60e547534620596 Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Tue, 23 Sep 2025 13:01:10 -0400 Subject: [PATCH 1/9] Added integration test, added functionality for quadratic mesh Jacobian, began work on fixing ordering problems --- CMakeLists.txt | 2 + src/MeshField.hpp | 28 +++++ src/MeshField_Element.hpp | 28 ++++- src/MeshField_Integrate.hpp | 8 +- test/testElementJacobian2d.cpp | 4 + test/testElementJacobian3d.cpp | 3 + test/testIntegrator.cpp | 190 +++++++++++++++++++++++++++++++++ 7 files changed, 256 insertions(+), 7 deletions(-) create mode 100644 test/testIntegrator.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 99ff8b3..305664b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -172,6 +172,7 @@ meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(PointMapping test/testPointMapping.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) +meshfields_add_exe(IntegratorTest test/testIntegrator.cpp) if(MeshFields_USE_Cabana) meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp) @@ -190,6 +191,7 @@ test_func(CountIntegrator ./CountIntegrator) test_func(OmegahTriTests ./OmegahTriTests) test_func(PointMapping ./PointMapping) test_func(OmegahTetTest, ./OmegahTetTest) +test_func(IntegratorTest ./IntegratorTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 78757b6..52d6221 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -75,6 +75,16 @@ struct LinearTriangleToVertexField { getTopology() { return {MeshField::Triangle}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO triNodeIdx) const { + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + (Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) + + 2) % + 3; + return localVtxIdx; + } KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx, @@ -106,7 +116,17 @@ struct LinearTetrahedronToVertexField { getTopology() { return {MeshField::Tetrahedron}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO tetNodeIdx) const { + const auto tetDim = 3; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + (Omega_h::simplex_down_template(tetDim, vtxDim, tetNodeIdx, ignored) + + 3) % + 4; + return localVtxIdx; + } KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx, MeshField::LO tet, MeshField::Mesh_Topology topo) const { @@ -141,6 +161,10 @@ struct QuadraticTriangleToField { return {MeshField::Triangle}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO triNodeIdx) const { + return triNodeIdx; + } + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx, MeshField::LO tri, MeshField::Mesh_Topology topo) const { @@ -203,6 +227,10 @@ struct QuadraticTetrahedronToField { return {MeshField::Tetrahedron}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO tetNodeIdx) const { + return tetNodeIdx; + } + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx, MeshField::LO tet, MeshField::Mesh_Topology topo) const { diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 2197a0a..2aa501e 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -303,6 +303,20 @@ struct FieldElement { return Kokkos::View("foo", J.extent(0)); } + template + KOKKOS_INLINE_FUNCTION auto getGradients(Kokkos::View lc) const { + if constexpr (ShapeType::Order == 1) { + return shapeFn.getLocalGradients(); + } + else { + Kokkos::Array coord; + for (int i = 0; i < MeshEntDim + 1; ++i) { + coord[i] = lc(i); + } + return shapeFn.getLocalGradients(coord); + } + } + /** * @brief * Given an array of parametric coordinates 'localCoords', one per mesh @@ -380,22 +394,30 @@ struct FieldElement { // one matrix per point Kokkos::View res("result", numPts, MeshEntDim, MeshEntDim); Kokkos::deep_copy(res, 0.0); // initialize all entries to zero - // fill the views of node coordinates and node gradients Kokkos::View nodeCoords( "nodeCoords", numPts); Kokkos::View nodalGradients( "nodalGradients", numPts); - const auto grad = shapeFn.getLocalGradients(); Kokkos::parallel_for( numMeshEnts, KOKKOS_CLASS_LAMBDA(const int ent) { const auto vals = getNodeValues(ent); assert(vals.size() == MeshEntDim * ShapeType::numNodes); for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { + auto ptCoords = Kokkos::subview(localCoords, pt, Kokkos::ALL()); + const auto grad = getGradients(ptCoords); for (size_t node = 0; node < ShapeType::numNodes; node++) { + auto idx = elm2dof(node); + Kokkos::printf("\n"); for (size_t d = 0; d < MeshEntDim; d++) { + Kokkos::printf("%lf:", vals[node * MeshEntDim + d]); nodeCoords(pt, node, d) = vals[node * MeshEntDim + d]; - nodalGradients(pt, node, d) = grad[node * MeshEntDim + d]; + if constexpr (ShapeType::Order == 1) { + nodalGradients(pt, node, d) = grad[idx * MeshEntDim + d]; + } + else { + nodalGradients(pt, node, d) = grad[idx][d]; + } } } } diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 0c16d68..849638b 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -92,16 +92,16 @@ class TetrahedronIntegration : public EntityIntegration<4> { virtual std::vector> getPoints() const { return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011, - 0.138196601125011, 0.585410196624969}, + 0.138196601125011, 0.585410196624967}, 0.25 / 6.0), - IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011, + IntegrationPoint(Vector4{0.585410196624967, 0.138196601125011, 0.138196601125011, 0.138196601125011}, 0.25 / 6.0), - IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969, + IntegrationPoint(Vector4{0.138196601125011, 0.585410196624967, 0.138196601125011, 0.138196601125011}, 0.25 / 6.0), IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011, - 0.585410196624969, 0.138196601125011}, + 0.585410196624967, 0.138196601125011}, 0.25 / 6.0)}; } virtual int getAccuracy() const { return 2; } diff --git a/test/testElementJacobian2d.cpp b/test/testElementJacobian2d.cpp index 733e56d..48ff3b3 100644 --- a/test/testElementJacobian2d.cpp +++ b/test/testElementJacobian2d.cpp @@ -15,6 +15,10 @@ struct LinearTriangleToVertexField { return {MeshField::Triangle}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO triNodeIdx) const { + return triNodeIdx; + } + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx, MeshField::LO tri, MeshField::Mesh_Topology topo) const { diff --git a/test/testElementJacobian3d.cpp b/test/testElementJacobian3d.cpp index 6803433..9ec5e2c 100644 --- a/test/testElementJacobian3d.cpp +++ b/test/testElementJacobian3d.cpp @@ -14,6 +14,9 @@ struct LinearTetrahedronToVertexField { getTopology() const { return {MeshField::Tetrahedron}; } + KOKKOS_FUNCTION MeshField::LO operator()(MeshField::LO tetNodeIdx) const { + return tetNodeIdx; + } KOKKOS_FUNCTION MeshField::ElementToDofHolderMap operator()(MeshField::LO tetNodeIdx, MeshField::LO tetCompIdx, diff --git a/test/testIntegrator.cpp b/test/testIntegrator.cpp new file mode 100644 index 0000000..fad462b --- /dev/null +++ b/test/testIntegrator.cpp @@ -0,0 +1,190 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" //remove? +#include "MeshField_Fail.hpp" //remove? +#include "MeshField_For.hpp" //remove? +#include "MeshField_ShapeField.hpp" //remove? +#include "Omega_h_build.hpp" +#include "Omega_h_file.hpp" +#include "Omega_h_simplex.hpp" +#include +#include +#include +#include + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; + +template +class testIntegrator : public MeshField::Integrator { +private: + testIntegrator(){}; + +protected: + MeshField::Real integral; + FieldElement &fes; + int m; + int n; + int l; + +public: + MeshField::Real getIntegral() { return integral; } + testIntegrator(FieldElement &fes_in, int _m, int _n, int _l, int order) + : Integrator(order), integral(0), fes(fes_in), m(_m), n(_n), l(_l){}; + void atPoints(Kokkos::View p, + Kokkos::View w, + Kokkos::View dV) { + const auto globalCoords = MeshField::evaluate(fes, p, p.extent(0)); + const int _m = m; + const int _n = n; + const int _l = l; + Kokkos::parallel_reduce( + "integrate", p.extent(0), + KOKKOS_LAMBDA(const int &ent, MeshField::Real &integ) { + const auto x = globalCoords(ent, 0); + const auto y = globalCoords(ent, 1); + const auto z = globalCoords.extent(1) == 3 ? globalCoords(ent, 2) : 1; + integ += w(ent) * pow(x, _m) * pow(y, _n) * pow(z, _l) * dV(ent); + }, + integral); + } +}; + +template