Skip to content

Commit 192472c

Browse files
cwsmithcwsmith
authored
Merge pull request #68 from SCOREC/jk/tetIntegrate
Adding functionality for integrating tetrahedron (linear and quadratic).
2 parents 164bd2b + 293589a commit 192472c

File tree

7 files changed

+245
-65
lines changed

7 files changed

+245
-65
lines changed

‎CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -166,6 +166,7 @@ meshfields_add_exe(SerializationTests test/testSerialize.cpp)
166166
meshfields_add_exe(ElementTests test/testElement.cpp)
167167
meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp)
168168
meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
169+
meshfields_add_exe(ElementJacobian3d test/testElementJacobian3d.cpp)
169170
meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
170171
meshfields_add_exe(OmegahTriTests test/testOmegahTri.cpp)
171172
meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
@@ -184,6 +185,7 @@ test_func(SerializationTests ./SerializationTests)
184185
test_func(ElementTests ./ElementTests)
185186
test_func(ElementJacobian1d ./ElementJacobian1d)
186187
test_func(ElementJacobian2d ./ElementJacobian2d)
188+
test_func(ElementJacobian3d ./ElementJacobian3d)
187189
test_func(CountIntegrator ./CountIntegrator)
188190
test_func(OmegahTriTests ./OmegahTriTests)
189191
test_func(PointMapping ./PointMapping)

‎src/MeshField.hpp‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -71,8 +71,8 @@ struct LinearTriangleToVertexField {
7171
}
7272
}
7373

74-
KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
75-
getTopology() const {
74+
static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
75+
getTopology() {
7676
return {MeshField::Triangle};
7777
}
7878

@@ -102,8 +102,8 @@ struct LinearTetrahedronToVertexField {
102102
__func__);
103103
}
104104
}
105-
KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
106-
getTopology() const {
105+
static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
106+
getTopology() {
107107
return {MeshField::Tetrahedron};
108108
}
109109

@@ -136,8 +136,8 @@ struct QuadraticTriangleToField {
136136
}
137137
}
138138

139-
KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
140-
getTopology() const {
139+
static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
140+
getTopology() {
141141
return {MeshField::Triangle};
142142
}
143143

@@ -198,8 +198,8 @@ struct QuadraticTetrahedronToField {
198198
}
199199
}
200200

201-
KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
202-
getTopology() const {
201+
static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
202+
getTopology() {
203203
return {MeshField::Tetrahedron};
204204
}
205205

‎src/MeshField_Element.hpp‎

Lines changed: 15 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -256,9 +256,17 @@ struct FieldElement {
256256
"tangent vectors" in 3D, the volume of their
257257
parallelpiped, which is the differential volume
258258
of the coordinate field */
259-
fail("getJacobianDeterminant doesn't yet support 3d. "
260-
"The given matrices have dimension %lu x %lu \n",
261-
J.extent(0), J.extent(1));
259+
Kokkos::View<Real *> determinants("3dJacobianDeterminants", J.extent(0));
260+
Kokkos::parallel_for(
261+
J.extent(0), KOKKOS_LAMBDA(const int i) {
262+
const auto Ji = Kokkos::subview(J, i, Kokkos::ALL(), Kokkos::ALL());
263+
const auto cofactorSum =
264+
Ji(0, 0) * (Ji(1, 1) * Ji(2, 2) - Ji(1, 2) * Ji(2, 1)) -
265+
Ji(0, 1) * (Ji(1, 0) * Ji(2, 2) - Ji(1, 2) * Ji(2, 0)) +
266+
Ji(0, 2) * (Ji(1, 0) * Ji(2, 1) - Ji(1, 1) * Ji(2, 0));
267+
determinants(i) = cofactorSum;
268+
});
269+
return determinants;
262270
}
263271
if (dimension == 2) {
264272
if (J.extent(1) != 2) {
@@ -349,8 +357,9 @@ struct FieldElement {
349357
fail("The input array of offsets must have size = %zu\n",
350358
numMeshEnts + 1);
351359
}
352-
if (MeshEntDim != 1 && MeshEntDim != 2) {
353-
fail("getJacobians only currently supports 1d and 2d meshes. Input mesh "
360+
if (MeshEntDim != 1 && MeshEntDim != 2 && MeshEntDim != 3) {
361+
fail("getJacobians only currently supports 1d, 2d, and 3d meshes. Input "
362+
"mesh "
354363
"has %zu dimensions.\n",
355364
numMeshEnts);
356365
}
@@ -366,7 +375,7 @@ struct FieldElement {
366375
}
367376
});
368377
return res;
369-
} else if constexpr (MeshEntDim == 2) {
378+
} else if constexpr (MeshEntDim == 2 || MeshEntDim == 3) {
370379
const auto numPts = MeshFieldUtil::getLastValue(offsets);
371380
// one matrix per point
372381
Kokkos::View<Real ***> res("result", numPts, MeshEntDim, MeshEntDim);

‎src/MeshField_Integrate.hpp‎

Lines changed: 73 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -10,50 +10,51 @@
1010

1111
namespace MeshField {
1212
// directly copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.[h|cc]
13-
struct IntegrationPoint {
14-
IntegrationPoint(Vector3 const &p, double w) : param(p), weight(w) {}
15-
Vector3 param;
13+
template <size_t pointSize> struct IntegrationPoint {
14+
// template parameter pointSize specifies the length of the integration point
15+
// array for one point
16+
IntegrationPoint(Kokkos::Array<Real, pointSize> const &p, double w)
17+
: param(p), weight(w) {}
18+
Kokkos::Array<Real, pointSize> param;
1619
double weight;
1720
};
18-
19-
class Integration {
21+
template <size_t pointSize> class Integration {
2022
public:
2123
virtual ~Integration() {}
2224
virtual int countPoints() const = 0;
23-
virtual std::vector<IntegrationPoint> getPoints() const = 0;
25+
virtual std::vector<IntegrationPoint<pointSize>> getPoints() const = 0;
2426
virtual int getAccuracy() const = 0;
2527
};
26-
27-
class EntityIntegration {
28+
template <size_t pointSize> class EntityIntegration {
2829
public:
2930
virtual ~EntityIntegration() {}
30-
Integration const *getAccurate(int minimumAccuracy) const {
31+
Integration<pointSize> const *getAccurate(int minimumAccuracy) const {
3132
int n = countIntegrations();
3233
for (int i = 0; i < n; ++i) {
33-
Integration const *integration = getIntegration(i);
34+
Integration<pointSize> const *integration = getIntegration(i);
3435
if (integration->getAccuracy() >= minimumAccuracy)
3536
return integration;
3637
}
3738
return NULL;
3839
}
3940
virtual int countIntegrations() const = 0;
40-
virtual Integration const *getIntegration(int i) const = 0;
41+
virtual Integration<pointSize> const *getIntegration(int i) const = 0;
4142
};
4243

43-
class TriangleIntegration : public EntityIntegration {
44+
class TriangleIntegration : public EntityIntegration<3> {
4445
public:
45-
class N1 : public Integration {
46+
class N1 : public Integration<3> {
4647
public:
4748
virtual int countPoints() const { return 1; }
48-
virtual std::vector<IntegrationPoint> getPoints() const {
49+
virtual std::vector<IntegrationPoint<3>> getPoints() const {
4950
return {IntegrationPoint(Vector3{1. / 3., 1. / 3., 1. / 3.}, 1.0 / 2.0)};
5051
}
5152
virtual int getAccuracy() const { return 1; }
5253
}; // end N1
53-
class N2 : public Integration {
54+
class N2 : public Integration<3> {
5455
public:
5556
virtual int countPoints() const { return 3; }
56-
virtual std::vector<IntegrationPoint> getPoints() const {
57+
virtual std::vector<IntegrationPoint<3>> getPoints() const {
5758
return {IntegrationPoint(Vector3{0.666666666666667, 0.166666666666667,
5859
0.166666666666667},
5960
1. / 3. / 2.0),
@@ -67,31 +68,69 @@ class TriangleIntegration : public EntityIntegration {
6768
virtual int getAccuracy() const { return 2; }
6869
}; // end N2
6970
virtual int countIntegrations() const { return 2; }
70-
virtual Integration const *getIntegration(int i) const {
71+
virtual Integration<3> const *getIntegration(int i) const {
7172
static N1 i1;
7273
static N2 i2;
73-
static Integration *integrations[2] = {&i1, &i2};
74+
static Integration<3> *integrations[2] = {&i1, &i2};
7475
return integrations[i];
7576
}
7677
};
7778

78-
std::shared_ptr<EntityIntegration> const getIntegration(Mesh_Topology topo) {
79-
if (topo != Triangle) {
80-
fail("getIntegration only supports triangles (Mesh_Topology::Triangle)\n");
79+
class TetrahedronIntegration : public EntityIntegration<4> {
80+
public:
81+
class N1 : public Integration<4> {
82+
public:
83+
virtual int countPoints() const { return 1; }
84+
virtual std::vector<IntegrationPoint<4>> getPoints() const {
85+
return {IntegrationPoint(Vector4{0.25, 0.25, 0.25, 0.25}, 1.0 / 6.0)};
86+
}
87+
virtual int getAccuracy() const { return 1; }
88+
};
89+
class N2 : public Integration<4> {
90+
public:
91+
virtual int countPoints() const { return 4; }
92+
virtual std::vector<IntegrationPoint<4>> getPoints() const {
93+
94+
return {IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
95+
0.138196601125011, 0.585410196624969},
96+
0.25 / 6.0),
97+
IntegrationPoint(Vector4{0.585410196624969, 0.138196601125011,
98+
0.138196601125011, 0.138196601125011},
99+
0.25 / 6.0),
100+
IntegrationPoint(Vector4{0.138196601125011, 0.585410196624969,
101+
0.138196601125011, 0.138196601125011},
102+
0.25 / 6.0),
103+
IntegrationPoint(Vector4{0.138196601125011, 0.138196601125011,
104+
0.585410196624969, 0.138196601125011},
105+
0.25 / 6.0)};
106+
}
107+
virtual int getAccuracy() const { return 2; }
108+
};
109+
virtual int countIntegrations() const { return 2; }
110+
virtual Integration<4> const *getIntegration(int i) const {
111+
static N1 i1;
112+
static N2 i2;
113+
static Integration<4> *integrations[2] = {&i1, &i2};
114+
return integrations[i];
81115
}
82-
return std::make_shared<TriangleIntegration>();
116+
};
117+
template <Mesh_Topology topo> auto const getIntegration() {
118+
if constexpr (topo == Triangle) {
119+
return std::make_shared<TriangleIntegration>();
120+
} else if constexpr (topo == Tetrahedron) {
121+
return std::make_shared<TetrahedronIntegration>();
122+
}
123+
fail("getIntegration does not support given topology\n");
83124
}
84-
85-
std::vector<IntegrationPoint> getIntegrationPoints(Mesh_Topology topo,
86-
int order) {
87-
auto ip = getIntegration(topo)->getAccurate(order)->getPoints();
125+
template <Mesh_Topology topo> auto getIntegrationPoints(int order) {
126+
auto ip = getIntegration<topo>()->getAccurate(order)->getPoints();
88127
return ip;
89128
}
90129

91130
template <typename FieldElement>
92-
Kokkos::View<MeshField::Real **>
93-
getIntegrationPointLocalCoords(FieldElement &fes,
94-
std::vector<IntegrationPoint> ip) {
131+
Kokkos::View<MeshField::Real **> getIntegrationPointLocalCoords(
132+
FieldElement &fes,
133+
std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
95134
const auto numPtsPerElm = ip.size();
96135
const auto numMeshEnts = fes.numMeshEnts;
97136
const auto meshEntDim = fes.MeshEntDim;
@@ -115,9 +154,9 @@ getIntegrationPointLocalCoords(FieldElement &fes,
115154
}
116155

117156
template <typename FieldElement>
118-
Kokkos::View<Real *>
119-
getIntegrationPointWeights(FieldElement &fes,
120-
std::vector<IntegrationPoint> ip) {
157+
Kokkos::View<Real *> getIntegrationPointWeights(
158+
FieldElement &fes,
159+
std::vector<IntegrationPoint<FieldElement::MeshEntDim + 1>> ip) {
121160
const auto numPtsPerElm = ip.size();
122161
const auto numMeshEnts = fes.numMeshEnts;
123162
const auto meshEntDim = fes.MeshEntDim;
@@ -194,13 +233,9 @@ class Integrator {
194233
* FIXME make the sensible
195234
* */
196235
template <typename FieldElement> void process(FieldElement &fes) {
197-
const auto topo = fes.elm2dof.getTopology();
198-
if (topo[0] != MeshField::Triangle) {
199-
throw std::invalid_argument(
200-
"Integrator::process only supports triangles.");
201-
}
236+
constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology();
202237
pre();
203-
auto ip = getIntegrationPoints(topo[0], order);
238+
auto ip = getIntegrationPoints<topo[0]>(order);
204239
auto localCoords = getIntegrationPointLocalCoords(fes, ip);
205240
auto weights = getIntegrationPointWeights(fes, ip);
206241
auto dV = getJacobianDeterminants(fes, localCoords, ip.size());

‎test/testCountIntegrator.cpp‎

Lines changed: 28 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -14,12 +14,12 @@
1414

1515
using ExecutionSpace = Kokkos::DefaultExecutionSpace;
1616
using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
17-
18-
Omega_h::Mesh createMeshTri18(Omega_h::Library &lib) {
17+
template <size_t dim> Omega_h::Mesh createMesh(Omega_h::Library &lib) {
1918
auto world = lib.world();
2019
const auto family = OMEGA_H_SIMPLEX;
2120
auto len = 1.0;
22-
return Omega_h::build_box(world, family, len, len, 0.0, 3, 3, 0);
21+
const auto numEnts3d = (dim == 3 ? 3 : 0);
22+
return Omega_h::build_box(world, family, len, len, len, 3, 3, numEnts3d);
2323
}
2424

2525
template <typename AnalyticFunction, typename ShapeField>
@@ -57,13 +57,19 @@ class CountIntegrator : public MeshField::Integrator {
5757
}
5858
};
5959

60-
template <template <typename...> typename Controller>
60+
template <template <typename...> typename Controller, size_t dim>
6161
void doRun(Omega_h::Mesh &mesh,
62-
MeshField::OmegahMeshField<ExecutionSpace, 2, Controller> &omf) {
62+
MeshField::OmegahMeshField<ExecutionSpace, dim, Controller> &omf) {
6363
const auto ShapeOrder = 1;
6464
auto field = omf.getCoordField();
65-
const auto [shp, map] =
66-
MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
65+
auto shapeSet = [&]() {
66+
if constexpr (dim == 3) {
67+
return MeshField::Omegah::getTetrahedronElement<ShapeOrder>(mesh);
68+
} else {
69+
return MeshField::Omegah::getTriangleElement<ShapeOrder>(mesh);
70+
}
71+
};
72+
const auto [shp, map] = shapeSet();
6773
MeshField::FieldElement fes(mesh.nelems(), field, shp, map);
6874

6975
CountIntegrator countInt(fes);
@@ -74,18 +80,28 @@ void doRun(Omega_h::Mesh &mesh,
7480
int main(int argc, char **argv) {
7581
Kokkos::initialize(argc, argv);
7682
auto lib = Omega_h::Library(&argc, &argv);
77-
auto mesh = createMeshTri18(lib);
7883
#ifdef MESHFIELDS_ENABLE_CABANA
7984
{
85+
auto mesh2D = createMesh<2>(lib);
86+
auto mesh3D = createMesh<3>(lib);
8087
MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::CabanaController>
81-
omf(mesh);
82-
doRun<MeshField::CabanaController>(mesh, omf);
88+
omf2D(mesh2D);
89+
doRun<MeshField::CabanaController>(mesh2D, omf2D);
90+
MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::CabanaController>
91+
omf3D(mesh3D);
92+
doRun<MeshField::CabanaController>(mesh3D, omf3D);
8393
}
8494
#endif
8595
{
96+
auto mesh2D = createMesh<2>(lib);
97+
auto mesh3D = createMesh<3>(lib);
8698
MeshField::OmegahMeshField<ExecutionSpace, 2, MeshField::KokkosController>
87-
omf(mesh);
88-
doRun<MeshField::KokkosController>(mesh, omf);
99+
omf2D(mesh2D);
100+
doRun<MeshField::KokkosController>(mesh2D, omf2D);
101+
MeshField::OmegahMeshField<ExecutionSpace, 3, MeshField::KokkosController>
102+
omf3D(mesh3D);
103+
doRun<MeshField::KokkosController>(mesh3D, omf3D);
89104
}
105+
Kokkos::finalize();
90106
return 0;
91107
}

‎test/testElementJacobian2d.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ void triJacobian() {
8484
LinearTriangleToVertexField());
8585

8686
Kokkos::View<MeshField::Real *[3]> lc("localCoords", 1);
87-
Kokkos::deep_copy(lc, 1.0 / 2);
87+
Kokkos::deep_copy(lc, 1.0 / 3);
8888
const auto numPtsPerElement = 1;
8989
const auto J = MeshField::getJacobians(f, lc, numPtsPerElement);
9090
const auto J_h =

0 commit comments

Comments
 (0)