add reduced quintic WIP

src/pcms/field/CMakeLists.txt

@@ -62,7 +62,9 @@ if (PCMS_ENABLE_MESHFIELDS)
             layout/mesh_fields.h
             evaluator/mesh_fields_backend.h
             evaluator/mesh_fields.h
+            evaluator/mesh_fields_reduced_quintic.h
             data/mesh_fields.h
+            data/mesh_fields_reduced_quintic.h
     )
 endif ()
 

src/pcms/field/data/mesh_fields_reduced_quintic.h

@@ -0,0 +1,105 @@
+#ifndef PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_FIELD_DATA_H
+#define PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_FIELD_DATA_H
+
+#include "pcms/field/layout/mesh_fields.h"
+#include "pcms/field/evaluator/mesh_fields_backend.h"
+#include "pcms/field/field_data.h"
+#include "pcms/field/field_metadata.h"
+#include "pcms/utility/assert.h"
+#include "pcms/utility/arrays.h"
+
+#include <Kokkos_Core.hpp>
+#include <memory>
+
+namespace pcms
+{
+
+// MeshFieldsReducedQuinticFieldData<T> is a near-duplicate of
+// MeshFieldsFieldData<T> that creates a reduced quintic backend instead of
+// a standard Lagrange backend.  The two classes are kept separate to
+// preserve the original API of MeshFieldsFieldData.
+template <typename T>
+class MeshFieldsReducedQuinticFieldData : public FieldData<T>
+{
+public:
+  MeshFieldsReducedQuinticFieldData(
+    std::shared_ptr<const MeshFieldsAdapterLayout> layout,
+    FieldMetadata metadata)
+    : layout_(std::move(layout)),
+      metadata_(metadata),
+      mesh_field_(MakeMeshFieldReducedQuinticBackend<T>(*layout_)),
+      host_data_("meshfields_reduced_quintic_field_data",
+                 static_cast<size_t>(layout_->OwnedSize())),
+      device_data_("meshfields_reduced_quintic_field_data_device",
+                   static_cast<size_t>(layout_->OwnedSize()))
+  {
+    if (!mesh_field_) {
+      throw pcms_error(
+        "MeshFieldsReducedQuinticFieldData does not support this layout/order");
+    }
+  }
+
+  const FieldMetadata& GetMetadata() const override { return metadata_; }
+
+  Rank1View<const T, HostMemorySpace> GetDOFHolderDataHost() const override
+  {
+    Kokkos::deep_copy(host_data_, device_data_);
+    return make_const_array_view(host_data_);
+  }
+
+  void SetDOFHolderDataHost(Rank1View<const T, HostMemorySpace> values) override
+  {
+    PCMS_ALWAYS_ASSERT(values.size() ==
+                       static_cast<size_t>(layout_->OwnedSize()));
+    CopyHostRank1ViewToDeviceView(device_data_, values);
+    SyncBackend(make_const_array_view(device_data_));
+  }
+
+  Rank1View<const T, DeviceMemorySpace> GetDOFHolderData() const override
+  {
+    return make_const_array_view(device_data_);
+  }
+
+  void SetDOFHolderData(Rank1View<const T, DeviceMemorySpace> values) override
+  {
+    PCMS_ALWAYS_ASSERT(values.size() ==
+                       static_cast<size_t>(layout_->OwnedSize()));
+    CopyDeviceRank1ViewToDeviceView(device_data_, values);
+    SyncBackend(make_const_array_view(device_data_));
+  }
+
+  std::shared_ptr<MeshFieldBackend<T>> GetMeshFieldBackend() const
+  {
+    return mesh_field_;
+  }
+
+private:
+  void SyncBackend(Rank1View<const T, DeviceMemorySpace> flat)
+  {
+    auto nodes_per_dim = layout_->GetNodesPerDim();
+    auto num_components = layout_->GetNumComponents();
+    auto& mesh = layout_->GetMesh();
+    size_t offset = 0;
+    for (int i = 0; i <= mesh.dim(); ++i) {
+      if (nodes_per_dim[i]) {
+        size_t len = static_cast<size_t>(mesh.nents(i)) *
+                     static_cast<size_t>(nodes_per_dim[i]) *
+                     static_cast<size_t>(num_components);
+        Rank1View<const T, DeviceMemorySpace> subspan{
+          flat.data_handle() + offset, len};
+        mesh_field_->SetData(subspan, nodes_per_dim[i], num_components, i);
+        offset += len;
+      }
+    }
+  }
+
+  std::shared_ptr<const MeshFieldsAdapterLayout> layout_;
+  FieldMetadata metadata_;
+  std::shared_ptr<MeshFieldBackend<T>> mesh_field_;
+  mutable Kokkos::View<T*, HostMemorySpace> host_data_;
+  Kokkos::View<T*, DeviceMemorySpace> device_data_;
+};
+
+} // namespace pcms
+
+#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_FIELD_DATA_H

src/pcms/field/evaluator/mesh_fields_backend.h

@@ -35,25 +35,70 @@ class MeshFieldBackend
 };
 
 // ---------------------------------------------------------------------------
-// Concrete backend implementation
+// Evaluation policies — each defines how to create the shape field and
+// evaluate at local coordinates.  This lets a single backend template serve
+// both standard Lagrange and reduced‑quintic elements.
 // ---------------------------------------------------------------------------
-template <typename T, int Dim, int Order>
+
+// Policy for standard Lagrange elements
+template <typename T, int Order>
+struct LagrangeEvalPolicy
+{
+  static constexpr int ShapeOrder = Order;
+  static constexpr int ShapeDOFs  = 1;
+
+  template <typename MeshFieldType, typename ShapeFieldType,
+            typename CoordsView, typename OffsetsView>
+  static auto Evaluate(MeshFieldType& mf, ShapeFieldType& sf,
+                       CoordsView coords, OffsetsView offsets)
+    -> decltype(mf.triangleLocalPointEval(coords, offsets, sf))
+  {
+    return mf.triangleLocalPointEval(coords, offsets, sf);
+  }
+};
+
+// Policy for reduced quintic elements (6 DOFs per vertex, C1 continuity)
+template <typename T>
+struct ReducedQuinticEvalPolicy
+{
+  static constexpr int ShapeOrder = 1;
+  static constexpr int ShapeDOFs  = 6;
+
+  template <typename MeshFieldType, typename ShapeFieldType,
+            typename CoordsView, typename OffsetsView>
+  static auto Evaluate(MeshFieldType& mf, ShapeFieldType& sf,
+                       CoordsView coords, OffsetsView offsets)
+    -> decltype(mf.template triangleReducedQuinticEval<CoordsView,
+                                                       ShapeFieldType>(
+      coords, offsets, sf))
+  {
+    return mf.template triangleReducedQuinticEval<CoordsView, ShapeFieldType>(
+      coords, offsets, sf);
+  }
+};
+
+// ---------------------------------------------------------------------------
+// Single backend implementation (policy‑based)
+// ---------------------------------------------------------------------------
+template <typename T, int Dim, typename EvalPolicy>
 class MeshFieldBackendImpl : public MeshFieldBackend<T>
 {
 public:
   MeshFieldBackendImpl(Omega_h::Mesh& mesh)
     : mesh_(mesh),
       mesh_field_(mesh),
-      shape_field_(mesh_field_.template CreateLagrangeField<T, Order, 1>())
+      shape_field_(
+        mesh_field_.template CreateLagrangeField<T, EvalPolicy::ShapeOrder,
+                                                  EvalPolicy::ShapeDOFs>())
   {
   }
 
   Kokkos::View<T* [1]> evaluate(Kokkos::View<T**> localCoords,
                                 Kokkos::View<LO*> offsets) const override
   {
-    auto self = const_cast<MeshFieldBackendImpl<T, Dim, Order>*>(this);
-    return self->mesh_field_.triangleLocalPointEval(localCoords, offsets,
-                                                    shape_field_);
+    auto self = const_cast<MeshFieldBackendImpl<T, Dim, EvalPolicy>*>(this);
+    return EvalPolicy::Evaluate(self->mesh_field_, self->shape_field_,
+                                localCoords, offsets);
   }
 
   void SetData(Rank1View<const T, DeviceMemorySpace> data, size_t num_nodes,
@@ -92,12 +137,13 @@ class MeshFieldBackendImpl : public MeshFieldBackend<T>
   Omega_h::Mesh& mesh_;
   MeshField::OmegahMeshField<DefaultExecutionSpace, Dim> mesh_field_;
   using ShapeField =
-    decltype(mesh_field_.template CreateLagrangeField<T, Order, 1>());
+    decltype(mesh_field_.template CreateLagrangeField<T, EvalPolicy::ShapeOrder,
+                                                       EvalPolicy::ShapeDOFs>());
   ShapeField shape_field_;
 };
 
 // ---------------------------------------------------------------------------
-// Factory function: create a MeshFieldBackend from a layout
+// Factory function: create a backend from a layout
 // ---------------------------------------------------------------------------
 template <typename T>
 std::shared_ptr<MeshFieldBackend<T>> MakeMeshFieldBackend(
@@ -118,21 +164,64 @@ std::shared_ptr<MeshFieldBackend<T>> MakeMeshFieldBackend(
   if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 0 && nodes_per_dim[2] == 0 &&
       nodes_per_dim[3] == 0) {
     switch (mesh.dim()) {
-      case 1: return std::make_shared<MeshFieldBackendImpl<T, 1, 1>>(mesh);
-      case 2: return std::make_shared<MeshFieldBackendImpl<T, 2, 1>>(mesh);
+      case 1:
+        return std::make_shared<MeshFieldBackendImpl<T, 1,
+                                                     LagrangeEvalPolicy<T, 1>>>(
+          mesh);
+      case 2:
+        return std::make_shared<MeshFieldBackendImpl<T, 2,
+                                                     LagrangeEvalPolicy<T, 1>>>(
+          mesh);
       default: break;
     }
   } else if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 1 &&
              nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) {
     switch (mesh.dim()) {
-      case 2: return std::make_shared<MeshFieldBackendImpl<T, 2, 2>>(mesh);
-      case 3: return std::make_shared<MeshFieldBackendImpl<T, 3, 2>>(mesh);
+      case 2:
+        return std::make_shared<MeshFieldBackendImpl<T, 2,
+                                                     LagrangeEvalPolicy<T, 2>>>(
+          mesh);
+      case 3:
+        return std::make_shared<MeshFieldBackendImpl<T, 3,
+                                                     LagrangeEvalPolicy<T, 2>>>(
+          mesh);
       default: break;
     }
   }
   return nullptr;
 }
 
+// ---------------------------------------------------------------------------
+// Factory function: create a reduced quintic backend from a layout
+// ---------------------------------------------------------------------------
+template <typename T>
+std::shared_ptr<MeshFieldBackend<T>> MakeMeshFieldReducedQuinticBackend(
+  const MeshFieldsAdapterLayout& layout)
+{
+  if constexpr (!std::is_same_v<T, MeshField::Real4> &&
+                !std::is_same_v<T, MeshField::Real8>) {
+    throw pcms_error(
+      "MeshFieldReducedQuinticBackend only supports the MeshFields scalar "
+      "types enabled in this build");
+  }
+
+  Omega_h::Mesh& mesh = layout.GetMesh();
+  if (mesh.dim() != 2) {
+    throw pcms_error(
+      "MeshFieldReducedQuinticBackend only supports 2D meshes");
+  }
+
+  // The reduced quintic backend is only valid for vertex-based layouts
+  // (nodes_per_dim[0] == 1, all others 0).
+  auto nodes_per_dim = layout.GetNodesPerDim();
+  if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 0 &&
+      nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) {
+    return std::make_shared<
+      MeshFieldBackendImpl<T, 2, ReducedQuinticEvalPolicy<T>>>(mesh);
+  }
+  return nullptr;
+}
+
 // ---------------------------------------------------------------------------
 // Helper functors used in MeshFieldsAdapter2LocalizationHint
 // ---------------------------------------------------------------------------
@@ -624,4 +713,4 @@ struct MeshFieldsAdapter2LocalizationHint
 
 } // namespace pcms
 
-#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H
+#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_BACKEND_H