diff --git a/src/pcms/field/CMakeLists.txt b/src/pcms/field/CMakeLists.txt index 49f9b8ff..923d6fdd 100644 --- a/src/pcms/field/CMakeLists.txt +++ b/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 () diff --git a/src/pcms/field/data/mesh_fields_reduced_quintic.h b/src/pcms/field/data/mesh_fields_reduced_quintic.h new file mode 100644 index 00000000..188fa1a0 --- /dev/null +++ b/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 +#include + +namespace pcms +{ + +// MeshFieldsReducedQuinticFieldData is a near-duplicate of +// MeshFieldsFieldData 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 +class MeshFieldsReducedQuinticFieldData : public FieldData +{ +public: + MeshFieldsReducedQuinticFieldData( + std::shared_ptr layout, + FieldMetadata metadata) + : layout_(std::move(layout)), + metadata_(metadata), + mesh_field_(MakeMeshFieldReducedQuinticBackend(*layout_)), + host_data_("meshfields_reduced_quintic_field_data", + static_cast(layout_->OwnedSize())), + device_data_("meshfields_reduced_quintic_field_data_device", + static_cast(layout_->OwnedSize())) + { + if (!mesh_field_) { + throw pcms_error( + "MeshFieldsReducedQuinticFieldData does not support this layout/order"); + } + } + + const FieldMetadata& GetMetadata() const override { return metadata_; } + + Rank1View GetDOFHolderDataHost() const override + { + Kokkos::deep_copy(host_data_, device_data_); + return make_const_array_view(host_data_); + } + + void SetDOFHolderDataHost(Rank1View values) override + { + PCMS_ALWAYS_ASSERT(values.size() == + static_cast(layout_->OwnedSize())); + CopyHostRank1ViewToDeviceView(device_data_, values); + SyncBackend(make_const_array_view(device_data_)); + } + + Rank1View GetDOFHolderData() const override + { + return make_const_array_view(device_data_); + } + + void SetDOFHolderData(Rank1View values) override + { + PCMS_ALWAYS_ASSERT(values.size() == + static_cast(layout_->OwnedSize())); + CopyDeviceRank1ViewToDeviceView(device_data_, values); + SyncBackend(make_const_array_view(device_data_)); + } + + std::shared_ptr> GetMeshFieldBackend() const + { + return mesh_field_; + } + +private: + void SyncBackend(Rank1View 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(mesh.nents(i)) * + static_cast(nodes_per_dim[i]) * + static_cast(num_components); + Rank1View subspan{ + flat.data_handle() + offset, len}; + mesh_field_->SetData(subspan, nodes_per_dim[i], num_components, i); + offset += len; + } + } + } + + std::shared_ptr layout_; + FieldMetadata metadata_; + std::shared_ptr> mesh_field_; + mutable Kokkos::View host_data_; + Kokkos::View device_data_; +}; + +} // namespace pcms + +#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_FIELD_DATA_H \ No newline at end of file diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index fdce6d6c..e382672e 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/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 + +// Policy for standard Lagrange elements +template +struct LagrangeEvalPolicy +{ + static constexpr int ShapeOrder = Order; + static constexpr int ShapeDOFs = 1; + + template + 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 +struct ReducedQuinticEvalPolicy +{ + static constexpr int ShapeOrder = 1; + static constexpr int ShapeDOFs = 6; + + template + static auto Evaluate(MeshFieldType& mf, ShapeFieldType& sf, + CoordsView coords, OffsetsView offsets) + -> decltype(mf.template triangleReducedQuinticEval( + coords, offsets, sf)) + { + return mf.template triangleReducedQuinticEval( + coords, offsets, sf); + } +}; + +// --------------------------------------------------------------------------- +// Single backend implementation (policy‑based) +// --------------------------------------------------------------------------- +template class MeshFieldBackendImpl : public MeshFieldBackend { public: MeshFieldBackendImpl(Omega_h::Mesh& mesh) : mesh_(mesh), mesh_field_(mesh), - shape_field_(mesh_field_.template CreateLagrangeField()) + shape_field_( + mesh_field_.template CreateLagrangeField()) { } Kokkos::View evaluate(Kokkos::View localCoords, Kokkos::View offsets) const override { - auto self = const_cast*>(this); - return self->mesh_field_.triangleLocalPointEval(localCoords, offsets, - shape_field_); + auto self = const_cast*>(this); + return EvalPolicy::Evaluate(self->mesh_field_, self->shape_field_, + localCoords, offsets); } void SetData(Rank1View data, size_t num_nodes, @@ -92,12 +137,13 @@ class MeshFieldBackendImpl : public MeshFieldBackend Omega_h::Mesh& mesh_; MeshField::OmegahMeshField mesh_field_; using ShapeField = - decltype(mesh_field_.template CreateLagrangeField()); + decltype(mesh_field_.template CreateLagrangeField()); ShapeField shape_field_; }; // --------------------------------------------------------------------------- -// Factory function: create a MeshFieldBackend from a layout +// Factory function: create a backend from a layout // --------------------------------------------------------------------------- template std::shared_ptr> MakeMeshFieldBackend( @@ -118,21 +164,64 @@ std::shared_ptr> 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>(mesh); - case 2: return std::make_shared>(mesh); + case 1: + return std::make_shared>>( + mesh); + case 2: + return std::make_shared>>( + 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>(mesh); - case 3: return std::make_shared>(mesh); + case 2: + return std::make_shared>>( + mesh); + case 3: + return std::make_shared>>( + mesh); default: break; } } return nullptr; } +// --------------------------------------------------------------------------- +// Factory function: create a reduced quintic backend from a layout +// --------------------------------------------------------------------------- +template +std::shared_ptr> MakeMeshFieldReducedQuinticBackend( + const MeshFieldsAdapterLayout& layout) +{ + if constexpr (!std::is_same_v && + !std::is_same_v) { + 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>>(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 \ No newline at end of file diff --git a/src/pcms/field/evaluator/mesh_fields_reduced_quintic.h b/src/pcms/field/evaluator/mesh_fields_reduced_quintic.h new file mode 100644 index 00000000..bd86a470 --- /dev/null +++ b/src/pcms/field/evaluator/mesh_fields_reduced_quintic.h @@ -0,0 +1,116 @@ +#ifndef PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_EVALUATOR_FACTORY_H +#define PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_EVALUATOR_FACTORY_H + +#include "pcms/field/evaluator/mesh_fields_backend.h" +#include "pcms/field/evaluator/mesh_fields.h" +#include "pcms/field/data/mesh_fields_reduced_quintic.h" +#include "pcms/field/field_evaluator_factory.h" +#include "pcms/field/out_of_bounds_policy.h" +#include "pcms/field/point_evaluator.h" +#include "pcms/field/field_data.h" +#include "pcms/utility/profile.h" +#include "pcms/utility/arrays.h" +#include + +namespace pcms +{ + +// MeshFieldsReducedQuinticEvaluatorFactory implements +// FieldEvaluatorFactory for MeshFields-backed reduced quintic elements on +// simplex meshes. It owns the spatial search structure and a +// MeshFieldBackend which internally creates the shape field with 6 DOFs +// per vertex for the reduced quintic element. +// +// The reduced quintic element uses 6 DOFs per vertex: +// DOF 0: function value +// DOF 1: df/dx +// DOF 2: df/dy +// DOF 3: d2f/dx2 +// DOF 4: d2f/dxy +// DOF 5: d2f/dy2 +// This provides C1 continuity on triangle meshes. +// +// Note: The PointEvaluator is shared with the standard Lagrange backend +// (MeshFieldsPointEvaluator) since the evaluation strategy is already +// encapsulated in the MeshFieldBackend via the policy-based design. +template +class MeshFieldsReducedQuinticEvaluatorFactory + : public FieldEvaluatorFactory +{ +public: + explicit MeshFieldsReducedQuinticEvaluatorFactory( + std::shared_ptr layout) + : layout_(std::move(layout)), + mesh_(layout_->GetMesh()), + search_(mesh_, 10, 10) + { + if (mesh_.dim() == 3) { + throw pcms_error( + "MeshFieldsReducedQuinticEvaluatorFactory does not support 3D meshes"); + } + if (layout_->GetNumComponents() != 1) { + throw pcms_error( + "MeshFieldsReducedQuinticEvaluatorFactory only supports " + "single-component fields"); + } + } + + const FieldLayout& GetLayout() const override { return *layout_; } + + CoordinateSystem GetCoordinateSystem() const override + { + return layout_->GetDOFHolderCoordinates().GetCoordinateSystem(); + } + + bool HasDOFHolderCoordinates() const override { return true; } + + bool SupportsNearestBoundary() const override { return false; } + + std::unique_ptr> CreatePointEvaluator( + const EvaluationRequest& request) const override + { + PCMS_FUNCTION_TIMER; + const auto coords = request.coords; + const auto policy = request.policy; + if (coords.GetCoordinateSystem() != GetCoordinateSystem()) { + throw pcms_error( + "MeshFieldsEvaluatorFactory: coordinate system mismatch"); + } + if (policy.mode == OutOfBoundsMode::NEAREST_BOUNDARY) { + throw pcms_error( + "MeshFieldsEvaluatorFactory: NearestBoundary is not supported"); + } + + auto coordinates = coords.GetCoordinates(); + Kokkos::View coords_search( + "coords_search", coordinates.extent(0)); + Kokkos::parallel_for( + "copy_coords", + Kokkos::RangePolicy( + 0, coordinates.extent(0)), + KOKKOS_LAMBDA(const int i) { + coords_search(i, 0) = coordinates(i, 0); + coords_search(i, 1) = coordinates(i, 1); + }); + auto results = search_(coords_search); + + MeshFieldsAdapter2LocalizationHint hint(mesh_, results, policy.mode); + + return std::make_unique>( + layout_, std::move(hint), policy.fill_value); + } + + CoordinateView GetDOFHolderCoordinates() const override + { + return layout_->GetDOFHolderCoordinates(); + } + +private: + std::shared_ptr layout_; + Omega_h::Mesh& mesh_; + mutable GridPointSearch2D search_; +}; + +} // namespace pcms + +#endif // PCMS_ADAPTER_MESHFIELDS_MESH_FIELDS_REDUCED_QUINTIC_EVALUATOR_FACTORY_H \ No newline at end of file