Skip to content

Conversation

@jvmespark
Copy link

No description provided.

Copy link
Contributor

@cwsmith cwsmith left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you @jvmespark .

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).

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));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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

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

ReducedQuinticImplicitIntegrator<FieldElement> integrator(fes);
Copy link
Collaborator

@Joshua-Kloepfer Joshua-Kloepfer Oct 31, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Collaborator

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

};

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@jacobmerson
Copy link
Contributor

jacobmerson commented Nov 7, 2025

@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
Copy link
Collaborator

@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?

@jvmespark jvmespark marked this pull request as draft November 14, 2025 11:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants