diff --git a/src/CabanaController.hpp b/src/CabanaController.hpp index 7fd3716..2f0b84a 100644 --- a/src/CabanaController.hpp +++ b/src/CabanaController.hpp @@ -94,8 +94,6 @@ class CabanaController { private: // all the type defenitions that are needed us to get the type of the slice // returned by the underlying AoSoA - using soa_t = Cabana::SoA; - template using member_data_t = typename Cabana::MemberTypeAtIndex::type; @@ -112,6 +110,15 @@ class CabanaController { template using wrapper_slice_t = CabanaSliceWrapper, T>; + template + static auto construct(std::integer_sequence, + std::vector &obj) { + return std::tuple< + Cabana::AoSoA, MemorySpace, vecLen>...>{ + Cabana::AoSoA, MemorySpace, vecLen>( + "sliceAoSoA", obj[I])...}; + } + template void construct_sizes() { // There are no dynamic ranks w/ cabana controller. // So we only need to know extent. @@ -134,26 +141,34 @@ class CabanaController { } break; } - extent_sizes[theta][0] = num_tuples; + extent_sizes[theta][0] = num_tuples[this->theta]; this->theta += 1; if constexpr (sizeof...(Tx) != 0) construct_sizes(); } // member vaiables - Cabana::AoSoA aosoa; - const int num_tuples; + std::tuple, MemorySpace, vecLen>...> + aosoa; + int num_tuples[sizeof...(Ts)]; unsigned short theta = 0; int extent_sizes[sizeof...(Ts)][MAX_RANK]; public: - CabanaController() : num_tuples(0) { + CabanaController() { static_assert(sizeof...(Ts) != 0); + for (int i = 0; i < sizeof...(Ts); ++i) { + num_tuples[i] = 0; + } construct_sizes(); } - - CabanaController(int n) : aosoa("sliceAoSoA", n), num_tuples(n) { + CabanaController(const std::initializer_list items) { static_assert(sizeof...(Ts) != 0); + std::vector obj(items); + for (std::size_t i = 0; i < obj.size(); ++i) { + num_tuples[i] = obj[i]; + } + aosoa = construct(std::make_integer_sequence{}, obj); construct_sizes(); } @@ -162,18 +177,19 @@ class CabanaController { assert(j >= 0); assert(j <= MAX_RANK); if (j == 0) - return this->tuples(); + return this->tuples(i); else return extent_sizes[i][j]; } - int tuples() const { return num_tuples; } + int tuples(int i) const { return num_tuples[i]; } template auto makeSlice() { // Creates wrapper object w/ meta data about the slice. using type = std::tuple_element_t; - const int stride = sizeof(soa_t) / sizeof(member_value_t); - auto slice = Cabana::slice(aosoa); + const int stride = sizeof(Cabana::SoA, vecLen>) / + sizeof(member_value_t); + auto slice = Cabana::slice<0>(std::get(aosoa)); int sizes[MAX_RANK]; for (int i = 0; i < MAX_RANK; i++) sizes[i] = this->size(index, i); diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 4c00df8..595457e 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -41,8 +41,8 @@ createCoordinateField(MeshField::MeshInfo mesh_info, Omega_h::Reals coords) { auto coordField = MeshField::CreateCoordinateField(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"); @@ -180,6 +180,7 @@ class OmegahMeshField { getMeshInfo(mesh_in), mesh_in.coords())) {} template + // Ordering of field indexing changed to 'entity, node, component' auto CreateLagrangeField() { return MeshField::CreateLagrangeField(meshInfo); diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index bde61f6..ed3c231 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -59,9 +59,20 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} - - using ValArray = Kokkos::Array; + /* general template for baseType which simply sets type + */ + template struct baseType { + using type = T; + }; + /* template specialization to recursively strip type to get base type + * Example: int[5][6] => int[6] => int + */ + template struct baseType { + using type = typename baseType::type; + }; + using ValArray = + Kokkos::Array::type, + ShapeType::numComponentsPerDof>; static const size_t NumComponents = ShapeType::numComponentsPerDof; /** @@ -91,7 +102,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]; } } diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index ca3a448..0e6541f 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -1,10 +1,10 @@ #ifndef MESHFIELD_SHAPEFIELD_HPP #define MESHFIELD_SHAPEFIELD_HPP -#include "KokkosController.hpp" #ifdef MESHFIELDS_ENABLE_CABANA #include "CabanaController.hpp" #endif +#include "KokkosController.hpp" #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype @@ -173,9 +173,29 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } - using Ctrlr = Controller; +#ifdef MESHFIELDS_ENABLE_CABANA + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; // 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, + MeshField::CabanaController>) { + return Ctrlr({numVtx}); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp}); + } + }; + Ctrlr kk_ctrl = createController(1, meshInfo.numVtx); +#else + using Ctrlr = Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; @@ -188,11 +208,32 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { if (meshInfo.numEdge <= 0) { fail("mesh has no edges\n"); } +#ifdef MESHFIELDS_ENABLE_CABANA + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; + // 1 dof with 1 comp per vtx/edge + auto createController = [](const int numComp, auto numVtx, auto numEdge) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr({numVtx, numEdge}); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp, + /*field 1*/ numEdge, 1, numComp}); + } + }; + Ctrlr kk_ctrl = createController(1, meshInfo.numVtx, meshInfo.numEdge); +#else using Ctrlr = Controller; - // 1 dof with 1 comp per vtx/edge - Ctrlr kk_ctrl({/*field 0*/ 1, 1, meshInfo.numVtx, - /*field 1*/ 1, 1, meshInfo.numEdge}); + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, 1, + /*field 1*/ meshInfo.numEdge, 1, 1}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); auto edgeField = MeshField::makeField(kk_ctrl); using QA = QuadraticAccessor; @@ -224,6 +265,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { * @param meshInfo defines on-process mesh metadata * @return a linear ShapeField */ + template typename Controller = MeshField::KokkosController> auto CreateCoordinateField(const MeshInfo &meshInfo) { @@ -232,9 +274,30 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { } using DataType = Real; using MemorySpace = typename ExecutionSpace::memory_space; - using Ctrlr = Controller; const int numComp = meshInfo.dim; - Ctrlr kk_ctrl({/*field 0*/ 1, numComp, meshInfo.numVtx}); +// FIXME Oversized cabana datatypes when numComp = 1|2 +#ifdef MESHFIELDS_ENABLE_CABANA + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; + auto createController = [](const int numComp, auto numVtx) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr({numVtx}); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp}); + } + }; + Ctrlr kk_ctrl = createController(numComp, meshInfo.numVtx); +#else + using Ctrlr = Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp}); +#endif auto vtxField = MeshField::makeField(kk_ctrl); using LA = LinearAccessor; using LinearLagrangeShapeField = ShapeField; diff --git a/test/testCabana.cpp b/test/testCabana.cpp index 02988ba..815a1b4 100644 --- a/test/testCabana.cpp +++ b/test/testCabana.cpp @@ -36,7 +36,7 @@ using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; void testMakeSliceCabana(int num_tuples) { printf("== START testMakeSliceCabana ==\n"); using Ctrl = MeshField::CabanaController; - Ctrl c(num_tuples); + Ctrl c({num_tuples}); auto field0 = MeshField::makeField(c); @@ -54,7 +54,7 @@ void testParallelReduceCabana() { printf("== START testParallelReduceCabana ==\n"); using Ctrl = MeshField::CabanaController; const int N = 9; - Ctrl c(N); + Ctrl c({N}); { double result = 0, verify = 0; @@ -132,7 +132,7 @@ void testCabanaControllerSize() { using simple = MeshField::CabanaController; - simple c1(a); + simple c1({a}); for (int i = 0; i < 2; i++) { assert(c1.size(0, i) == psi[i]); } @@ -140,7 +140,7 @@ void testCabanaControllerSize() { using multi = MeshField::CabanaController; - multi c2(a); + multi c2({a, a, a}); for (int i = 0; i < 3; i++) { for (int j = 0; j < 4; j++) { assert(c2.size(i, j) == psi[j]); @@ -150,7 +150,7 @@ void testCabanaControllerSize() { using varied = MeshField::CabanaController; - varied c3(a); + varied c3({a, a, a, a}); for (int i = 0; i < 3; i++) assert(c3.size(0, i) == psi[i]); assert(c3.size(1, 0) == psi[0]); @@ -164,6 +164,29 @@ void testCabanaControllerSize() { for (int i = 0; i < 4; i++) { assert(c4.size(0, i) == 0); } + + multi diffc2({a, b, c}); + assert(diffc2.size(0, 0) == psi[0]); + assert(diffc2.size(1, 0) == psi[1]); + assert(diffc2.size(2, 0) == psi[2]); + for (int i = 0; i < 3; i++) { + for (int j = 1; j < 4; j++) { + assert(c2.size(i, j) == psi[j]); + } + } + + varied diffc3({a, b, c, d}); + assert(diffc3.size(0, 0) == psi[0]); + assert(diffc3.size(1, 0) == psi[1]); + assert(diffc3.size(2, 0) == psi[2]); + assert(diffc3.size(3, 0) == psi[3]); + for (int i = 1; i < 3; i++) + assert(diffc3.size(0, i) == psi[i]); + for (int i = 1; i < 4; i++) + assert(diffc3.size(2, i) == psi[i]); + for (int i = 1; i < 2; i++) + assert(diffc3.size(3, i) == psi[i]); + printf("== END testCabanaControllerSize ==\n"); } @@ -182,9 +205,9 @@ void testCabanaFieldSize() { MeshField::CabanaController; using empty = MeshField::CabanaController; - simple c1(a); - multi c2(a); - varied c3(a); + simple c1({a}); + multi c2({a, a, a}); + varied c3({a, a, a, a}); empty c4; const int MAX_RANK = 4; @@ -231,6 +254,43 @@ void testCabanaFieldSize() { assert(field0.size(i) == 0); } } + multi diffc2({a, a + 1, a + 2}); + varied diffc3({a, a + 3, a + 4, a + 5}); + { + auto field0 = MeshField::makeField(diffc2); + auto field1 = MeshField::makeField(diffc2); + auto field2 = MeshField::makeField(diffc2); + assert(field0.size(0) == a); + assert(field1.size(0) == a + 1); + assert(field2.size(0) == a + 2); + for (int i = 1; i < MAX_RANK; ++i) { + assert(field0.size(i) == psi[i]); + assert(field1.size(i) == psi[i]); + assert(field2.size(i) == psi[i]); + } + } + { + auto field0 = MeshField::makeField(diffc3); + auto field1 = MeshField::makeField(diffc3); + auto field2 = MeshField::makeField(diffc3); + auto field3 = MeshField::makeField(diffc3); + assert(field0.size(0) == a); + assert(field1.size(0) == a + 3); + assert(field2.size(0) == a + 4); + assert(field3.size(0) == a + 5); + for (int i = 1; i < 3; i++) { + assert(field0.size(i) == psi[i]); + } + for (int i = 1; i < 1; i++) { + assert(field1.size(i) == psi[i]); + } + for (int i = 1; i < 4; i++) { + assert(field2.size(i) == psi[i]); + } + for (int i = 1; i < 2; i++) { + assert(field3.size(i) == psi[i]); + } + } printf("== END testCabanaFieldSize ==\n"); } @@ -241,7 +301,7 @@ void testCabanaParallelFor() { using simd_ctrlr = MeshField::CabanaController; - simd_ctrlr c1(x); + simd_ctrlr c1({x, x, x, x}); auto field0 = MeshField::makeField(c1); auto field1 = MeshField::makeField(c1); auto field2 = MeshField::makeField(c1); @@ -276,6 +336,45 @@ void testCabanaParallelFor() { MeshField::simd_parallel_for(c1, {0, 0, 0, 0}, {x, y, z, a}, vectorKernel4, "simple_loop"); } + { + using simd_ctrlr = + MeshField::CabanaController; + simd_ctrlr c1({x, x + 1, x + 2, x + 3}); + auto field0 = MeshField::makeField(c1); + auto field1 = MeshField::makeField(c1); + auto field2 = MeshField::makeField(c1); + auto field3 = MeshField::makeField(c1); + + auto vectorKernel = KOKKOS_LAMBDA(const int &i) { + field0(i) = i; + assert(field0(i) == i); + }; + MeshField::simd_parallel_for(c1, {0}, {x}, vectorKernel, "simple_loop"); + + auto vectorKernel2 = KOKKOS_LAMBDA(const int &i, const int &j) { + field1(i, j) = i + j; + assert(field1(i, j) == i + j); + }; + MeshField::simd_parallel_for(c1, {0, 0}, {x + 1, y}, vectorKernel2, + "simple_loop"); + + auto vectorKernel3 = + KOKKOS_LAMBDA(const int &i, const int &j, const int &k) { + field2(i, j, k) = i + j + k; + assert(field2(i, j, k) == i + j + k); + }; + MeshField::simd_parallel_for(c1, {0, 0, 0}, {x + 2, y, z}, vectorKernel3, + "simple_loop"); + + auto vectorKernel4 = + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + field3(i, j, k, l) = i + j + k + l; + assert(field3(i, j, k, l) == i + j + k + l); + }; + MeshField::simd_parallel_for(c1, {0, 0, 0, 0}, {x + 3, y, z, a}, + vectorKernel4, "simple_loop"); + } printf("== END testCabanaParallelFor() ==\n"); } @@ -285,7 +384,7 @@ void testParallelScan() { printf("== START testParallelScan ==\n"); using s_cab = MeshField::CabanaController; - s_cab c1(N); + s_cab c1({N, N}); auto pre = MeshField::makeField(c1); auto post = MeshField::makeField(c1); @@ -312,43 +411,91 @@ void testSetField() { const int N = 10; using cab1 = MeshField::CabanaController; - cab1 c1(N); - auto f1 = MeshField::makeField(c1); - auto f2 = MeshField::makeField(c1); - auto f3 = MeshField::makeField(c1); - auto f4 = MeshField::makeField(c1); - - Kokkos::View v1("1", N); - Kokkos::View v2("2", N, N); - Kokkos::View v3("3", N, N, N); - Kokkos::View v4("4", N, N, N, N); - - Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); - Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N, N, N, N}); - Kokkos::MDRangePolicy> p(start, end); - - Kokkos::parallel_for( - "", p, - KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { - v1(i) += i; - v2(i, j) += i + j; - v3(i, j, k) += i + j + k; - v4(i, j, k, l) += i + j + k + l; - }); - - f1.set(v1); - f2.set(v2); - f3.set(v3); - f4.set(v4); - - Kokkos::parallel_for( - "", p, - KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { - assert(f1(i) == v1(i)); - assert(f2(i, j) == v2(i, j)); - assert(f3(i, j, k) == v3(i, j, k)); - assert(f4(i, j, k, l) == v4(i, j, k, l)); - }); + { + cab1 c1({N, N, N, N}); + auto f1 = MeshField::makeField(c1); + auto f2 = MeshField::makeField(c1); + auto f3 = MeshField::makeField(c1); + auto f4 = MeshField::makeField(c1); + + Kokkos::View v1("1", N); + Kokkos::View v2("2", N, N); + Kokkos::View v3("3", N, N, N); + Kokkos::View v4("4", N, N, N, N); + + Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); + Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N, N, N, N}); + Kokkos::MDRangePolicy> p(start, end); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + v1(i) += i; + v2(i, j) += i + j; + v3(i, j, k) += i + j + k; + v4(i, j, k, l) += i + j + k + l; + }); + + f1.set(v1); + f2.set(v2); + f3.set(v3); + f4.set(v4); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + assert(f1(i) == v1(i)); + assert(f2(i, j) == v2(i, j)); + assert(f3(i, j, k) == v3(i, j, k)); + assert(f4(i, j, k, l) == v4(i, j, k, l)); + }); + } + { + cab1 diffc1({N, N + 1, N + 2, N + 3}); + auto f1 = MeshField::makeField(diffc1); + auto f2 = MeshField::makeField(diffc1); + auto f3 = MeshField::makeField(diffc1); + auto f4 = MeshField::makeField(diffc1); + + Kokkos::View v1("1", N); + Kokkos::View v2("2", N + 1, N); + Kokkos::View v3("3", N + 2, N, N); + Kokkos::View v4("4", N + 3, N, N, N); + + Kokkos::Array start = MeshFieldUtil::to_kokkos_array<4>({0, 0, 0, 0}); + Kokkos::Array end = MeshFieldUtil::to_kokkos_array<4>({N + 3, N, N, N}); + Kokkos::MDRangePolicy> p(start, end); + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + if (i < N) + v1(i) += i; + if (i < N + 1) + v2(i, j) += i + j; + if (i < N + 2) + v3(i, j, k) += i + j + k; + if (i < N + 3) + v4(i, j, k, l) += i + j + k + l; + }); + + f1.set(v1); + f2.set(v2); + f3.set(v3); + f4.set(v4); + + Kokkos::parallel_for( + "", p, + KOKKOS_LAMBDA(const int &i, const int &j, const int &k, const int &l) { + if (i < N) + assert(f1(i) == v1(i)); + if (i < N + 1) + assert(f2(i, j) == v2(i, j)); + if (i < N + 2) + assert(f3(i, j, k) == v3(i, j, k)); + if (i < N + 3) + assert(f4(i, j, k, l) == v4(i, j, k, l)); + }); + } printf("== END testSetField ==\n"); } diff --git a/test/testOmegahCoordField.cpp b/test/testOmegahCoordField.cpp index 4a8fb6d..79997bf 100644 --- a/test/testOmegahCoordField.cpp +++ b/test/testOmegahCoordField.cpp @@ -79,8 +79,8 @@ bool triangleLocalToGlobal(Omega_h::Mesh mesh) { MeshField::CreateCoordinateField(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"); diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index b86bae5..1107151 100644 --- a/test/testOmegahElement.cpp +++ b/test/testOmegahElement.cpp @@ -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" @@ -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"); @@ -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"); @@ -138,25 +141,20 @@ void doFail(std::string_view order, std::string_view function, MeshField::fail(msg); } -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 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( - mesh.nfaces(), {1 / 3.0, 1 / 3.0, 1 / 3.0}); - auto interior = - createElmAreaCoords(mesh.nfaces(), {0.1, 0.4, 0.5}); - auto vertex = - createElmAreaCoords(mesh.nfaces(), {0.0, 0.0, 1.0}); - // clang-format off +template