Skip to content

Commit 0db40d2

Browse files
jacobmersonjacobmerson
authored
Merge pull request #368 from Sichao25/evaluation-library
Use meshFields parametric coordinates for dof coordinate construction
2 parents 4c92c7d + 69f3ef1 commit 0db40d2

2 files changed

Lines changed: 216 additions & 66 deletions

File tree

‎src/pcms/field/layout/mesh_fields.cpp‎

Lines changed: 118 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "pcms/utility/inclusive_scan.h"
44
#include "pcms/utility/profile.h"
55
#include <Omega_h_for.hpp>
6+
#include <MeshField.hpp>
67
#include <memory>
78

89
namespace pcms
@@ -36,59 +37,133 @@ Omega_h::Write<Omega_h::GO> GetGidsHelper(LO total_ents,
3637
return owned_gids;
3738
}
3839

39-
// this is a workaround to specify the parametric coordinates for MeshFields to
40-
// be replaced when https://github.com/SCOREC/meshFields/issues/70 is resolved
41-
struct ComputeVertexCoordsFunctor
40+
namespace
4241
{
43-
Kokkos::View<Real**> dof_holder_coords_;
44-
Omega_h::Reals coords_;
45-
size_t offset_;
46-
47-
ComputeVertexCoordsFunctor(Kokkos::View<Real**> dof_holder_coords,
48-
Omega_h::Reals coords, size_t offset)
49-
: dof_holder_coords_(dof_holder_coords), coords_(coords), offset_(offset)
50-
{
42+
// Maps a meshFields topology to the dimension.
43+
KOKKOS_INLINE_FUNCTION int TopologyToDim(MeshField::Mesh_Topology topo)
44+
{
45+
switch (topo) {
46+
case MeshField::Vertex: return 0;
47+
case MeshField::Edge: return 1;
48+
case MeshField::Triangle: return 2;
49+
case MeshField::Tetrahedron: return 3;
50+
default: return -1;
5151
}
52+
}
5253

53-
KOKKOS_INLINE_FUNCTION
54-
void operator()(LO i) const
55-
{
56-
dof_holder_coords_(offset_ + i, 0) = coords_[2 * i + 0];
57-
dof_holder_coords_(offset_ + i, 1) = coords_[2 * i + 1];
54+
// Gets the appropriate MeshField element
55+
template <int Dim, int Order>
56+
auto GetMeshFieldElement(Omega_h::Mesh& mesh)
57+
{
58+
if constexpr (Dim == 2) {
59+
return MeshField::Omegah::getTriangleElement<Order>(mesh);
60+
} else {
61+
return MeshField::Omegah::getTetrahedronElement<Order>(mesh);
5862
}
59-
};
63+
}
6064

61-
// this is a workaround to specify the parametric coordinates for MeshFields to
62-
// be replaced when https://github.com/SCOREC/meshFields/issues/70 is resolved
63-
struct ComputeEdgeCoordsFunctor
65+
template <int Dim, int Order>
66+
void BuildDofHolderCoordsFromMeshFieldImpl(
67+
Omega_h::Mesh& mesh, Kokkos::View<Real**> holder_coords,
68+
const std::array<int, 4>& nodes_per_dim)
6469
{
65-
Kokkos::View<Real**> dof_holder_coords_;
66-
Omega_h::Reals coords_;
67-
Omega_h::LOs edge_verts_;
68-
size_t offset_;
70+
const auto elem = GetMeshFieldElement<Dim, Order>(mesh);
6971

70-
ComputeEdgeCoordsFunctor(Kokkos::View<Real**> dof_holder_coords,
71-
Omega_h::Reals coords, Omega_h::LOs edge_verts,
72-
size_t offset)
73-
: dof_holder_coords_(dof_holder_coords),
74-
coords_(coords),
75-
edge_verts_(edge_verts),
76-
offset_(offset)
72+
using ShapeT = std::decay_t<decltype(elem.shp)>;
73+
constexpr size_t numNodes = ShapeT::numNodes;
74+
constexpr size_t meshDim = ShapeT::meshEntDim;
75+
constexpr MeshField::Mesh_Topology elemTopo =
76+
(Dim == 2) ? MeshField::Triangle : MeshField::Tetrahedron;
77+
const auto map = elem.map;
78+
79+
// Compute the starting offset of DOF holders for each entity dimension.
80+
Kokkos::Array<LO, 4> dof_holder_offset;
7781
{
82+
LO cur_offset = 0;
83+
for (int d = 0; d < 4; ++d) {
84+
dof_holder_offset[d] = cur_offset;
85+
if (d <= mesh.dim() && nodes_per_dim[d] != 0)
86+
cur_offset += mesh.nents(d);
87+
}
7888
}
7989

80-
KOKKOS_INLINE_FUNCTION
81-
void operator()(LO i) const
82-
{
83-
auto verts = Omega_h::gather_verts<2>(edge_verts_, i);
84-
Real x0 = coords_[2 * verts[0] + 0];
85-
Real y0 = coords_[2 * verts[0] + 1];
86-
Real x1 = coords_[2 * verts[1] + 0];
87-
Real y1 = coords_[2 * verts[1] + 1];
88-
dof_holder_coords_(offset_ + i, 0) = (x0 + x1) / 2;
89-
dof_holder_coords_(offset_ + i, 1) = (y0 + y1) / 2;
90+
const LO nelems = mesh.nelems();
91+
const auto vtx_coords = mesh.coords();
92+
const auto param_coords = elem.shp.getNodeParametricCoords();
93+
94+
Kokkos::parallel_for(
95+
"MeshFieldsDofHolderCoords",
96+
Kokkos::RangePolicy<DeviceMemorySpace::execution_space>(0, nelems),
97+
KOKKOS_LAMBDA(const LO e) {
98+
constexpr size_t NVerts = Dim + 1;
99+
100+
Real Xv[NVerts * Dim];
101+
for (size_t b = 0; b < NVerts; ++b) {
102+
const auto hv = map(static_cast<MeshField::LO>(b), 0, e, elemTopo);
103+
const LO v = hv.entity;
104+
for (size_t d = 0; d < Dim; ++d)
105+
Xv[b * Dim + d] = vtx_coords[static_cast<size_t>(v) * Dim + d];
106+
}
107+
108+
for (size_t n = 0; n < numNodes; ++n) {
109+
const auto h = map(static_cast<MeshField::LO>(n), 0, e, elemTopo);
110+
111+
// Compute the barycentric coordinates L[b] at the node's parametric
112+
// coordinates.
113+
Real L[NVerts];
114+
L[0] = 1;
115+
for (size_t d = 0; d < Dim; ++d) {
116+
const Real xi_d = param_coords[n * Dim + d];
117+
L[d + 1] = xi_d;
118+
L[0] -= xi_d;
119+
}
120+
121+
Real X[Dim] = {0};
122+
for (size_t b = 0; b < NVerts; ++b)
123+
for (size_t d = 0; d < Dim; ++d)
124+
X[d] += L[b] * Xv[b * Dim + d];
125+
126+
const int dim_of_holder = TopologyToDim(h.topo);
127+
const LO row = dof_holder_offset[dim_of_holder] + h.entity;
128+
for (size_t d = 0; d < Dim; ++d)
129+
holder_coords(row, d) = X[d];
130+
}
131+
});
132+
}
133+
134+
void BuildDofHolderCoordsFromMeshField(Omega_h::Mesh& mesh,
135+
Kokkos::View<Real**> holder_coords,
136+
const std::array<int, 4>& nodes_per_dim)
137+
{
138+
int dim = mesh.dim();
139+
int order = 0;
140+
for (int i = 0; i <= dim; ++i) {
141+
if (nodes_per_dim[i] == 1)
142+
++order;
143+
else if (nodes_per_dim[i] != 0) {
144+
std::cerr << "Unsupported" << std::endl;
145+
std::abort();
146+
}
90147
}
91-
};
148+
149+
if (dim == 2 && order == 1)
150+
BuildDofHolderCoordsFromMeshFieldImpl<2, 1>(mesh, holder_coords,
151+
nodes_per_dim);
152+
else if (dim == 2 && order == 2)
153+
BuildDofHolderCoordsFromMeshFieldImpl<2, 2>(mesh, holder_coords,
154+
nodes_per_dim);
155+
else if (dim == 3 && order == 1)
156+
BuildDofHolderCoordsFromMeshFieldImpl<3, 1>(mesh, holder_coords,
157+
nodes_per_dim);
158+
else if (dim == 3 && order == 2)
159+
BuildDofHolderCoordsFromMeshFieldImpl<3, 2>(mesh, holder_coords,
160+
nodes_per_dim);
161+
else {
162+
std::cerr << "Unsupported element/order combination" << std::endl;
163+
std::abort();
164+
}
165+
}
166+
} // namespace
92167

93168
struct CopyClassInfoFunctor
94169
{
@@ -152,32 +227,9 @@ MeshFieldsAdapterLayout::MeshFieldsAdapterLayout(
152227
std::abort();
153228
}
154229

155-
auto coords = mesh_.coords();
230+
BuildDofHolderCoordsFromMeshField(mesh_, dof_holder_coords_, nodes_per_dim);
156231

157232
size_t offset = 0;
158-
for (int i = 0; i <= mesh_.dim(); ++i) {
159-
if (nodes_per_dim[i] == 1) {
160-
if (i == 0) {
161-
ComputeVertexCoordsFunctor functor(dof_holder_coords_, coords, offset);
162-
Kokkos::parallel_for(mesh_.nents(0), functor);
163-
} else if (i == 1) {
164-
auto edge_verts = mesh_.ask_verts_of(1);
165-
ComputeEdgeCoordsFunctor functor(dof_holder_coords_, coords, edge_verts,
166-
offset);
167-
Kokkos::parallel_for(mesh_.nents(1), functor);
168-
} else {
169-
std::cerr << "Unsupported" << std::endl;
170-
std::abort();
171-
}
172-
} else if (nodes_per_dim[i] != 0) {
173-
std::cerr << "Unsupported" << std::endl;
174-
std::abort();
175-
}
176-
177-
offset += mesh.nents(i);
178-
}
179-
180-
offset = 0;
181233
for (int i = 0; i <= mesh_.dim(); ++i) {
182234
if (nodes_per_dim_[i]) {
183235
auto ids = mesh_.get_array<Omega_h::ClassId>(i, "class_id");

‎test/test_field_evaluation.cpp‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,9 @@
66
#include "pcms/field/field_metadata.h"
77
#include "pcms/utility/assert.h"
88
#include "field_test_utils.h"
9+
#include <array>
910
#include <cmath>
11+
#include <vector>
1012

1113
using pcms::Real;
1214

@@ -63,6 +65,102 @@ TEST_CASE("evaluate quadratic 2d meshfields_field")
6365
}
6466
#endif
6567

68+
#ifdef PCMS_ENABLE_MESHFIELDS
69+
// Recomputes the expected DOF-holder coordinates from mesh geometry exactly as
70+
// the layout's isoparametric map does for affine simplex elements, accumulated
71+
// per entity dimension in the same block order the layout stores:
72+
// dim 0 (vertex) holders -> the vertex's own coordinate,
73+
// dim 1 (edge) holders -> the midpoint of the edge's two endpoint vertices.
74+
// Returns a flat node-major vector matching GetDOFHolderCoordinates() layout.
75+
static std::vector<Real> ComputeManualHolderCoords(
76+
Omega_h::Mesh& mesh, const std::array<int, 4>& nodes_per_dim)
77+
{
78+
const int dim = mesh.dim();
79+
auto mesh_coords_h = Omega_h::HostRead<Omega_h::Real>(mesh.coords());
80+
auto edge2vtx_h = Omega_h::HostRead<Omega_h::LO>(mesh.ask_verts_of(1));
81+
82+
std::vector<Real> expected;
83+
for (int e_dim = 0; e_dim <= dim; ++e_dim) {
84+
if (nodes_per_dim[static_cast<size_t>(e_dim)] != 1)
85+
continue;
86+
const int num_ents = mesh.nents(e_dim);
87+
for (int ent = 0; ent < num_ents; ++ent) {
88+
if (e_dim == 0) {
89+
for (int d = 0; d < dim; ++d)
90+
expected.push_back(mesh_coords_h[ent * dim + d]);
91+
} else { // e_dim == 1: edge midpoint
92+
const int a = edge2vtx_h[2 * ent + 0];
93+
const int b = edge2vtx_h[2 * ent + 1];
94+
for (int d = 0; d < dim; ++d)
95+
expected.push_back(
96+
0.5 * (mesh_coords_h[a * dim + d] + mesh_coords_h[b * dim + d]));
97+
}
98+
}
99+
}
100+
return expected;
101+
}
102+
103+
// Builds a MeshFields layout of the given order on a small known box and checks
104+
// that GetDOFHolderCoordinates() matches the manual reference computation.
105+
static void CheckMeshFieldsDofHolderCoordsAgainstManual(Omega_h::Library& lib,
106+
int order)
107+
{
108+
auto mesh =
109+
Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, 2, 3, 0, false);
110+
111+
std::array<int, 4> nodes_per_dim{};
112+
if (order == 1)
113+
nodes_per_dim = {1, 0, 0, 0};
114+
else if (order == 2)
115+
nodes_per_dim = {1, 1, 0, 0};
116+
117+
auto factory = pcms::LagrangeFunctionSpace::FromMesh(
118+
mesh, order, 1, pcms::CoordinateSystem::Cartesian, "global",
119+
pcms::LagrangeFunctionSpace::Backend::MeshFields);
120+
auto layout = factory->GetLayout();
121+
122+
// Sanity-check the block layout: vertices then edges.
123+
if (order == 2) {
124+
const auto offsets = layout->GetEntOffsets();
125+
const int nverts = mesh.nents(0);
126+
REQUIRE(offsets[0] == 0);
127+
REQUIRE(offsets[1] == static_cast<size_t>(nverts));
128+
}
129+
130+
const std::vector<Real> expected =
131+
ComputeManualHolderCoords(mesh, nodes_per_dim);
132+
133+
auto coords = layout->GetDOFHolderCoordinates().GetValues();
134+
const int dim = mesh.dim();
135+
REQUIRE(coords.extent(1) == static_cast<size_t>(dim));
136+
REQUIRE(coords.extent(0) == expected.size() / static_cast<size_t>(dim));
137+
138+
auto coords_h = pcms::test::CopyCoordinatesToHost(
139+
coords, static_cast<int>(coords.extent(0)), dim);
140+
141+
size_t idx = 0;
142+
for (size_t r = 0; r < coords_h.extent(0); ++r) {
143+
for (int d = 0; d < dim; ++d) {
144+
INFO("row " << r << " dim " << d);
145+
REQUIRE(coords_h(r, d) == Catch::Approx(expected[idx]));
146+
++idx;
147+
}
148+
}
149+
}
150+
151+
TEST_CASE("MeshFields: linear dof-holder coordinates match manual reference")
152+
{
153+
auto lib = Omega_h::Library{};
154+
CheckMeshFieldsDofHolderCoordsAgainstManual(lib, 1);
155+
}
156+
157+
TEST_CASE("MeshFields: quadratic dof-holder coordinates match manual reference")
158+
{
159+
auto lib = Omega_h::Library{};
160+
CheckMeshFieldsDofHolderCoordsAgainstManual(lib, 2);
161+
}
162+
#endif
163+
66164
TEST_CASE("evaluate quadratic 2d omega_h_field throws")
67165
{
68166
auto lib = Omega_h::Library{};

0 commit comments

Comments
 (0)