Skip to content

tetrahedron integration #68

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
cwsmith merged 6 commits into from
Aug 18, 2025
Merged

tetrahedron integration #68

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
655c5f3
tetrahedron integration
Aug 15, 2025
01f3e08
formatting
Aug 15, 2025
3d90884
format
Aug 15, 2025
ed8c4e2
removed constexpr
Aug 15, 2025
7acfb1f
changed templating to be topology and changed template paramater name
Aug 18, 2025
293589a
changed template name
Aug 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,6 +166,7 @@ meshfields_add_exe(SerializationTests test/testSerialize.cpp)
meshfields_add_exe(ElementTests test/testElement.cpp)
meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp)
meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
meshfields_add_exe(ElementJacobian3d test/testElementJacobian3d.cpp)
meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
Expand All @@ -184,6 +185,7 @@ test_func(SerializationTests ./SerializationTests)
test_func(ElementTests ./ElementTests)
test_func(ElementJacobian1d ./ElementJacobian1d)
test_func(ElementJacobian2d ./ElementJacobian2d)
test_func(ElementJacobian3d ./ElementJacobian3d)
test_func(CountIntegrator ./CountIntegrator)
test_func(OmegahTriTests ./OmegahTriTests)
test_func(PointMapping ./PointMapping)
Expand Down
21 changes: 15 additions & 6 deletions src/MeshField_Element.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,9 +256,17 @@ struct FieldElement {
"tangent vectors" in 3D, the volume of their
parallelpiped, which is the differential volume
of the coordinate field */
fail("getJacobianDeterminant doesn't yet support 3d. "
"The given matrices have dimension %lu x %lu \n",
J.extent(0), J.extent(1));
Kokkos::View<Real *> determinants("3dJacobianDeterminants", J.extent(0));
Kokkos::parallel_for(
J.extent(0), KOKKOS_LAMBDA(const int i) {
const auto Ji = Kokkos::subview(J, i, Kokkos::ALL(), Kokkos::ALL());
const auto cofactorSum =
Ji(0, 0) * (Ji(1, 1) * Ji(2, 2) - Ji(1, 2) * Ji(2, 1)) -
Ji(0, 1) * (Ji(1, 0) * Ji(2, 2) - Ji(1, 2) * Ji(2, 0)) +
Ji(0, 2) * (Ji(1, 0) * Ji(2, 1) - Ji(1, 1) * Ji(2, 0));
determinants(i) = cofactorSum;
});
return determinants;
}
if (dimension == 2) {
if (J.extent(1) != 2) {
Expand Down Expand Up @@ -349,8 +357,9 @@ struct FieldElement {
fail("The input array of offsets must have size = %zu\n",
numMeshEnts + 1);
}
if (MeshEntDim != 1 && MeshEntDim != 2) {
fail("getJacobians only currently supports 1d and 2d meshes. Input mesh "
if (MeshEntDim != 1 && MeshEntDim != 2 && MeshEntDim != 3) {
fail("getJacobians only currently supports 1d, 2d, and 3d meshes. Input "
"mesh "
"has %zu dimensions.\n",
numMeshEnts);
}
Expand All @@ -366,7 +375,7 @@ struct FieldElement {
}
});
return res;
} else if constexpr (MeshEntDim == 2) {
} else if constexpr (MeshEntDim == 2 || MeshEntDim == 3) {
const auto numPts = MeshFieldUtil::getLastValue(offsets);
// one matrix per point
Kokkos::View<Real ***> res("result", numPts, MeshEntDim, MeshEntDim);
Expand Down
117 changes: 80 additions & 37 deletions src/MeshField_Integrate.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,50 +10,49 @@

namespace MeshField {
// directly copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.[h|cc]
struct IntegrationPoint {
IntegrationPoint(Vector3 const &p, double w) : param(p), weight(w) {}
Vector3 param;
template <size_t dim> struct IntegrationPoint {
IntegrationPoint(Kokkos::Array<Real, dim> const &p, double w)
: param(p), weight(w) {}
Kokkos::Array<Real, dim> param;
double weight;
};

class Integration {
template <size_t dim> class Integration {
public:
virtual ~Integration() {}
virtual int countPoints() const = 0;
virtual std::vector<IntegrationPoint> getPoints() const = 0;
virtual std::vector<IntegrationPoint<dim>> getPoints() const = 0;
virtual int getAccuracy() const = 0;
};

class EntityIntegration {
template <size_t dim> class EntityIntegration {
public:
virtual ~EntityIntegration() {}
Integration const *getAccurate(int minimumAccuracy) const {
Integration<dim> const *getAccurate(int minimumAccuracy) const {
int n = countIntegrations();
for (int i = 0; i < n; ++i) {
Integration const *integration = getIntegration(i);
Integration<dim> const *integration = getIntegration(i);
if (integration->getAccuracy() >= minimumAccuracy)
return integration;
}
return NULL;
}
virtual int countIntegrations() const = 0;
virtual Integration const *getIntegration(int i) const = 0;
virtual Integration<dim> const *getIntegration(int i) const = 0;
};

class TriangleIntegration : public EntityIntegration {
class TriangleIntegration : public EntityIntegration<3> {
public:
class N1 : public Integration {
class N1 : public Integration<3> {
public:
virtual int countPoints() const { return 1; }
virtual std::vector<IntegrationPoint> getPoints() const {
virtual std::vector<IntegrationPoint<3>> getPoints() const {
return {IntegrationPoint(Vector3{1. / 3., 1. / 3., 1. / 3.}, 1.0 / 2.0)};
}
virtual int getAccuracy() const { return 1; }
}; // end N1
class N2 : public Integration {
class N2 : public Integration<3> {
public:
virtual int countPoints() const { return 3; }
virtual std::vector<IntegrationPoint> getPoints() const {
virtual std::vector<IntegrationPoint<3>> getPoints() const {
return {IntegrationPoint(Vector3{0.666666666666667, 0.166666666666667,
0.166666666666667},
1. / 3. / 2.0),
Expand All @@ -67,31 +66,78 @@ class TriangleIntegration : public EntityIntegration {
virtual int getAccuracy() const { return 2; }
}; // end N2
virtual int countIntegrations() const { return 2; }
virtual Integration const *getIntegration(int i) const {
virtual Integration<3> const *getIntegration(int i) const {
static N1 i1;
static N2 i2;
static Integration *integrations[2] = {&i1, &i2};
static Integration<3> *integrations[2] = {&i1, &i2};
return integrations[i];
}
};

std::shared_ptr<EntityIntegration> const getIntegration(Mesh_Topology topo) {
if (topo != Triangle) {
fail("getIntegration only supports triangles (Mesh_Topology::Triangle)\n");
class TetrahedronIntegration : public EntityIntegration<4> {
public:
class N1 : public Integration<4> {
public:
virtual int countPoints() const { return 1; }
virtual std::vector<IntegrationPoint<4>> getPoints() const {
return {IntegrationPoint(Vector4{0.25, 0.25, 0.25, 0.25}, 1.0 / 6.0)};
}
virtual int getAccuracy() const { return 1; }
};
class N2 : public Integration<4> {
public:
virtual int countPoints() const { return 4; }
virtual std::vector<IntegrationPoint<4>> getPoints() const {

return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
0.138196601125011, 0.585410196624969},
0.25 / 6.0),
IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011,
0.138196601125011, 0.138196601125011},
0.25 / 6.0),
IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969,
0.138196601125011, 0.138196601125011},
0.25 / 6.0),
IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
0.585410196624969, 0.138196601125011},
0.25 / 6.0)};
}
virtual int getAccuracy() const { return 2; }
};
virtual int countIntegrations() const { return 2; }
virtual Integration<4> const *getIntegration(int i) const {
static N1 i1;
static N2 i2;
static Integration<4> *integrations[2] = {&i1, &i2};
return integrations[i];
}
};
template <size_t dim>
std::shared_ptr<EntityIntegration<dim>> const
getIntegration(Mesh_Topology topo) {
if constexpr (dim == 3) {
if (topo == Triangle) {
return std::make_shared<TriangleIntegration>();
}
} else if constexpr (dim == 4) {
if (topo == Tetrahedron) {
return std::make_shared<TetrahedronIntegration>();
}
}
return std::make_shared<TriangleIntegration>();
fail("getIntegration does not support given topology\n");
return nullptr;
}

std::vector<IntegrationPoint> getIntegrationPoints(Mesh_Topology topo,
int order) {
auto ip = getIntegration(topo)->getAccurate(order)->getPoints();
template <size_t dim>
std::vector<IntegrationPoint<dim>> getIntegrationPoints(Mesh_Topology topo,
int order) {
auto ip = getIntegration<dim>(topo)->getAccurate(order)->getPoints();
return ip;
}

template <typename FieldElement>
Kokkos::View<MeshField::Real **>
getIntegrationPointLocalCoords(FieldElement &fes,
std::vector<IntegrationPoint> ip) {
Kokkos::View<MeshField::Real **> getIntegrationPointLocalCoords(
FieldElement &fes,
std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
const auto numPtsPerElm = ip.size();
const auto numMeshEnts = fes.numMeshEnts;
const auto meshEntDim = fes.MeshEntDim;
Expand All @@ -115,9 +161,9 @@ getIntegrationPointLocalCoords(FieldElement &fes,
}

template <typename FieldElement>
Kokkos::View<Real *>
getIntegrationPointWeights(FieldElement &fes,
std::vector<IntegrationPoint> ip) {
Kokkos::View<Real *> getIntegrationPointWeights(
FieldElement &fes,
std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
const auto numPtsPerElm = ip.size();
const auto numMeshEnts = fes.numMeshEnts;
const auto meshEntDim = fes.MeshEntDim;
Expand Down Expand Up @@ -195,12 +241,9 @@ class Integrator {
* */
template <typename FieldElement> void process(FieldElement &fes) {
const auto topo = fes.elm2dof.getTopology();
if (topo[0] != MeshField::Triangle) {
throw std::invalid_argument(
"Integrator::process only supports triangles.");
}
pre();
auto ip = getIntegrationPoints(topo[0], order);
auto ip =
getIntegrationPoints<FieldElement::MeshEntDim + 1>(topo[0], order);
auto localCoords = getIntegrationPointLocalCoords(fes, ip);
auto weights = getIntegrationPointWeights(fes, ip);
auto dV = getJacobianDeterminants(fes, localCoords, ip.size());
Expand Down
40 changes: 28 additions & 12 deletions test/testCountIntegrator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,12 @@

using ExecutionSpace = Kokkos::DefaultExecutionSpace;
using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;

Omega_h::Mesh createMeshTri18(Omega_h::Library &lib) {
template <size_t dim> Omega_h::Mesh createMesh(Omega_h::Library &lib) {
auto world = lib.world();
const auto family = OMEGA_H_SIMPLEX;
auto len = 1.0;
return Omega_h::build_box(world, family, len, len, 0.0, 3, 3, 0);
const auto numEnts3d = (dim == 3 ? 3 : 0);
return Omega_h::build_box(world, family, len, len, len, 3, 3, numEnts3d);
}

template <typename AnalyticFunction, typename ShapeField>
Expand Down Expand Up @@ -57,13 +57,19 @@ class CountIntegrator : public MeshField::Integrator {
}
};

template <template <typename...> typename Controller>
template <template <typename...> typename Controller, size_t dim>
void doRun(Omega_h::Mesh &mesh,
MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf) {
MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
const auto ShapeOrder = 1;
auto field = omf.getCoordField();
const auto [shp, map] =
MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
auto shapeSet = [&]() {
if constexpr (dim == 3) {
return MeshField::Omegah::getTetrahedronElement<ShapeOrder>(mesh);
} else {
return MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
}
};
const auto [shp, map] = shapeSet();
MeshField::FieldElement fes(mesh.nelems(), field, shp, map);

CountIntegrator countInt(fes);
Expand All @@ -74,18 +80,28 @@ void doRun(Omega_h::Mesh &mesh,
int main(int argc, char **argv) {
Kokkos::initialize(argc, argv);
auto lib = Omega_h::Library(&argc, &argv);
auto mesh = createMeshTri18(lib);
#ifdef MESHFIELDS_ENABLE_CABANA
{
auto mesh2D = createMesh<2>(lib);
auto mesh3D = createMesh<3>(lib);
MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::CabanaController>
omf(mesh);
doRun<MeshField::CabanaController>(mesh, omf);
omf2D(mesh2D);
doRun<MeshField::CabanaController>(mesh2D, omf2D);
MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::CabanaController>
omf3D(mesh3D);
doRun<MeshField::CabanaController>(mesh3D, omf3D);
}
#endif
{
auto mesh2D = createMesh<2>(lib);
auto mesh3D = createMesh<3>(lib);
MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController>
omf(mesh);
doRun<MeshField::KokkosController>(mesh, omf);
omf2D(mesh2D);
doRun<MeshField::KokkosController>(mesh2D, omf2D);
MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::KokkosController>
omf3D(mesh3D);
doRun<MeshField::KokkosController>(mesh3D, omf3D);
}
Kokkos::finalize();
return 0;
}
2 changes: 1 addition & 1 deletion test/testElementJacobian2d.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ void triJacobian() {
LinearTriangleToVertexField());

Kokkos::View<MeshField::Real *[3]> lc("localCoords", 1);
Kokkos::deep_copy(lc, 1.0 / 2);
Kokkos::deep_copy(lc, 1.0 / 3);
const auto numPtsPerElement = 1;
const auto J = MeshField::getJacobians(f, lc, numPtsPerElement);
const auto J_h =
Expand Down
Loading
Loading