Skip to content
Merged
Show file tree
Hide file tree
Changes from 8 commits
Commits
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
4 changes: 2 additions & 2 deletions src/MeshField.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,8 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) {
auto coordField =
MeshField::CreateCoordinateField<ExecutionSpace, Controller>(mesh_info);
auto setCoordField = KOKKOS_LAMBDA(const int &i) {
coordField(0, 0, i, MeshField::Vertex) = coords[i * meshDim];
coordField(0, 1, i, MeshField::Vertex) = coords[i * meshDim + 1];
coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim];
coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1];
};
MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx},
setCoordField, "setCoordField");
Expand Down
14 changes: 10 additions & 4 deletions src/MeshField_Element.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -59,9 +59,15 @@ struct FieldElement {
const ElementDofHolderAccessor elm2dofIn)
: numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn),
elm2dof(elm2dofIn) {}

using ValArray = Kokkos::Array<typename FieldAccessor::BaseType,
ShapeType::numComponentsPerDof>;
template <typename T> struct baseType {
using type = T;
};
template <typename T, size_t N> struct baseType<T[N]> {
using type = typename baseType<T>::type;
};
using ValArray =
Kokkos::Array<typename baseType<typename FieldAccessor::BaseType>::type,
ShapeType::numComponentsPerDof>;
static const size_t NumComponents = ShapeType::numComponentsPerDof;

/**
Expand Down Expand Up @@ -91,7 +97,7 @@ struct FieldElement {
for (int ci = 0; ci < shapeFn.numComponentsPerDof; ++ci) {
auto map = elm2dof(ni, ci, ent, topo);
const auto fval =
field(map.node, map.component, map.entity, map.topo);
field(map.entity, map.node, map.component, map.topo);
c[ci] += fval * shapeValues[ni];
}
}
Expand Down
62 changes: 54 additions & 8 deletions src/MeshField_ShapeField.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -173,9 +173,24 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
if (meshInfo.numVtx <= 0) {
fail("mesh has no vertices\n");
}
using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
using Ctrlr = std::conditional_t<
std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
Controller<ExecutionSpace, MemorySpace, DataType[1][1]>,
Controller<MemorySpace, ExecutionSpace, DataType ***>>;
// 1 dof with 1 component per vtx
Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx});
auto createController = [](const int numComp, auto numVtx) {
if constexpr (std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace,
DataType>>) {
return Ctrlr(numVtx);
} else {
return Ctrlr({/*field 0*/ numVtx, 1, numComp});
}
};
Ctrlr kk_ctrl = createController(1, meshInfo.numVtx);
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
using LA = LinearAccessor<decltype(vtxField)>;
using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;
Expand All @@ -188,11 +203,25 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
if (meshInfo.numEdge <= 0) {
fail("mesh has no edges\n");
}
using Ctrlr =
Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>;
using Ctrlr = std::conditional_t<
std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
Controller<ExecutionSpace, MemorySpace, DataType[1][1], DataType[1][1]>,
Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>>;
// 1 dof with 1 comp per vtx/edge
Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx,
/*field 1*/ 1, 1, meshInfo.numEdge});
auto createController = [](const int numComp, auto numVtx, auto numEdge) {
if constexpr (std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace,
DataType>>) {
return Ctrlr(std::max(numVtx, numEdge));
} else {
return Ctrlr({/*field 0*/ numVtx, 1, numComp,
/*field 1*/ numEdge, 1, numComp});
}
};
Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge);
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
using QA = QuadraticAccessor<decltype(vtxField), decltype(edgeField)>;
Expand Down Expand Up @@ -224,6 +253,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) {
* @param meshInfo defines on-process mesh metadata
* @return a linear ShapeField
*/

template <typename ExecutionSpace, template <typename...> typename Controller =
MeshField::KokkosController>
auto CreateCoordinateField(const MeshInfo &meshInfo) {
Expand All @@ -232,9 +262,25 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) {
}
using DataType = Real;
using MemorySpace = typename ExecutionSpace::memory_space;
using Ctrlr = Controller<MemorySpace, ExecutionSpace, DataType ***>;
const int numComp = meshInfo.dim;
Ctrlr kk_ctrl({/*field 0*/ 1, numComp, meshInfo.numVtx});
// FIXME Oversized cabana datatypes when numComp = 1|2
using Ctrlr = std::conditional_t<
std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
Controller<ExecutionSpace, MemorySpace, DataType[1][3]>,
Controller<MemorySpace, ExecutionSpace, DataType ***>>;
auto createController = [](const int numComp, auto numVtx) {
if constexpr (std::is_same_v<
Controller<ExecutionSpace, MemorySpace, DataType>,
MeshField::CabanaController<ExecutionSpace, MemorySpace,
DataType>>) {
return Ctrlr(numVtx);
} else {
return Ctrlr({/*field 0*/ numVtx, 1, numComp});
}
};
Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx);
auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
using LA = LinearAccessor<decltype(vtxField)>;
using LinearLagrangeShapeField = ShapeField<Ctrlr, LinearTriangleShape, LA>;
Expand Down
4 changes: 2 additions & 2 deletions test/testOmegahCoordField.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,8 @@ bool triangleLocalToGlobal(Omega_h::Mesh mesh) {
MeshField::CreateCoordinateField<ExecutionSpace,
MeshField::KokkosController>(meshInfo);
auto setCoordField = KOKKOS_LAMBDA(const int &i) {
coordField(0, 0, i, MeshField::Vertex) = coords[i * MeshDim];
coordField(0, 1, i, MeshField::Vertex) = coords[i * MeshDim + 1];
coordField(i, 0, 0, MeshField::Vertex) = coords[i * MeshDim];
coordField(i, 0, 1, MeshField::Vertex) = coords[i * MeshDim + 1];
};
MeshField::parallel_for(ExecutionSpace(), {0}, {meshInfo.numVtx},
setCoordField, "setCoordField");
Expand Down
89 changes: 87 additions & 2 deletions test/testOmegahElement.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@
#include "MeshField_Fail.hpp" //remove?
#include "MeshField_For.hpp" //remove?
#include "MeshField_ShapeField.hpp" //remove?
#ifdef MESHFIELDS_ENABLE_CABANA
#include "CabanaController.hpp"
#endif
#include "Omega_h_build.hpp"
#include "Omega_h_file.hpp"
#include "Omega_h_simplex.hpp"
Expand Down Expand Up @@ -85,7 +88,7 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
// - TODO should be encoded in the field?
const auto x = coords[vtx * MeshDim];
const auto y = coords[vtx * MeshDim + 1];
field(0, 0, vtx, MeshField::Vertex) = func(x, y);
field(vtx, 0, 0, MeshField::Vertex) = func(x, y);
};
MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()},
setFieldAtVertices, "setFieldAtVertices");
Expand All @@ -106,7 +109,7 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) {
const auto x = (coords[left * MeshDim] + coords[right * MeshDim]) / 2.0;
const auto y =
(coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0;
field(0, 0, edge, MeshField::Edge) = func(x, y);
field(edge, 0, 0, MeshField::Edge) = func(x, y);
};
MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()},
setFieldAtEdges, "setFieldAtEdges");
Expand Down Expand Up @@ -142,6 +145,88 @@ int main(int argc, char **argv) {
Kokkos::initialize(argc, argv);
auto lib = Omega_h::Library(&argc, &argv);
MeshField::Debug = true;
{
auto mesh = createMeshTri18(lib);
MeshField::OmegahMeshField<ExecutionSpace, MeshField::CabanaController> omf(
mesh);

// setup field with values from the analytic function
static const size_t OnePtPerElem = 1;
static const size_t ThreePtsPerElem = 3;
auto centroids = createElmAreaCoords<OnePtPerElem>(
mesh.nfaces(), {1 / 3.0, 1 / 3.0, 1 / 3.0});
auto interior =
createElmAreaCoords<OnePtPerElem>(mesh.nfaces(), {0.1, 0.4, 0.5});
auto vertex =
createElmAreaCoords<OnePtPerElem>(mesh.nfaces(), {0.0, 0.0, 1.0});
// clang-format off
auto allVertices = createElmAreaCoords<ThreePtsPerElem>(mesh.nfaces(),
{1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0});
const auto cases = {TestCoords{centroids, OnePtPerElem, "centroids"},
TestCoords{interior, OnePtPerElem, "interior"},
TestCoords{vertex, OnePtPerElem, "vertex"},
TestCoords{allVertices, ThreePtsPerElem, "allVertices"}};
// clang-format on

auto coords = mesh.coords();
const auto MeshDim = 2;
for (auto testCase : cases) {
using ViewType = decltype(testCase.coords);
{
const auto ShapeOrder = 1;
auto field =
omf.CreateLagrangeField<MeshField::Real, ShapeOrder, MeshDim>();
auto func = LinearFunction();
setVertices(mesh, func, field);
using FieldType = decltype(field);
auto result =
omf.triangleLocalPointEval<LinearFunction, ViewType, FieldType>(
testCase.coords, testCase.NumPtsPerElem, LinearFunction{},
field);
auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
LinearFunction{});
if (failed)
doFail("linear", "linear", testCase.name);
}
{
const auto ShapeOrder = 2;
auto field =
omf.CreateLagrangeField<MeshField::Real, ShapeOrder, MeshDim>();
auto func = QuadraticFunction();
setVertices(mesh, func, field);
setEdges(mesh, func, field);
using FieldType = decltype(field);
auto result =
omf.triangleLocalPointEval<QuadraticFunction, ViewType, FieldType>(
testCase.coords, testCase.NumPtsPerElem, QuadraticFunction{},
field);
auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
QuadraticFunction{});
if (failed)
doFail("quadratic", "quadratic", testCase.name);
}
{
const auto ShapeOrder = 2;
auto field =
omf.CreateLagrangeField<MeshField::Real, ShapeOrder, MeshDim>();
auto func = LinearFunction();
auto coords = mesh.coords();
setVertices(mesh, func, field);
setEdges(mesh, func, field);
using FieldType = decltype(field);
auto result =
omf.triangleLocalPointEval<LinearFunction, ViewType, FieldType>(
testCase.coords, testCase.NumPtsPerElem, LinearFunction{},
field);
auto failed = checkResult(mesh, result, omf.getCoordField(), testCase,
LinearFunction{});
if (failed)
doFail("quadratic", "linear", testCase.name);
}
}
}
{
auto mesh = createMeshTri18(lib);
MeshField::OmegahMeshField<ExecutionSpace, MeshField::KokkosController> omf(
Expand Down
Loading