Templated tests

CMakeLists.txt

@@ -26,6 +26,8 @@ set(MESHFIELD_HEADERS
   src/MeshField_Macros.hpp
   src/KokkosController.hpp
   src/MeshField_Element.hpp
+  src/MeshField_LinearAlgebra.hpp
+  src/MeshField_MappedElement.hpp
   src/MeshField_Integrate.hpp
   src/MeshField_Shape.hpp
   src/MeshField_ShapeField.hpp
@@ -201,40 +203,30 @@ function(meshfields_add_exe EXE_NAME EXE_SRC)
   endif()
 endfunction()
 
-meshfields_add_smoke_exe(KokkosTests test/testKokkos.cpp)
+meshfields_add_smoke_exe(ControllerTests test/testController.cpp)
 meshfields_add_exe(SerializationTests test/testSerialize.cpp)
-meshfields_add_exe(ElementTests test/testElement.cpp)
 meshfields_add_exe(ControllerLifetimeTest test/testControllerLifetime.cpp)
-meshfields_add_exe(ElementJacobian1d test/testElementJacobian1d.cpp)
-meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
-meshfields_add_exe(ElementJacobian3d test/testElementJacobian3d.cpp)
-meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
-meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp)
+meshfields_add_exe(ElementJacobian test/testElementJacobian.cpp)
+meshfields_add_exe(Reproduction test/testReproduction.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
-meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp)
-meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp)
+meshfields_add_exe(Integrator test/testIntegrator.cpp)
+meshfields_add_exe(Interpolate test/testInterpolate.cpp)
+meshfields_add_exe(OmegahMappedElement test/testOmegahMappedElement.cpp)
 meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
-  meshfields_add_smoke_exe(CabanaTests test/testCabana.cpp)
-  meshfields_add_exe(CabanaControllerLifetimeTest test/testCabanaControllerLifetime.cpp)
-  smoke_test_func(CabanaTests ./CabanaTests)
   test_func(ControllerPerformance ./ControllerPerformance)
-  test_func(CabanaControllerLifetimeTest ./CabanaControllerLifetimeTest)
 endif()
 
-smoke_test_func(KokkosTests ./KokkosTests)
+smoke_test_func(ControllerTests ./ControllerTests)
 test_func(SerializationTests ./SerializationTests)
-test_func(ElementTests ./ElementTests)
 test_func(ControllerLifetimeTest ./ControllerLifetimeTest)
-test_func(ElementJacobian1d ./ElementJacobian1d)
-test_func(ElementJacobian2d ./ElementJacobian2d)
-test_func(ElementJacobian3d ./ElementJacobian3d)
-test_func(CountIntegrator ./CountIntegrator)
-smoke_test_func(OmegahTriTests ./OmegahTriTests)
-test_func(OmegahTetTest ./OmegahTetTest)
-test_func(IntegratorMeasure ./IntegratorMeasure)
+test_func(ElementJacobian ./ElementJacobian)
+test_func(Reproduction ./Reproduction)
+test_func(Integrator ./Integrator)
+test_func(Interpolate ./Interpolate)
+test_func(OmegahMappedElement ./OmegahMappedElement)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

src/MeshField.hpp

@@ -5,6 +5,7 @@
 #include "MeshField_Element.hpp"
 #include "MeshField_Fail.hpp"
 #include "MeshField_For.hpp"
+#include "MeshField_MappedElement.hpp"
 #include "MeshField_ShapeField.hpp"
 #include "Omega_h_file.hpp"    //move
 #include "Omega_h_mesh.hpp"    //move
@@ -31,30 +32,27 @@ MeshField::MeshInfo getMeshInfo(Omega_h::Mesh &mesh) {
   }
   return meshInfo;
 }
-
-template <typename ExecutionSpace, size_t dim,
-          template <typename...>
-          typename Controller = MeshField::KokkosController>
-decltype(MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(
-    MeshField::MeshInfo()))
-createCoordinateField(const MeshField::MeshInfo &mesh_info,
-                      Omega_h::Reals coords) {
-  const auto meshDim = mesh_info.dim;
-  auto coordFieldWithCtrlr =
-      MeshField::CreateCoordinateField<ExecutionSpace, Controller, dim>(
-          mesh_info);
-  auto coordField = coordFieldWithCtrlr.field;
-  auto setCoordField = KOKKOS_LAMBDA(const int &i) {
-    coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim];
-    coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1];
-    if constexpr (dim == 3) {
-      coordField(i, 0, 2, MeshField::Vertex) = coords[i * meshDim + 2];
-    }
-  };
-  MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx},
-                          setCoordField, "setCoordField");
-  return coordFieldWithCtrlr;
-}
+// we use functors because cannot use auto return type with KOKKOS_LAMBDAS
+template <int dim, typename Field, typename Coords, typename CoordinateMap>
+struct SetCoordFieldFunctor
+{
+  SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim,
+                       const CoordinateMap& map) :
+    field_(field), coords_(coords), meshDim_(meshDim), map_(map) {}
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const size_t vtx) const {
+    Kokkos::Array<MeshField::Real, dim> x{};
+    for (size_t d = 0; d < meshDim_; ++d)
+      x[d] = coords_[vtx * meshDim_ + d];
+    const auto mapped = map_(x);
+    for (size_t d = 0; d < meshDim_; ++d)
+      field_(vtx, 0, d, MeshField::Vertex) = mapped[d];
+  }
+  Field field_;
+  Coords coords_;
+  int meshDim_;
+  const CoordinateMap map_;
+};
 
 } // anonymous namespace
 
@@ -284,6 +282,16 @@ template <int ShapeOrder> auto getTetrahedronElement(Omega_h::Mesh &mesh) {
 
 } // namespace Omegah
 
+/**
+ * @brief the CoordinateMap that leaves the mesh geometry unchanged
+ */
+struct IdentityCoordinateMap {
+  template <typename Point>
+  KOKKOS_INLINE_FUNCTION Point operator()(Point const &x) const {
+    return x;
+  }
+};
+
 template <typename ExecutionSpace, size_t dim,
           template <typename...> typename Controller =
               MeshField::KokkosController>
@@ -292,15 +300,14 @@ class OmegahMeshField {
   Omega_h::Mesh &mesh;
   const MeshField::MeshInfo meshInfo;
   using CoordField =
-      decltype(createCoordinateField<ExecutionSpace, dim, Controller>(
-          MeshField::MeshInfo(), Omega_h::Reals()));
+      decltype(MeshField::CreateCoordinateField<ExecutionSpace, Controller,
+                                                dim>(MeshField::MeshInfo()));
   CoordField coordField;
 
 public:
   OmegahMeshField(Omega_h::Mesh &mesh_in)
       : mesh(mesh_in), meshInfo(getMeshInfo(mesh)),
-        coordField(createCoordinateField<ExecutionSpace, dim, Controller>(
-            getMeshInfo(mesh_in), mesh_in.coords())) {
+        coordField(CreateLagrangeCoordinateField<1>()) {
     static_assert(dim == 1 || dim == 2 || dim == 3);
   }
 
@@ -311,8 +318,81 @@ class OmegahMeshField {
                                           order, dim, numComp>(meshInfo);
   }
 
+  /**
+   * @brief create a Lagrange coordinate field of the given order
+   *
+   * @details
+   * Each node is placed at map(x), where x is the node's position on the
+   * linear (straight-sided) mesh geometry.  The linear field is set from the
+   * mesh coordinates and higher order fields are interpolated from the
+   * linear coordinate field, so with the default IdentityCoordinateMap the
+   * geometry is the mesh's own and is straight-sided.  A nonlinear map gives
+   * curved geometry for orders above one.
+   *
+   * @tparam order the order of the coordinate field
+   * @param map a device callable taking and returning a
+   * Kokkos::Array<Real, dim> of physical coordinates
+   */
+  template <size_t order, typename CoordinateMap = IdentityCoordinateMap>
+  auto CreateLagrangeCoordinateField(
+      const CoordinateMap &map = CoordinateMap()) const {
+    auto fieldWithCtrlr = CreateLagrangeField<Real, order, dim>();
+    auto field = fieldWithCtrlr.field;
+    // note order 1 is called in the constructor to initialize the coordinate field
+    // it must be called before higher order coordinate fields are constructed
+    if constexpr (order == 1) {
+      const auto meshDim = meshInfo.dim;
+      const auto coords = mesh.coords();
+      MeshField::parallel_for(
+          ExecutionSpace(), {0}, {meshInfo.numVtx},
+          SetCoordFieldFunctor<dim, decltype(field), decltype(coords),
+                               CoordinateMap>(field, coords, meshDim, map),
+          "setCoordField");
+    } else {
+      const auto linear = CreateLagrangeElement(coordField.field);
+      auto target = CreateLagrangeElement(field, linear);
+      MeshField::setDofByInterpolation(target, map);
+    }
+    return fieldWithCtrlr;
+  }
+
   auto getCoordField() { return coordField; }
 
+  /**
+   * @brief create a FieldElement over the mesh elements for a Lagrange field
+   *
+   * @details
+   * The shape and dof map are chosen from the field's order.
+   *
+   * @param coordinate field used to create the element
+   */
+  template <typename Field>
+  auto CreateLagrangeElement(const Field &field) const {
+    static_assert(dim == 2 || dim == 3,
+                  "CreateLagrangeElement supports triangles and tetrahedra");
+    constexpr int order = Field::Order;
+    if constexpr (dim == 2) {
+      const auto [shp, map] = Omegah::getTriangleElement<order>(mesh);
+      return MeshField::FieldElement(meshInfo.numTri, field, shp, map);
+    } else {
+      const auto [shp, map] = Omegah::getTetrahedronElement<order>(mesh);
+      return MeshField::FieldElement(meshInfo.numTet, field, shp, map);
+    }
+  }
+
+  /**
+   * @brief create a MappedElement for a field on the given geometry
+   *
+   * @param field the Lagrange field to evaluate
+   * @param geomElem the geometry element: CreateLagrangeElement of a
+   * coordinate field
+   */
+  template <typename ShapeField, typename GeomElem>
+  auto CreateLagrangeElement(const ShapeField &field,
+                             const GeomElem &geomElem) const {
+    return MeshField::MappedElement(CreateLagrangeElement(field), geomElem);
+  }
+
   // FIXME support 2d and 3d and fields with order>1
   template <typename Field> void writeVtk(Field &field) const {
     using FieldDataType = typename decltype(field.vtxField)::BaseType;