Skip to content

Commit b9319c5

Browse files
cwsmithcwsmith
authored
Merge pull request #61 from SCOREC/jk/numCompLagrange
Jk/num comp lagrange
2 parents 7590f4a + 0c5d459 commit b9319c5

File tree

6 files changed

+115
-81
lines changed

6 files changed

+115
-81
lines changed

‎src/MeshField.hpp‎

Lines changed: 7 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -189,11 +189,11 @@ class OmegahMeshField {
189189
static_assert(dim == 1 || dim == 2 || dim == 3);
190190
}
191191

192-
template <typename DataType, size_t order>
192+
template <typename DataType, size_t order, size_t numComp>
193193
// Ordering of field indexing changed to 'entity, node, component'
194194
auto CreateLagrangeField() {
195195
return MeshField::CreateLagrangeField<ExecutionSpace, Controller, DataType,
196-
order, dim>(meshInfo);
196+
order, dim, numComp>(meshInfo);
197197
}
198198

199199
auto getCoordField() { return coordField; }
@@ -219,12 +219,13 @@ class OmegahMeshField {
219219
return offsets;
220220
}
221221

222-
// evaluate a field at the specified local coordinates for each triangle
222+
// evaluate a field at the specified local coordinate for each triangle
223223
template <typename ViewType, typename ShapeField>
224224
auto triangleLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
225225
ShapeField field) {
226226
auto offsets = createOffsets(meshInfo.numTri, NumPtsPerElem);
227-
auto eval = triangleLocalPointEval(localCoords, offsets, field);
227+
auto eval = triangleLocalPointEval<ViewType, ShapeField>(localCoords,
228+
offsets, field);
228229
return eval;
229230
}
230231

@@ -243,7 +244,8 @@ class OmegahMeshField {
243244

244245
const auto [shp, map] = Omegah::getTriangleElement<ShapeOrder>(mesh);
245246

246-
MeshField::FieldElement f(meshInfo.numTri, field, shp, map);
247+
MeshField::FieldElement<ShapeField, decltype(shp), decltype(map)> f(
248+
meshInfo.numTri, field, shp, map);
247249
auto eval = MeshField::evaluate(f, localCoords, offsets);
248250
return eval;
249251
}

‎src/MeshField_Element.hpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -150,8 +150,8 @@ struct FieldElement {
150150
};
151151
using ValArray =
152152
Kokkos::Array<typename baseType<typename FieldAccessor::BaseType>::type,
153-
ShapeType::numComponentsPerDof>;
154-
static const size_t NumComponents = ShapeType::numComponentsPerDof;
153+
FieldAccessor::numComp>;
154+
static const size_t NumComponents = FieldAccessor::numComp;
155155

156156
/**
157157
* @brief
@@ -173,11 +173,11 @@ struct FieldElement {
173173
assert(ent < numMeshEnts);
174174
ValArray c;
175175
const auto shapeValues = shapeFn.getValues(localCoord);
176-
for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci)
176+
for (int ci = 0; ci < NumComponents; ++ci)
177177
c[ci] = 0;
178178
for (auto topo : elm2dof.getTopology()) { // element topology
179179
for (int ni = 0; ni < shapeFn.numNodes; ++ni) {
180-
for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci) {
180+
for (int ci = 0; ci < NumComponents; ++ci) {
181181
auto map = elm2dof(ni, ci, ent, topo);
182182
const auto fval =
183183
field(map.entity, map.node, map.component, map.topo);

‎src/MeshField_Shape.hpp‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,6 @@ using Vector4 = Kokkos::Array<Real, 4>;
3636

3737
struct LinearEdgeShape {
3838
static const size_t numNodes = 2;
39-
static const size_t numComponentsPerDof = 1;
4039
static const size_t meshEntDim = 1;
4140
constexpr static Mesh_Topology DofHolders[1] = {Vertex};
4241
constexpr static size_t Order = 1;
@@ -90,7 +89,6 @@ struct LinearTriangleShape {
9089

9190
struct LinearTriangleCoordinateShape {
9291
static const size_t numNodes = 3;
93-
static const size_t numComponentsPerDof = 2;
9492
static const size_t meshEntDim = 2;
9593
constexpr static Mesh_Topology DofHolders[1] = {Vertex};
9694
constexpr static size_t Order = 1;
@@ -109,7 +107,6 @@ struct LinearTriangleCoordinateShape {
109107

110108
struct QuadraticTriangleShape {
111109
static const size_t numNodes = 6;
112-
static const size_t numComponentsPerDof = 1;
113110
static const size_t meshEntDim = 2;
114111
constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
115112
constexpr static size_t NumDofHolders[2] = {3, 3};
@@ -149,7 +146,6 @@ struct QuadraticTriangleShape {
149146

150147
struct QuadraticTetrahedronShape {
151148
static const size_t numNodes = 10;
152-
static const size_t numComponentsPerDof = 1;
153149
static const size_t meshEntDim = 3;
154150
constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge};
155151
constexpr static size_t NumDofHolders[2] = {4, 6};

‎src/MeshField_ShapeField.hpp‎

Lines changed: 19 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -55,10 +55,12 @@ struct MeshInfo {
5555
* @param meshInfoIn defines on-process mesh metadata
5656
* @param mixins object(s) needed to construct the Accessor
5757
*/
58-
template <typename MeshFieldType, typename Shape, typename... Mixins>
58+
template <size_t numCompIn, typename MeshFieldType, typename Shape,
59+
typename... Mixins>
5960
struct ShapeField : public Mixins... {
6061
MeshFieldType meshField;
6162
Shape shape;
63+
static const size_t numComp = numCompIn;
6264
const MeshInfo meshInfo;
6365
constexpr static auto Order = Shape::Order;
6466
ShapeField(MeshFieldType &meshFieldIn, const MeshInfo &meshInfoIn,
@@ -158,7 +160,7 @@ template <typename VtxAccessor> struct LinearAccessor {
158160
template <typename ExecutionSpace,
159161
template <typename...>
160162
typename Controller = MeshField::KokkosController,
161-
typename DataType, size_t order, size_t dim>
163+
typename DataType, size_t order, size_t dim, size_t numComp>
162164
auto CreateLagrangeField(const MeshInfo &meshInfo) {
163165
static_assert((std::is_same_v<Real4, DataType> == true ||
164166
std::is_same_v<Real8, DataType> == true),
@@ -179,10 +181,10 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
179181
std::is_same_v<
180182
Controller<ExecutionSpace, MemorySpace, DataType>,
181183
MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
182-
Controller<ExecutionSpace, MemorySpace, DataType[1][1]>,
184+
Controller<ExecutionSpace, MemorySpace, DataType[1][numComp]>,
183185
Controller<MemorySpace, ExecutionSpace, DataType ***>>;
184186
// 1 dof with 1 component per vtx
185-
auto createController = [](const int numComp, auto numVtx) {
187+
auto createController = [](auto numVtx) {
186188
if constexpr (std::is_same_v<
187189
Controller<ExecutionSpace, MemorySpace, DataType>,
188190
MeshField::CabanaController<ExecutionSpace, MemorySpace,
@@ -192,14 +194,15 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
192194
return Ctrlr({/*field 0*/ numVtx, 1, numComp});
193195
}
194196
};
195-
Ctrlr kk_ctrl = createController(1, meshInfo.numVtx);
197+
Ctrlr kk_ctrl = createController(meshInfo.numVtx);
196198
#else
197199
using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
198-
Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1});
200+
Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp});
199201
#endif
200202
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
201203
using LA = LinearAccessor<decltype(vtxField)>;
202-
using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;
204+
using LinearLagrangeShapeField =
205+
ShapeField<numComp, Ctrlr, LinearTriangleShape, LA>;
203206
LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField});
204207
return llsf;
205208
} else if constexpr (order == 2 && (dim == 2 || dim == 3)) {
@@ -214,10 +217,11 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
214217
std::is_same_v<
215218
Controller<ExecutionSpace, MemorySpace, DataType>,
216219
MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
217-
Controller<ExecutionSpace, MemorySpace, DataType[1][1], DataType[1][1]>,
220+
Controller<ExecutionSpace, MemorySpace, DataType[1][numComp],
221+
DataType[1][numComp]>,
218222
Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>>;
219223
// 1 dof with 1 comp per vtx/edge
220-
auto createController = [](const int numComp, auto numVtx, auto numEdge) {
224+
auto createController = [](auto numVtx, auto numEdge) {
221225
if constexpr (std::is_same_v<
222226
Controller<ExecutionSpace, MemorySpace, DataType>,
223227
MeshField::CabanaController<ExecutionSpace, MemorySpace,
@@ -228,18 +232,18 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
228232
/*field 1*/ numEdge, 1, numComp});
229233
}
230234
};
231-
Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge);
235+
Ctrlr kk_ctrl = createController(meshInfo.numVtx, meshInfo.numEdge);
232236
#else
233237
using Ctrlr =
234238
Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>;
235-
Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1,
236-
/*field 1*/ meshInfo.numEdge, 1, 1});
239+
Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp,
240+
/*field 1*/ meshInfo.numEdge, 1, numComp});
237241
#endif
238242
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
239243
auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
240244
using QA = QuadraticAccessor<decltype(vtxField), decltype(edgeField)>;
241245
using QuadraticLagrangeShapeField =
242-
ShapeField<Ctrlr, QuadraticTriangleShape, QA>;
246+
ShapeField<numComp, Ctrlr, QuadraticTriangleShape, QA>;
243247
QuadraticLagrangeShapeField qlsf(kk_ctrl, meshInfo, {vtxField, edgeField});
244248
return qlsf;
245249
} else {
@@ -302,7 +306,8 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) {
302306
#endif
303307
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
304308
using LA = LinearAccessor<decltype(vtxField)>;
305-
using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;
309+
using LinearLagrangeShapeField =
310+
ShapeField<dim, Ctrlr, LinearTriangleShape, LA>;
306311
LinearLagrangeShapeField llsf(kk_ctrl, meshInfo, {vtxField});
307312
return llsf;
308313
};

‎test/testElement.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ void triangleLocalPointEval() {
4141
const auto numElms = 3; // provided by the mesh
4242
const MeshField::MeshInfo meshInfo{.numVtx = 5, .numTri = 3};
4343
auto field = MeshField::CreateLagrangeField<
44-
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 2>(
44+
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 2, 1>(
4545
meshInfo);
4646

4747
MeshField::FieldElement f(numElms, field, MeshField::LinearTriangleShape(),
@@ -85,7 +85,7 @@ struct LinearEdgeToVertexField {
8585
void edgeLocalPointEval() {
8686
const MeshField::MeshInfo meshInfo{.numVtx = 5, .numEdge = 7, .dim = 1};
8787
auto field = MeshField::CreateLagrangeField<
88-
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 1>(
88+
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 1, 1, 1>(
8989
meshInfo);
9090

9191
MeshField::FieldElement f(meshInfo.numEdge, field,
@@ -137,7 +137,7 @@ void quadraticTriangleLocalPointEval() {
137137
const MeshField::MeshInfo meshInfo{
138138
.numVtx = 3, .numEdge = 3, .numTri = 1, .dim = 2};
139139
auto field = MeshField::CreateLagrangeField<
140-
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 2, 2>(
140+
ExecutionSpace, MeshField::KokkosController, MeshField::Real, 2, 2, 1>(
141141
meshInfo);
142142

143143
MeshField::FieldElement f(meshInfo.numTri, field,
@@ -192,7 +192,7 @@ void quadraticTetrahedronLocalPointEval() {
192192
auto field =
193193
MeshField::CreateLagrangeField<ExecutionSpace,
194194
MeshField::KokkosController,
195-
MeshField::Real, ShapeOrder, MeshDim>(
195+
MeshField::Real, ShapeOrder, MeshDim, 1>(
196196
meshInfo);
197197

198198
MeshField::FieldElement f(meshInfo.numTet, field,

0 commit comments

Comments
 (0)