Reduced Quintic Implicit by jvmespark · Pull Request #77 · SCOREC/meshFields

jvmespark · 2025-10-30T07:34:00Z

No description provided.

cwsmith

In addition to the comments below, I have a a few requests:

Please add a description in the PR. If you got the shape function/gradient/quadrature points from literature please reference it.
The CMakeLists.txt file needs to be modified to add the new test.
Please add a test to exercise evaluate(...) (given user defined parametric coords in each element, return the field value at those locations). As with the integration test, I'd strongly prefer a test that we know the result for and can check against it (as opposed to a simpler 'it doesn't crash' test).

cwsmith · 2025-10-30T16:20:52Z

test/testReducedQuinticImplicitIntegrator.cpp

+  const auto family = OMEGA_H_SIMPLEX;
+  auto len = 1.0;
+  return Omega_h::build_box(world, family, len, len, len,
+                            2, 2, (dim == 3 ? 2 : 0));


I think the template dim can be removed and the third dimension hardcoded as zero.

cwsmith · 2025-10-30T16:23:21Z

test/testReducedQuinticImplicitIntegrator.cpp

+  {
+    auto mesh2D = createMesh<2>(lib);
+    MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController> omf2D(mesh2D);
+    runReducedQuinticImplicit<MeshField::KokkosController>(mesh2D, omf2D);


Is there an expected result from the integration? I'd much prefer a test that we know the result for.

There is probably some test cases in M3DC1 that I can dig up and share with James.

Joshua-Kloepfer · 2025-10-31T19:14:44Z

src/MeshField_Integrate.hpp

+                           0.1167862757264 / 2.0),
+          IntegrationPoint(Vector3{0.249286745, 0.249286745, 0.501426509},
+                           0.1167862757264 / 2.0),
+


I think some of the points here do not have enough precision and fail summing to 1.

@jvmespark This looks like it is the standard Gauss-Legendre points. You can copy them from here with more precision: https://github.com/SCOREC/core/blob/6a972af6daceec626e589b4d4f63acc22e27317f/apf/apfIntegrate.cc#L216

Alternatively, it looks like we may be able to use this script in APF to calculate to machine precision: https://github.com/SCOREC/core/blob/6a972af6daceec626e589b4d4f63acc22e27317f/apf/apfPolyBasis1D.cc

Joshua-Kloepfer · 2025-10-31T19:17:21Z

test/testReducedQuinticImplicitIntegrator.cpp

+  const auto [shp, map] = MeshField::Omegah::getReducedQuinticImplicitElement(mesh);
+  MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
+
+  ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);


Suggested change

ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);

ReducedQuinticImplicitIntegrator integrator(fes);

Either do not put in template or you have to include the templating of FieldElement. ie MeshField::FieldElement<templating>.

Joshua-Kloepfer · 2025-10-31T19:19:44Z

Can you add
meshfields_add_exe(ReducedQuinticImplicitIntegratorTest test/testReducedQuinticImplicitIntegrator.cpp)
test_func(ReducedQuinticImplicitIntegratorTest ./ReducedQuinticImplicitIntegratorTest)
to CMakeLists.txt

Joshua-Kloepfer · 2025-10-31T19:20:38Z

src/MeshField.hpp

+};
+
 template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
  static_assert(ShapeOrder == 1 || ShapeOrder == 2);


Suggested change

static_assert(ShapeOrder == 1 || ShapeOrder == 2);

static_assert(ShapeOrder == 1 || ShapeOrder == 2 || ShapeOrder == 5);

jacobmerson · 2025-11-07T18:25:15Z

@cwsmith @Joshua-Kloepfer One thing that we noticed is that in core the integration only stores n-1 coordinates where n is the number of vertices. This is why Vector3 can be used for Tets (it has 4 parametric coordinates). Do we want to make the integration point definitions consistent between triangles/tets? If so, we should either ignore the last coordinate and compute it with $$\xi_n = 1-\sum_{i=1}^{num nodes-1} \xi_i$$ or the storage needs to be able to accommodate 4 entries.

Joshua-Kloepfer · 2025-11-12T19:05:29Z

@jacobmerson Currently in our code, we store all n coordinates so for Tets we decided to store all 4 parametric coordinates (https://github.com/SCOREC/meshFields/blob/jk/FixIntegration/src/MeshField_Integrate.hpp#L85) and in triangles we store 3 (https://github.com/SCOREC/meshFields/blob/jk/FixIntegration/src/MeshField_Integrate.hpp#L85). We felt that this was the right approach. Is there somewhere else where we are not storing all the entries?

jacobmerson · 2026-02-13T17:01:23Z

test/testReducedQuinticImplicitEvaluate.cpp

+  auto& rqShape = element.shp;

-  MeshField::Vector3 xi = {1.0/3.0, 1.0/3.0, 1.0/3.0};
+  constexpr int nn = static_cast<int>(MeshField::ReducedQuinticImplicitShape::numNodes);


why does this require a static cast to int?

jacobmerson · 2026-02-13T17:01:56Z

test/testReducedQuinticImplicitEvaluate.cpp

-  auto dN = rqShape.getLocalGradients(xi);
+  // Evaluate shape functions and local gradients at the centroid.
+  const auto N  = rqShape.getValues(xi);
+  (void)rqShape.getLocalGradients(xi); // computed for coverage; not asserted here


? what do you mean here?

jacobmerson · 2026-02-13T17:02:50Z

test/testReducedQuinticImplicitEvaluate.cpp

-  for (int i = 0; i < nn; ++i) {
-    sumN += N[i];
-  }
+  for (int i = 0; i < nn; ++i) sumN += static_cast<double>(N[i]);


why are you static casting to double? Since this is on CPU, you can use std::accumulate

jacobmerson · 2026-02-13T17:04:08Z

test/testReducedQuinticImplicitIntegrator.cpp

-    printf("[ReducedQuinticImplicit] Integrated Value = %.12e\n", total);
+                Kokkos::View<MeshField::Real *> /*dV*/) override
+  {
+    const int npts = static_cast<int>(w.extent(0));


can you use const auto here? why do you need to cast to int?

jacobmerson · 2026-02-13T17:05:04Z

test/testReducedQuinticImplicitIntegrator.cpp

+      "Integrate_f_bary_to_phys", npts,
+      KOKKOS_LAMBDA(const int ip, double& acc) {
+        // Barycentric weights (l0,l1,l2) corresponding to the element's vertices (v0,v1,v2).
+        const double l0 = static_cast<double>(p(ip, 0));


I expect the result of p to be double? why is static cast neeeded? Use const auto instead

jacobmerson · 2026-02-13T17:05:37Z

test/testReducedQuinticImplicitIntegrator.cpp

+        const double y = l0 * y0 + l1 * y1 + l2 * y2;
+
+        // Accumulate f(x,y) * w. We assume w already includes dV.
+        acc += f_analytic(x, y) * static_cast<double>(w(ip)); // w includes dV


static_cast?

jacobmerson · 2026-02-13T17:05:43Z

test/testReducedQuinticImplicitIntegrator.cpp

+    Kokkos::parallel_reduce(
+      "Integrate_1_bary_to_phys", npts,
+      KOKKOS_LAMBDA(const int ip, double& acc) {
+        acc += static_cast<double>(w(ip));


static_cast?

jacobmerson · 2026-02-13T17:08:47Z

src/MeshField.hpp

+    if (triNodeIdx == 0) {
+      osh_ent = triVerts[3*tri + 0];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 1) {
+      osh_ent = triVerts[3*tri + 1];
+      nodeInHolder = 0;
+    }
+    else if (triNodeIdx == 2) {
+      osh_ent = triVerts[3*tri + 2];
+      nodeInHolder = 0;
+    }


can't all of this be one statement? i.e., osh_ent = triVerts[3*tri + triNodeIdx]

jacobmerson · 2026-02-13T17:11:27Z

src/MeshField.hpp

+    }
+    else if (triNodeIdx >= 3 && triNodeIdx <= 6) {
+      const MeshField::LO v0 = triVerts[3*tri + 0];
+      const MeshField::LO v1 = triVerts[3*tri + 1];


why is this needed in this branch of the if statement?

jacobmerson · 2026-02-13T17:14:21Z

src/MeshField_Shape.hpp

+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);


If you need to do this you should use an IIFE, don't bother naming the lambda.

jacobmerson · 2026-02-13T17:15:33Z

src/MeshField_Shape.hpp

+    for (int i = 0; i <= 5; ++i) {
+      for (int j = 0; j <= 5 - i; ++j) {
+        int k = 5 - i - j;
+        double coeff = f5 / (fact(i) * fact(j) * fact(k));


this calculation needs to be explained in the comment. One issue I see is that you are recomputing the factorials many times. Would it make more sense to have a table of the coefficients?

jacobmerson · 2026-02-13T17:16:13Z

src/MeshField_Shape.hpp

+    auto fact = [](int n) {
+      double r = 1.0;
+      for (int i = 2; i <= n; ++i) r *= double(i);
+      return r;
+    };
+    const double f5 = fact(5);


why have you repeated the code?

jacobmerson · 2026-02-13T17:42:33Z

test/testReducedQuinticImplicitIntegrator.cpp

+    Omega_h::Mesh mesh(&lib);
+    Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, verts, coords);
+
+    runTest<MeshField::KokkosController>(mesh);


have run test return a bool or error code that you can use to set the return value of main. This can be used with ctest to know if the test passed/failed.

jacobmerson · 2026-02-13T17:48:39Z

test/testReducedQuinticImplicitEvaluate.cpp

+  }
+
+  std::cout << "[PASS] ReducedQuinticImplicit evaluator test (partition of unity)\n";
+}


You are missing a test case where you evaluate the shape function for example at 1/3, 1/3, 1/3 for up to a 4th order polynomial. You should be able to get an exact value.

init

01a16db

cwsmith requested changes Oct 30, 2025

View reviewed changes

cwsmith assigned Joshua-Kloepfer Oct 30, 2025

Joshua-Kloepfer reviewed Oct 31, 2025

View reviewed changes

test evaluation

5dad2c9

jvmespark added 3 commits November 14, 2025 00:37

updated precision vals

ae07b45

correctly sums to one

3622fa0

integral test

06d337d

jvmespark marked this pull request as draft November 14, 2025 11:47

jvmespark added 4 commits November 14, 2025 06:47

test quintic eval

01da95c

test eval & integrator, changes to field and shape

c5d418b

fix integrator test

ee9e7a0

updated evaluate and integrate tests for reduced quintic implicit

de30493

jacobmerson suggested changes Feb 13, 2026

View reviewed changes

jacobmerson reviewed Feb 13, 2026

View reviewed changes

jvmespark added 3 commits February 17, 2026 15:36

new integrate test, tabled factorials, iife, simplified troNodeIdx leq 2

4b36ffa

new eval test, precomputed table

3a508d9

4th order eval test

6526e8e

	ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);
	ReducedQuinticImplicitIntegrator integrator(fes);

	static_assert(ShapeOrder == 1 \|\| ShapeOrder == 2);
	static_assert(ShapeOrder == 1 \|\| ShapeOrder == 2 \|\| ShapeOrder == 5);

Comments

Conversation

jvmespark commented Oct 30, 2025

Uh oh!

cwsmith left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Joshua-Kloepfer Oct 31, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Joshua-Kloepfer commented Oct 31, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jacobmerson commented Nov 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Joshua-Kloepfer commented Nov 12, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Joshua-Kloepfer Oct 31, 2025 •

edited

Loading

jacobmerson commented Nov 7, 2025 •

edited

Loading