From f49d13e2a13244f2e7dae27ca655b06b4f074ad3 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 14:43:51 -0400 Subject: [PATCH 01/10] add a degree-4 triangle integration rule --- src/MeshField_Integrate.hpp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index da80ba4..56164ef 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -74,11 +74,33 @@ class TriangleIntegration : public EntityIntegration<2> { } virtual int getAccuracy() const { return 2; } }; // end N2 + + // 6-point rule with positive weights, exact for polynomials of degree 4 + // (D. A. Dunavant, Int. J. Numer. Meth. Eng. 21, 1985). The points are the + // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: + // the only symmetric 4-point degree-3 rule has a negative weight. + class N4 : public Integration<2> { + public: + virtual int countPoints() const { return 6; } + virtual std::vector> getPoints() const { + const Real a1 = 0.44594849091596488632, w1 = 0.22338158967801146570; + const Real a2 = 0.091576213509770743460, w2 = 0.10995174365532186764; + return {IntegrationPoint(Vector2{a1, a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{1 - 2 * a1, a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{a1, 1 - 2 * a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{a2, a2}, w2 / 2.0, 2, 0), + IntegrationPoint(Vector2{1 - 2 * a2, a2}, w2 / 2.0, 2, 0), + IntegrationPoint(Vector2{a2, 1 - 2 * a2}, w2 / 2.0, 2, 0)}; + } + virtual int getAccuracy() const { return 4; } + }; // end N4 + virtual int countIntegrations() const { return 2; } 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]; } }; From a1c09552fb92f44e57a2b45caee42f3852968c3d Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:00:22 -0400 Subject: [PATCH 02/10] add TriangleIntegrationTest for the triangle quadrature rules --- test/testTriangleIntegration.cpp | 139 +++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 test/testTriangleIntegration.cpp diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp new file mode 100644 index 0000000..fd6c6b9 --- /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. +// 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); +} + +// 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 { + const int numPtsPerElem = p.extent(0) / fe.numMeshEnts; + const auto x = MeshField::evaluate(fe, p, numPtsPerElem); + const Real a = px, b = py; + Real sum = 0; + Kokkos::parallel_reduce( + "integrateMonomial", w.extent(0), + KOKKOS_LAMBDA(const int i, Real &lsum) { + lsum += + w(i) * 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; +} From 9af321daafcdcd32dfc964685ce2575b673faaf0 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:01:03 -0400 Subject: [PATCH 03/10] remove comment --- test/testTriangleIntegration.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp index fd6c6b9..5749068 100644 --- a/test/testTriangleIntegration.cpp +++ b/test/testTriangleIntegration.cpp @@ -15,8 +15,6 @@ // for order 5, // 3. Integrator::process with orders 3 and 4 integrates monomials of that // degree exactly over a triangulated unit square. -// Failures are counted explicitly so the test also fails in builds where -// assert() is disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; From 8205b1a53a73de60f98a01a073b213f25cce78fc Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:20:32 -0400 Subject: [PATCH 04/10] Register TriangleIntegrationTest in CMake --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 45dc805..e8efcf2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -213,6 +213,7 @@ meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) +meshfields_add_exe(TriangleIntegrationTest test/testTriangleIntegration.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -235,6 +236,7 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) +test_func(TriangleIntegrationTest ./TriangleIntegrationTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) From 995db181f36093acb9cccd0b08db5e37b600a413 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:40:49 -0400 Subject: [PATCH 05/10] set countIntegrations() to return 3 --- src/MeshField_Integrate.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 56164ef..5dcb2b1 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -95,7 +95,7 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 4; } }; // end N4 - virtual int countIntegrations() const { return 2; } + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; static N2 i2; From 59e41b120035cd338c58d778a56e823380b64556 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 14:30:42 -0400 Subject: [PATCH 06/10] add reference --- src/MeshField_Integrate.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 5dcb2b1..f425889 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,10 +75,13 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 - // 6-point rule with positive weights, exact for polynomials of degree 4 - // (D. A. Dunavant, Int. J. Numer. Meth. Eng. 21, 1985). The points are the + // 6-point rule with positive weights, exact for polynomials of degree 4. The points are the // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: // the only symmetric 4-point degree-3 rule has a negative weight. + // + // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature + // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). + // class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } From e542ba04cbc5d2729386199186e1b6c7afe53666 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 15:26:30 -0400 Subject: [PATCH 07/10] copied N4 rule from apfIntegrate.cc --- src/MeshField_Integrate.hpp | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index f425889..c2abdec 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,29 +75,33 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 - // 6-point rule with positive weights, exact for polynomials of degree 4. The points are the - // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: - // the only symmetric 4-point degree-3 rule has a negative weight. - // - // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature - // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). - // +// 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 p = 4 rules of from the Appendix II (p. 1140), +// copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). +// weights sum to one, so 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 { - const Real a1 = 0.44594849091596488632, w1 = 0.22338158967801146570; - const Real a2 = 0.091576213509770743460, w2 = 0.10995174365532186764; - return {IntegrationPoint(Vector2{a1, a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{1 - 2 * a1, a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{a1, 1 - 2 * a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{a2, a2}, w2 / 2.0, 2, 0), - IntegrationPoint(Vector2{1 - 2 * a2, a2}, w2 / 2.0, 2, 0), - IntegrationPoint(Vector2{a2, 1 - 2 * a2}, w2 / 2.0, 2, 0)}; + 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; From 7ef97e15b8f06a820232007fc53c2b80efde2139 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:01:16 -0400 Subject: [PATCH 08/10] Clean up the triangle N4 reference comment --- src/MeshField_Integrate.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index c2abdec..ceee260 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -82,6 +82,7 @@ class TriangleIntegration : public EntityIntegration<2> { // copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). // weights sum to one, so 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; } @@ -101,7 +102,7 @@ class TriangleIntegration : public EntityIntegration<2> { } 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; From 7fc9e565e62c2c34d4f3963934156333fb49bca5 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:33:02 -0400 Subject: [PATCH 09/10] Fix the N4 comment wording and indentation --- src/MeshField_Integrate.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index ceee260..e370fe0 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,14 +75,13 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 -// 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 p = 4 rules of from the Appendix II (p. 1140), -// copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). -// weights sum to one, so divided by 2, the area of the reference -// triangle. The p = 3 rule has a negative weight, so order 3 uses N4. - + // 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; } From 82b6e1b569e47c1d434a4d5c6be445f9949ec59f Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 11:07:42 -0400 Subject: [PATCH 10/10] Adapt TriangleIntegrationTest to per-element integration points --- test/testTriangleIntegration.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp index 5749068..8aacc8d 100644 --- a/test/testTriangleIntegration.cpp +++ b/test/testTriangleIntegration.cpp @@ -76,15 +76,17 @@ class MonomialIntegrator : public MeshField::Integrator { : MeshField::Integrator(order), fe(feIn), px(pxIn), py(pyIn) {} void atPoints(Kokkos::View p, Kokkos::View w, Kokkos::View dV) override { - const int numPtsPerElem = p.extent(0) / fe.numMeshEnts; - const auto x = MeshField::evaluate(fe, p, numPtsPerElem); + // 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", w.extent(0), + "integrateMonomial", dV.extent(0), KOKKOS_LAMBDA(const int i, Real &lsum) { - lsum += - w(i) * dV(i) * Kokkos::pow(x(i, 0), a) * Kokkos::pow(x(i, 1), b); + lsum += w(i % numPtsPerElem) * dV(i) * Kokkos::pow(x(i, 0), a) * + Kokkos::pow(x(i, 1), b); }, sum); integral += sum;