1511 lines
54 KiB
C++
1511 lines
54 KiB
C++
module;
|
|
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <concepts>
|
|
#include <cstddef>
|
|
#include <memory>
|
|
#include <optional>
|
|
#include <stdexcept>
|
|
#include <utility>
|
|
|
|
#include <mfem.hpp>
|
|
|
|
export module mean_field:field.mfem;
|
|
|
|
export import :field.registry;
|
|
|
|
namespace mean_field::field::detail {
|
|
template <typename T> inline constexpr bool alwaysFalse = false;
|
|
|
|
// -------------------------------------------------------------------------
|
|
// MFEM finite-element collection construction
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <typename SpaceT> struct FecFor;
|
|
|
|
template <> struct FecFor<L2> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make(
|
|
int familyOrder,
|
|
int dimension
|
|
) {
|
|
return std::make_unique<mfem::L2_FECollection>(familyOrder, dimension);
|
|
}
|
|
};
|
|
|
|
template <> struct FecFor<H1> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make(
|
|
int familyOrder,
|
|
int dimension
|
|
) {
|
|
return std::make_unique<mfem::H1_FECollection>(familyOrder, dimension);
|
|
}
|
|
};
|
|
|
|
template <> struct FecFor<RT> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make(
|
|
int familyOrder,
|
|
int dimension
|
|
) {
|
|
return std::make_unique<mfem::RT_FECollection>(familyOrder, dimension);
|
|
}
|
|
};
|
|
|
|
template <> struct FecFor<ND> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make(
|
|
int familyOrder,
|
|
int dimension
|
|
) {
|
|
return std::make_unique<mfem::ND_FECollection>(familyOrder, dimension);
|
|
}
|
|
};
|
|
|
|
// -------------------------------------------------------------------------
|
|
// MFEM polynomial-order interpretation
|
|
//
|
|
// familyOrder is the collection constructor argument.
|
|
//
|
|
// For RT_p:
|
|
// value order = p + 1
|
|
// divergence order = p
|
|
// normal-trace order = p
|
|
//
|
|
// This distinction is what allows Disc<RT, p> and Disc<L2, p> to form a
|
|
// compatible pair while still giving different value-shape orders.
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <typename OperandT> struct MfemOperandOrder;
|
|
|
|
template <RegisteredQuantity QuantityT, FieldOperationTag OperationT>
|
|
struct MfemOperandOrder<Operand<QuantityT, OperationT>> {
|
|
static constexpr int orderValue = []() consteval {
|
|
if constexpr (GlobalScalarQuantity<QuantityT>) {
|
|
static_assert(
|
|
std::same_as<OperationT, FieldOperation::Value>, "Global scalars support only the value operation."
|
|
);
|
|
|
|
return 0;
|
|
} else {
|
|
using Space = typename QuantityT::Space;
|
|
constexpr int familyOrder = QuantityT::familyOrder;
|
|
|
|
if constexpr (std::same_as<OperationT, FieldOperation::Value>) {
|
|
if constexpr (std::same_as<Space, RT>) {
|
|
return familyOrder + 1;
|
|
} else {
|
|
return familyOrder;
|
|
}
|
|
} else if constexpr (std::same_as<OperationT, FieldOperation::Divergence>) {
|
|
static_assert(
|
|
std::same_as<Space, RT>, "Only RT quantities currently support the divergence "
|
|
"polynomial-order rule."
|
|
);
|
|
|
|
return familyOrder;
|
|
} else if constexpr (std::same_as<OperationT, FieldOperation::Gradient>) {
|
|
static_assert(
|
|
std::same_as<Space, H1>, "Only H1 quantities currently support the gradient "
|
|
"polynomial-order rule."
|
|
);
|
|
|
|
return familyOrder > 0 ? familyOrder - 1 : 0;
|
|
} else if constexpr (std::same_as<OperationT, FieldOperation::Curl>) {
|
|
static_assert(
|
|
std::same_as<Space, ND>, "Only ND quantities currently support the curl "
|
|
"polynomial-order rule."
|
|
);
|
|
|
|
return familyOrder > 0 ? familyOrder - 1 : 0;
|
|
} else if constexpr (std::same_as<OperationT, FieldOperation::NormalTrace>) {
|
|
static_assert(
|
|
std::same_as<Space, RT>, "Only RT quantities currently support the normal-trace "
|
|
"polynomial-order rule."
|
|
);
|
|
|
|
return familyOrder;
|
|
} else {
|
|
static_assert(alwaysFalse<OperationT>, "Unsupported MFEM field operation.");
|
|
}
|
|
}
|
|
}();
|
|
};
|
|
|
|
// -------------------------------------------------------------------------
|
|
// Static polynomial-order contribution of an entire form
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <typename FormT> struct MfemFormOrder;
|
|
|
|
template <auto PolicyKeyV, std::size_t DynamicOrderCountV, FieldOperand... OperandTs>
|
|
struct MfemFormOrder<FormSpec<PolicyKeyV, DynamicOrderCountV, OperandTs...>> {
|
|
static constexpr int staticOrder = (MfemOperandOrder<OperandTs>::orderValue + ... + 0);
|
|
};
|
|
|
|
// -------------------------------------------------------------------------
|
|
// MFEM vector-dimension and ordering rules
|
|
//
|
|
// Vector H1/L2 fields are represented using multiple copies of a scalar
|
|
// finite-element space. RT and ND elements are intrinsically vector-valued
|
|
// and therefore use vdim = 1.
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <FieldQuantity QuantityT> int get_vdim(int spaceDimension) {
|
|
if (spaceDimension <= 0) {
|
|
throw std::invalid_argument("Space dimension must be positive.");
|
|
}
|
|
|
|
if constexpr (QuantityT::rankValue == 0) {
|
|
return 1;
|
|
} else if constexpr (
|
|
std::same_as<typename QuantityT::Space, H1> || std::same_as<typename QuantityT::Space, L2>
|
|
) {
|
|
return spaceDimension;
|
|
} else {
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
template <FieldQuantity QuantityT> constexpr mfem::Ordering::Type get_ordering() {
|
|
if constexpr (
|
|
QuantityT::rankValue == 1 &&
|
|
(std::same_as<typename QuantityT::Space, H1> || std::same_as<typename QuantityT::Space, L2>)
|
|
) {
|
|
return mfem::Ordering::byVDIM;
|
|
} else {
|
|
return mfem::Ordering::byNODES;
|
|
}
|
|
}
|
|
|
|
// -------------------------------------------------------------------------
|
|
// Quantity-specific MFEM realization
|
|
//
|
|
// Backend choices that are part of a field definition live here rather
|
|
// than leaking into FEM setup or call sites.
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <FieldQuantity QuantityT> struct MfemQuantityTraits {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make_fec(int dimension) {
|
|
return FecFor<typename QuantityT::Space>::make(QuantityT::familyOrder, dimension);
|
|
}
|
|
|
|
static constexpr mfem::Ordering::Type ordering = get_ordering<QuantityT>();
|
|
};
|
|
|
|
template <> struct MfemQuantityTraits<Gravity::Flux> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make_fec(int dimension) {
|
|
return std::make_unique<mfem::RT_FECollection>(
|
|
Gravity::Flux::familyOrder, dimension, mfem::BasisType::GaussLobatto, mfem::BasisType::IntegratedGLL
|
|
);
|
|
}
|
|
|
|
static constexpr mfem::Ordering::Type ordering = mfem::Ordering::byNODES;
|
|
};
|
|
|
|
template <> struct MfemQuantityTraits<Displacement::Vector> {
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make_fec(int dimension) {
|
|
return FecFor<H1>::make(Displacement::Vector::familyOrder, dimension);
|
|
}
|
|
|
|
static constexpr mfem::Ordering::Type ordering = mfem::Ordering::byNODES;
|
|
};
|
|
} // namespace mean_field::field::detail
|
|
|
|
export namespace mean_field::field {
|
|
// -------------------------------------------------------------------------
|
|
// User-facing field type
|
|
//
|
|
// The object itself is currently a zero-cost compile-time descriptor:
|
|
//
|
|
// Field<Gravity> gravityField;
|
|
//
|
|
// MFEM construction and typed quadrature-query generation are provided as
|
|
// static operations. Runtime ownership can later be added without changing
|
|
// Gravity, Displacement, or their form definitions.
|
|
// -------------------------------------------------------------------------
|
|
|
|
template <FieldTag TagT> class Field {
|
|
public:
|
|
using Tag = TagT;
|
|
|
|
// ---------------------------------------------------------------------
|
|
// MFEM finite-element collection construction
|
|
// ---------------------------------------------------------------------
|
|
|
|
template <FieldQuantity QuantityT>
|
|
requires typeListContains<
|
|
QuantityT,
|
|
typename TagT::Quantities>
|
|
static std::unique_ptr<mfem::FiniteElementCollection> make_fec(int dimension) {
|
|
if (dimension <= 0) {
|
|
throw std::invalid_argument("Mesh dimension must be positive.");
|
|
}
|
|
|
|
return detail::MfemQuantityTraits<QuantityT>::make_fec(dimension);
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// MFEM parallel finite-element space construction
|
|
//
|
|
// The finite-element collection must outlive the returned space.
|
|
// ---------------------------------------------------------------------
|
|
|
|
template <FieldQuantity QuantityT>
|
|
requires typeListContains<
|
|
QuantityT,
|
|
typename TagT::Quantities>
|
|
static std::unique_ptr<mfem::ParFiniteElementSpace> make_fespace(
|
|
mfem::ParMesh &mesh,
|
|
mfem::FiniteElementCollection &finiteElementCollection
|
|
) {
|
|
return std::make_unique<mfem::ParFiniteElementSpace>(
|
|
&mesh, &finiteElementCollection, detail::get_vdim<QuantityT>(mesh.SpaceDimension()),
|
|
detail::MfemQuantityTraits<QuantityT>::ordering
|
|
);
|
|
}
|
|
|
|
// ---------------------------------------------------------------------
|
|
// Typed quadrature-query construction
|
|
//
|
|
// geometryWeightOrder is supplied at runtime because it depends on the
|
|
// actual element transformation.
|
|
//
|
|
// dynamicOrders contains the form-specific polynomial orders that are
|
|
// not represented by registered compile-time quantities.
|
|
//
|
|
// Examples:
|
|
//
|
|
// Density::Form::CenterOfMass:
|
|
// { positionOrder }
|
|
//
|
|
// Gravity source forms need no dynamic orders because density and
|
|
// potential are both registered quantities.
|
|
//
|
|
// The completed base order is stored in Query::base_order, so Policy
|
|
// does not need to understand divergence, RT conventions, or individual
|
|
// field layouts.
|
|
// ---------------------------------------------------------------------
|
|
|
|
template <FieldForm FormT>
|
|
requires typeListContains<
|
|
FormT,
|
|
typename TagT::FormList>
|
|
static constexpr quadrature::Query make_query(
|
|
quadrature::QuadratureRole role,
|
|
int geometryWeightOrder,
|
|
std::array<
|
|
int,
|
|
FormT::dynamicOrderCount> dynamicOrders = {},
|
|
utils::DOMAINS domain = utils::DOMAINS::ALL,
|
|
quadrature::MappingKind mapping = quadrature::MappingKind::none
|
|
) {
|
|
if (geometryWeightOrder < 0) {
|
|
throw std::invalid_argument("Geometry weight order cannot be negative.");
|
|
}
|
|
|
|
int baseOrder = detail::MfemFormOrder<FormT>::staticOrder + geometryWeightOrder;
|
|
|
|
for (const int dynamicOrder : dynamicOrders) {
|
|
if (dynamicOrder < 0) {
|
|
throw std::invalid_argument("Dynamic polynomial orders cannot be negative.");
|
|
}
|
|
|
|
baseOrder += dynamicOrder;
|
|
}
|
|
|
|
return {
|
|
.term = FormT::policyKey,
|
|
.role = role,
|
|
.domain = domain,
|
|
.mapping = mapping,
|
|
.trial_order = 0,
|
|
.test_order = 0,
|
|
.coefficient_order = 0,
|
|
.geometry_weight_order = geometryWeightOrder,
|
|
.base_order = baseOrder
|
|
};
|
|
}
|
|
};
|
|
|
|
static_assert(FieldTag<Gravity>);
|
|
static_assert(FieldTag<Displacement>);
|
|
static_assert(FieldTag<Density>);
|
|
static_assert(FieldTag<BarotropicConstant>);
|
|
|
|
/*
|
|
* Field-support realization onto MFEM element and DOF indices.
|
|
*
|
|
* A field's compile-time Support is declared in field.registry.
|
|
* These utilities resolve that semantic support through a DomainSchema
|
|
* onto a concrete MFEM finite-element space.
|
|
*
|
|
* Important:
|
|
*
|
|
* active DOFs = union of DOFs touched by supported elements
|
|
*
|
|
* This is deliberately NOT implemented as "remove every DOF touched by
|
|
* an unsupported element". For continuous spaces such as H1, a DOF on
|
|
* the Stellar/Vacuum interface is shared by elements on both sides and
|
|
* remains an active stellar-field DOF.
|
|
*/
|
|
|
|
struct FieldLocalDofSupport {
|
|
/*
|
|
* Marker in local/vector-DOF numbering.
|
|
*
|
|
* Size == finiteElementSpace.GetVSize().
|
|
* Entries are 1 for active DOFs and 0 otherwise.
|
|
*/
|
|
mfem::Array<int> activeVDofMarker;
|
|
|
|
/*
|
|
* Sorted MFEM local/vector DOF indices.
|
|
*/
|
|
mfem::Array<int> activeVDofs;
|
|
mfem::Array<int> inactiveVDofs;
|
|
};
|
|
|
|
struct FieldDofSupport {
|
|
/*
|
|
* Local/vector-DOF information.
|
|
*
|
|
* For a ParFiniteElementSpace the marker is synchronized across
|
|
* neighboring ranks before these lists are constructed, so a shared
|
|
* DOF is active on every rank carrying it if any rank has a supported
|
|
* element touching it.
|
|
*/
|
|
mfem::Array<int> activeVDofMarker;
|
|
mfem::Array<int> activeVDofs;
|
|
mfem::Array<int> inactiveVDofs;
|
|
|
|
/*
|
|
* True-DOF information owned by this MPI rank.
|
|
*
|
|
* Size of activeTrueDofMarker == GetTrueVSize().
|
|
*/
|
|
mfem::Array<int> activeTrueDofMarker;
|
|
mfem::Array<int> activeTrueDofs;
|
|
mfem::Array<int> inactiveTrueDofs;
|
|
};
|
|
|
|
template <typename FieldT>
|
|
concept MfemDomainField = FieldTag<FieldT> && DomainSupportedField<FieldT>;
|
|
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
bool element_is_in_field_support(
|
|
const mfem::Mesh &mesh,
|
|
const int elementId
|
|
) {
|
|
using DomainT = FieldDomainT<FieldT>;
|
|
|
|
static_assert(
|
|
SchemaT::template contains_domain<DomainT>(), "The field support is not completely registered in the "
|
|
"supplied DomainSchema."
|
|
);
|
|
|
|
MFEM_VERIFY(
|
|
elementId >= 0 && elementId < mesh.GetNE(), "The requested field-support element ID is outside the mesh."
|
|
);
|
|
|
|
return SchemaT::template attribute_belongs_to<DomainT>(mesh.GetAttribute(elementId));
|
|
}
|
|
|
|
namespace detail {
|
|
inline void build_marker_lists(
|
|
const mfem::Array<int> &activeMarker,
|
|
mfem::Array<int> &activeDofs,
|
|
mfem::Array<int> &inactiveDofs
|
|
) {
|
|
mfem::FiniteElementSpace::MarkerToList(activeMarker, activeDofs);
|
|
|
|
mfem::Array<int> inactiveMarker(activeMarker.Size());
|
|
|
|
for (int dofId = 0; dofId < activeMarker.Size(); ++dofId) {
|
|
inactiveMarker[dofId] = activeMarker[dofId] == 0 ? 1 : 0;
|
|
}
|
|
|
|
mfem::FiniteElementSpace::MarkerToList(inactiveMarker, inactiveDofs);
|
|
}
|
|
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
mfem::Array<int> build_local_active_vdof_marker(const mfem::FiniteElementSpace &finiteElementSpace) {
|
|
using DomainT = FieldDomainT<FieldT>;
|
|
|
|
static_assert(
|
|
SchemaT::template contains_domain<DomainT>(), "The field support is not completely registered in the "
|
|
"supplied DomainSchema."
|
|
);
|
|
|
|
const mfem::Mesh *mesh = finiteElementSpace.GetMesh();
|
|
|
|
MFEM_VERIFY(mesh != nullptr, "Field-support DOF resolution requires an MFEM mesh.");
|
|
|
|
MFEM_VERIFY(
|
|
finiteElementSpace.GetNE() == mesh->GetNE(), "The finite-element space and mesh have incompatible "
|
|
"element counts."
|
|
);
|
|
|
|
mfem::Array<int> activeMarker(finiteElementSpace.GetVSize());
|
|
|
|
activeMarker = 0;
|
|
|
|
mfem::Array<int> elementVDofs;
|
|
|
|
for (int elementId = 0; elementId < mesh->GetNE(); ++elementId) {
|
|
const int materialId = mesh->GetAttribute(elementId);
|
|
|
|
if (!SchemaT::template attribute_belongs_to<DomainT>(materialId)) {
|
|
continue;
|
|
}
|
|
|
|
finiteElementSpace.GetElementVDofs(elementId, elementVDofs);
|
|
|
|
for (int localIndex = 0; localIndex < elementVDofs.Size(); ++localIndex) {
|
|
/*
|
|
* MFEM can encode orientation in a DOF index by using a
|
|
* negative value. DecodeDof removes that orientation sign
|
|
* and returns the actual local/vector DOF index.
|
|
*/
|
|
const int vdof = mfem::FiniteElementSpace::DecodeDof(elementVDofs[localIndex]);
|
|
|
|
MFEM_VERIFY(
|
|
vdof >= 0 && vdof < finiteElementSpace.GetVSize(),
|
|
"MFEM returned an invalid element vector DOF."
|
|
);
|
|
|
|
activeMarker[vdof] = 1;
|
|
}
|
|
}
|
|
|
|
return activeMarker;
|
|
}
|
|
} // namespace detail
|
|
|
|
/*
|
|
* Serial/local support resolution.
|
|
*
|
|
* This works with any mfem::FiniteElementSpace and is particularly
|
|
* useful for topology/unit tests.
|
|
*
|
|
* The returned indices use MFEM local/vector-DOF numbering, not
|
|
* true-DOF numbering.
|
|
*/
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
FieldLocalDofSupport resolve_field_local_dof_support(const mfem::FiniteElementSpace &finiteElementSpace) {
|
|
FieldLocalDofSupport result;
|
|
|
|
result.activeVDofMarker = detail::build_local_active_vdof_marker<FieldT, SchemaT>(finiteElementSpace);
|
|
|
|
detail::build_marker_lists(result.activeVDofMarker, result.activeVDofs, result.inactiveVDofs);
|
|
|
|
return result;
|
|
}
|
|
|
|
/*
|
|
* Parallel production support resolution.
|
|
*
|
|
* This additionally converts the field support to the locally-owned
|
|
* true-DOF numbering used by nonlinear vectors and operators.
|
|
*
|
|
* For now this intentionally requires a conforming ParFiniteElementSpace.
|
|
* MFEM's nonconforming spaces require an additional constraint/conforming-
|
|
* DOF projection step; silently treating their local DOFs as ordinary
|
|
* true DOFs would be incorrect.
|
|
*/
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
FieldDofSupport resolve_field_dof_support(const mfem::ParFiniteElementSpace &finiteElementSpace) {
|
|
FieldDofSupport result;
|
|
|
|
MFEM_VERIFY(
|
|
!finiteElementSpace.Nonconforming(), "Field-support true-DOF resolution currently requires a "
|
|
"conforming mfem::ParFiniteElementSpace."
|
|
);
|
|
|
|
result.activeVDofMarker = detail::build_local_active_vdof_marker<FieldT, SchemaT>(finiteElementSpace);
|
|
|
|
/*
|
|
* Shared H1/RT DOFs can lie on an MPI partition boundary.
|
|
*
|
|
* If a supported element exists on one rank and the shared DOF also
|
|
* exists on a neighboring rank whose local elements are unsupported,
|
|
* that DOF must nevertheless be active globally.
|
|
*
|
|
* MFEM Synchronize performs the required OR-like synchronization of
|
|
* the marker across shared local DOFs.
|
|
*/
|
|
finiteElementSpace.Synchronize(result.activeVDofMarker);
|
|
|
|
detail::build_marker_lists(result.activeVDofMarker, result.activeVDofs, result.inactiveVDofs);
|
|
|
|
result.activeTrueDofMarker.SetSize(finiteElementSpace.GetTrueVSize());
|
|
|
|
result.activeTrueDofMarker = 0;
|
|
|
|
for (int vdof = 0; vdof < result.activeVDofMarker.Size(); ++vdof) {
|
|
if (result.activeVDofMarker[vdof] == 0) {
|
|
continue;
|
|
}
|
|
|
|
/*
|
|
* GetLocalTDofNumber returns the locally-owned true-DOF index
|
|
* for this local/vector DOF, or -1 when this rank does not own
|
|
* the shared true DOF.
|
|
*
|
|
* Because activeVDofMarker was synchronized first, the owning
|
|
* rank will also see the active marker.
|
|
*/
|
|
const int trueDof = finiteElementSpace.GetLocalTDofNumber(vdof);
|
|
|
|
if (trueDof < 0) {
|
|
continue;
|
|
}
|
|
|
|
MFEM_VERIFY(trueDof < result.activeTrueDofMarker.Size(), "MFEM returned an invalid local true DOF.");
|
|
|
|
result.activeTrueDofMarker[trueDof] = 1;
|
|
}
|
|
|
|
detail::build_marker_lists(result.activeTrueDofMarker, result.activeTrueDofs, result.inactiveTrueDofs);
|
|
|
|
return result;
|
|
}
|
|
|
|
/*
|
|
* Canonical correspondence between a dense reduced field vector and
|
|
* the selected MFEM true DOFs representing that field.
|
|
*
|
|
* The map contains no field, domain, mesh, or solver policy. It is an
|
|
* immutable indexing object once constructed:
|
|
*
|
|
* reduced index i
|
|
* |
|
|
* v
|
|
* reducedToTrue[i]
|
|
* |
|
|
* v
|
|
* MFEM true DOF
|
|
*
|
|
* trueToReduced supplies the inverse map. Unsupported true DOFs carry
|
|
* the sentinel -1.
|
|
*
|
|
* The reduced-to-true list is required to be strictly increasing.
|
|
* This makes reduced ordering deterministic and agrees with the
|
|
* canonical ordering produced by MFEM MarkerToList().
|
|
*/
|
|
class FieldDofMap {
|
|
public:
|
|
FieldDofMap() = default;
|
|
|
|
FieldDofMap(
|
|
const int fullTrueDofSize,
|
|
const mfem::Array<int> &reducedToTrue
|
|
) {
|
|
if (fullTrueDofSize < 0) {
|
|
throw std::invalid_argument("FieldDofMap requires a non-negative full true-DOF size.");
|
|
}
|
|
|
|
m_fullTrueDofSize = fullTrueDofSize;
|
|
|
|
m_reducedToTrue.SetSize(reducedToTrue.Size());
|
|
|
|
m_trueToReduced.SetSize(m_fullTrueDofSize);
|
|
|
|
m_trueToReduced = -1;
|
|
|
|
int previousTrueDof = -1;
|
|
|
|
for (int reducedDof = 0; reducedDof < reducedToTrue.Size(); ++reducedDof) {
|
|
const int trueDof = reducedToTrue[reducedDof];
|
|
|
|
if (trueDof < 0 || trueDof >= m_fullTrueDofSize) {
|
|
throw std::invalid_argument(
|
|
"FieldDofMap contains a true DOF outside the full "
|
|
"true-DOF space."
|
|
);
|
|
}
|
|
|
|
if (reducedDof > 0 && trueDof <= previousTrueDof) {
|
|
throw std::invalid_argument(
|
|
"FieldDofMap reduced-to-true indices must be "
|
|
"strictly increasing and unique."
|
|
);
|
|
}
|
|
|
|
m_reducedToTrue[reducedDof] = trueDof;
|
|
|
|
m_trueToReduced[trueDof] = reducedDof;
|
|
|
|
previousTrueDof = trueDof;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Construct directly from the support result produced by
|
|
* resolve_field_dof_support().
|
|
*
|
|
* The marker is checked against the active true-DOF list so that
|
|
* an internally inconsistent FieldDofSupport cannot silently
|
|
* produce a solver map.
|
|
*/
|
|
explicit FieldDofMap(const FieldDofSupport &support)
|
|
: FieldDofMap(
|
|
support.activeTrueDofMarker.Size(),
|
|
support.activeTrueDofs
|
|
) {
|
|
for (int trueDof = 0; trueDof < m_fullTrueDofSize; ++trueDof) {
|
|
const bool markerSaysActive = support.activeTrueDofMarker[trueDof] != 0;
|
|
|
|
const bool mapSaysActive = m_trueToReduced[trueDof] >= 0;
|
|
|
|
if (markerSaysActive != mapSaysActive) {
|
|
throw std::invalid_argument(
|
|
"FieldDofSupport active marker and active true-DOF "
|
|
"list are inconsistent."
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
[[nodiscard]]
|
|
int full_size() const noexcept {
|
|
return m_fullTrueDofSize;
|
|
}
|
|
|
|
[[nodiscard]]
|
|
int reduced_size() const noexcept {
|
|
return m_reducedToTrue.Size();
|
|
}
|
|
|
|
[[nodiscard]]
|
|
int inactive_size() const noexcept {
|
|
return full_size() - reduced_size();
|
|
}
|
|
|
|
/*
|
|
* Because reducedToTrue is strictly increasing, a map containing
|
|
* every true DOF necessarily has
|
|
*
|
|
* reducedToTrue[i] == i.
|
|
*/
|
|
[[nodiscard]]
|
|
bool is_identity() const noexcept {
|
|
return reduced_size() == full_size();
|
|
}
|
|
|
|
[[nodiscard]]
|
|
const mfem::Array<int> &reduced_to_true() const noexcept {
|
|
return m_reducedToTrue;
|
|
}
|
|
|
|
/*
|
|
* Values are:
|
|
*
|
|
* >= 0 reduced DOF index
|
|
* -1 unsupported/inactive true DOF
|
|
*/
|
|
[[nodiscard]]
|
|
const mfem::Array<int> &true_to_reduced() const noexcept {
|
|
return m_trueToReduced;
|
|
}
|
|
|
|
[[nodiscard]]
|
|
bool contains_true_dof(const int trueDof) const {
|
|
validate_true_dof(trueDof);
|
|
|
|
return m_trueToReduced[trueDof] >= 0;
|
|
}
|
|
|
|
[[nodiscard]]
|
|
int true_dof(const int reducedDof) const {
|
|
if (reducedDof < 0 || reducedDof >= reduced_size()) {
|
|
throw std::out_of_range("Reduced DOF index is outside FieldDofMap.");
|
|
}
|
|
|
|
return m_reducedToTrue[reducedDof];
|
|
}
|
|
|
|
[[nodiscard]]
|
|
std::optional<int> reduced_dof(const int trueDof) const {
|
|
validate_true_dof(trueDof);
|
|
|
|
const int reducedDof = m_trueToReduced[trueDof];
|
|
|
|
if (reducedDof < 0) {
|
|
return std::nullopt;
|
|
}
|
|
|
|
return reducedDof;
|
|
}
|
|
|
|
/*
|
|
* Gather:
|
|
*
|
|
* full MFEM true vector
|
|
* |
|
|
* v
|
|
* dense reduced solver vector
|
|
*/
|
|
void gather(
|
|
const mfem::Vector &full,
|
|
mfem::Vector &reduced
|
|
) const {
|
|
require_full_size(full);
|
|
|
|
require_reduced_size(reduced);
|
|
|
|
if (is_identity()) {
|
|
reduced = full;
|
|
return;
|
|
}
|
|
|
|
for (int reducedDof = 0; reducedDof < reduced_size(); ++reducedDof) {
|
|
reduced(reducedDof) = full(m_reducedToTrue[reducedDof]);
|
|
}
|
|
}
|
|
|
|
[[nodiscard]]
|
|
mfem::Vector gather(const mfem::Vector &full) const {
|
|
mfem::Vector reduced(reduced_size());
|
|
|
|
gather(full, reduced);
|
|
|
|
return reduced;
|
|
}
|
|
|
|
/*
|
|
* Scatter with projection semantics.
|
|
*
|
|
* All unsupported true DOFs are explicitly zeroed.
|
|
*
|
|
* This is the normal operation for constructing a complete MFEM
|
|
* representation of a supported field from the reduced nonlinear
|
|
* state.
|
|
*
|
|
* The output vector is NOT resized. This is intentional: callers
|
|
* may provide an mfem::Vector view into an mfem::BlockVector.
|
|
*/
|
|
void scatter(
|
|
const mfem::Vector &reduced,
|
|
mfem::Vector &full
|
|
) const {
|
|
require_reduced_size(reduced);
|
|
|
|
require_full_size(full);
|
|
|
|
if (is_identity()) {
|
|
full = reduced;
|
|
return;
|
|
}
|
|
|
|
full = 0.0;
|
|
|
|
scatter_into(reduced, full);
|
|
}
|
|
|
|
[[nodiscard]]
|
|
mfem::Vector scatter(const mfem::Vector &reduced) const {
|
|
mfem::Vector full(full_size());
|
|
|
|
scatter(reduced, full);
|
|
|
|
return full;
|
|
}
|
|
|
|
/*
|
|
* Scatter while preserving unsupported values already present in
|
|
* the full vector.
|
|
*
|
|
* This is distinct from scatter() because future constrained field
|
|
* representations may need to preserve prescribed values outside
|
|
* the current reduced/free set.
|
|
*/
|
|
void scatter_into(
|
|
const mfem::Vector &reduced,
|
|
mfem::Vector &full
|
|
) const {
|
|
require_reduced_size(reduced);
|
|
|
|
require_full_size(full);
|
|
|
|
if (is_identity()) {
|
|
full = reduced;
|
|
return;
|
|
}
|
|
|
|
for (int reducedDof = 0; reducedDof < reduced_size(); ++reducedDof) {
|
|
full(m_reducedToTrue[reducedDof]) = reduced(reducedDof);
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Add a reduced vector into the selected true DOFs.
|
|
*
|
|
* Unsupported true DOFs are untouched.
|
|
*/
|
|
void scatter_add(
|
|
const mfem::Vector &reduced,
|
|
mfem::Vector &full,
|
|
const double scale = 1.0
|
|
) const {
|
|
require_reduced_size(reduced);
|
|
|
|
require_full_size(full);
|
|
|
|
if (is_identity()) {
|
|
full.Add(scale, reduced);
|
|
return;
|
|
}
|
|
|
|
for (int reducedDof = 0; reducedDof < reduced_size(); ++reducedDof) {
|
|
full(m_reducedToTrue[reducedDof]) += scale * reduced(reducedDof);
|
|
}
|
|
}
|
|
|
|
private:
|
|
void validate_true_dof(const int trueDof) const {
|
|
if (trueDof < 0 || trueDof >= full_size()) {
|
|
throw std::out_of_range("True DOF index is outside FieldDofMap.");
|
|
}
|
|
}
|
|
|
|
void require_full_size(const mfem::Vector &vector) const {
|
|
if (vector.Size() != full_size()) {
|
|
throw std::invalid_argument("FieldDofMap full vector has an incompatible size.");
|
|
}
|
|
}
|
|
|
|
void require_reduced_size(const mfem::Vector &vector) const {
|
|
if (vector.Size() != reduced_size()) {
|
|
throw std::invalid_argument("FieldDofMap reduced vector has an incompatible size.");
|
|
}
|
|
}
|
|
|
|
int m_fullTrueDofSize{0};
|
|
|
|
/*
|
|
* Canonical forward mapping:
|
|
*
|
|
* reduced -> MFEM true
|
|
*/
|
|
mfem::Array<int> m_reducedToTrue;
|
|
|
|
/*
|
|
* Inverse mapping:
|
|
*
|
|
* MFEM true -> reduced
|
|
*
|
|
* Unsupported true DOFs are -1.
|
|
*/
|
|
mfem::Array<int> m_trueToReduced;
|
|
};
|
|
|
|
/*
|
|
* Canonical coordinates on a scalar finite-element boundary.
|
|
*
|
|
* Unlike FieldBoundaryDofMap, this is not a row selection inside an
|
|
* existing physical field. It defines an independent dense coordinate
|
|
* vector whose entries are the locally owned scalar true DOFs on a
|
|
* semantic boundary.
|
|
*
|
|
* Local coordinates follow increasing MFEM true-DOF order. Global
|
|
* coordinates use the distributed-vector convention: ranks are ordered by
|
|
* communicator rank and each rank contributes its locally sorted block.
|
|
*/
|
|
class ScalarBoundaryDofMap final {
|
|
public:
|
|
ScalarBoundaryDofMap() = default;
|
|
|
|
ScalarBoundaryDofMap(
|
|
const int volumeTrueDofSize,
|
|
const mfem::Array<int> &boundaryTrueDofs,
|
|
const long long globalOffset,
|
|
const long long globalSize
|
|
)
|
|
: m_boundaryDofs(
|
|
volumeTrueDofSize,
|
|
boundaryTrueDofs
|
|
),
|
|
m_globalOffset(globalOffset),
|
|
m_globalSize(globalSize) {
|
|
if (m_globalOffset < 0) {
|
|
throw std::invalid_argument("ScalarBoundaryDofMap requires a non-negative global offset.");
|
|
}
|
|
if (m_globalSize < 0) {
|
|
throw std::invalid_argument("ScalarBoundaryDofMap requires a non-negative global size.");
|
|
}
|
|
if (m_globalOffset + local_size() > m_globalSize) {
|
|
throw std::invalid_argument(
|
|
"ScalarBoundaryDofMap local coordinates lie outside the global boundary coordinate vector."
|
|
);
|
|
}
|
|
}
|
|
|
|
[[nodiscard]] int volume_true_dof_size() const noexcept {
|
|
return m_boundaryDofs.full_size();
|
|
}
|
|
|
|
[[nodiscard]] int local_size() const noexcept {
|
|
return m_boundaryDofs.reduced_size();
|
|
}
|
|
|
|
[[nodiscard]] long long global_size() const noexcept {
|
|
return m_globalSize;
|
|
}
|
|
|
|
[[nodiscard]] long long global_offset() const noexcept {
|
|
return m_globalOffset;
|
|
}
|
|
|
|
[[nodiscard]] bool empty() const noexcept {
|
|
return local_size() == 0;
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Array<int> &boundary_true_dofs() const noexcept {
|
|
return m_boundaryDofs.reduced_to_true();
|
|
}
|
|
|
|
[[nodiscard]] int volume_true_dof(const int localBoundaryDof) const {
|
|
return m_boundaryDofs.true_dof(localBoundaryDof);
|
|
}
|
|
|
|
[[nodiscard]] std::optional<int> local_boundary_dof(const int volumeTrueDof) const {
|
|
return m_boundaryDofs.reduced_dof(volumeTrueDof);
|
|
}
|
|
|
|
[[nodiscard]] bool contains_volume_true_dof(const int volumeTrueDof) const {
|
|
return m_boundaryDofs.contains_true_dof(volumeTrueDof);
|
|
}
|
|
|
|
[[nodiscard]] long long global_boundary_dof(const int localBoundaryDof) const {
|
|
if (localBoundaryDof < 0 || localBoundaryDof >= local_size()) {
|
|
throw std::out_of_range("Local boundary DOF is outside ScalarBoundaryDofMap.");
|
|
}
|
|
return m_globalOffset + localBoundaryDof;
|
|
}
|
|
|
|
void gather(
|
|
const mfem::Vector &volumeTrueValues,
|
|
mfem::Vector &boundaryValues
|
|
) const {
|
|
m_boundaryDofs.gather(volumeTrueValues, boundaryValues);
|
|
}
|
|
|
|
[[nodiscard]] mfem::Vector gather(const mfem::Vector &volumeTrueValues) const {
|
|
return m_boundaryDofs.gather(volumeTrueValues);
|
|
}
|
|
|
|
void scatter(
|
|
const mfem::Vector &boundaryValues,
|
|
mfem::Vector &volumeTrueValues
|
|
) const {
|
|
m_boundaryDofs.scatter(boundaryValues, volumeTrueValues);
|
|
}
|
|
|
|
[[nodiscard]] mfem::Vector scatter(const mfem::Vector &boundaryValues) const {
|
|
return m_boundaryDofs.scatter(boundaryValues);
|
|
}
|
|
|
|
private:
|
|
FieldDofMap m_boundaryDofs;
|
|
long long m_globalOffset{0};
|
|
long long m_globalSize{0};
|
|
};
|
|
|
|
template <
|
|
utils::domain::IsBoundary BoundaryT,
|
|
utils::domain::IsSchema SchemaT>
|
|
requires(SchemaT::template contains_boundary<BoundaryT>())
|
|
[[nodiscard]] ScalarBoundaryDofMap
|
|
make_scalar_boundary_dof_map(const mfem::ParFiniteElementSpace &finiteElementSpace) {
|
|
MFEM_VERIFY(
|
|
!finiteElementSpace.Nonconforming(),
|
|
"Scalar boundary true-DOF resolution currently requires a conforming mfem::ParFiniteElementSpace."
|
|
);
|
|
MFEM_VERIFY(finiteElementSpace.GetVDim() == 1, "ScalarBoundaryDofMap requires a scalar finite-element space.");
|
|
|
|
const mfem::Mesh *mesh = finiteElementSpace.GetMesh();
|
|
MFEM_VERIFY(mesh != nullptr, "Scalar boundary DOF resolution requires an MFEM mesh.");
|
|
|
|
mfem::Array<int> boundaryVDofMarker(finiteElementSpace.GetVSize());
|
|
boundaryVDofMarker = 0;
|
|
|
|
mfem::Array<int> boundaryElementVDofs;
|
|
for (int boundaryElement = 0; boundaryElement < mesh->GetNBE(); ++boundaryElement) {
|
|
if (!SchemaT::template boundary_attribute_matches<BoundaryT>(mesh->GetBdrAttribute(boundaryElement))) {
|
|
continue;
|
|
}
|
|
|
|
finiteElementSpace.GetBdrElementVDofs(boundaryElement, boundaryElementVDofs);
|
|
for (const int encodedVDof : boundaryElementVDofs) {
|
|
const int vdof = mfem::FiniteElementSpace::DecodeDof(encodedVDof);
|
|
MFEM_VERIFY(
|
|
vdof >= 0 && vdof < finiteElementSpace.GetVSize(),
|
|
"MFEM returned an invalid scalar boundary vector DOF."
|
|
);
|
|
boundaryVDofMarker[vdof] = 1;
|
|
}
|
|
}
|
|
|
|
finiteElementSpace.Synchronize(boundaryVDofMarker);
|
|
|
|
mfem::Array<int> boundaryTrueDofMarker(finiteElementSpace.GetTrueVSize());
|
|
boundaryTrueDofMarker = 0;
|
|
|
|
for (int vdof = 0; vdof < boundaryVDofMarker.Size(); ++vdof) {
|
|
if (boundaryVDofMarker[vdof] == 0) {
|
|
continue;
|
|
}
|
|
|
|
const int trueDof = finiteElementSpace.GetLocalTDofNumber(vdof);
|
|
if (trueDof < 0) {
|
|
continue;
|
|
}
|
|
|
|
MFEM_VERIFY(trueDof < boundaryTrueDofMarker.Size(), "MFEM returned an invalid scalar boundary true DOF.");
|
|
boundaryTrueDofMarker[trueDof] = 1;
|
|
}
|
|
|
|
mfem::Array<int> boundaryTrueDofs;
|
|
mfem::FiniteElementSpace::MarkerToList(boundaryTrueDofMarker, boundaryTrueDofs);
|
|
|
|
const long long localSize = boundaryTrueDofs.Size();
|
|
long long globalSize = 0;
|
|
long long globalOffset = 0;
|
|
|
|
MPI_Allreduce(&localSize, &globalSize, 1, MPI_LONG_LONG, MPI_SUM, finiteElementSpace.GetComm());
|
|
MPI_Exscan(&localSize, &globalOffset, 1, MPI_LONG_LONG, MPI_SUM, finiteElementSpace.GetComm());
|
|
|
|
if (finiteElementSpace.GetMyRank() == 0) {
|
|
globalOffset = 0;
|
|
}
|
|
|
|
MFEM_VERIFY(globalSize > 0, "The requested semantic boundary has no scalar true DOFs.");
|
|
|
|
return ScalarBoundaryDofMap(finiteElementSpace.GetTrueVSize(), boundaryTrueDofs, globalOffset, globalSize);
|
|
}
|
|
|
|
template <utils::domain::IsSchema SchemaT = utils::domain::CoreEnvelopeVacuumDomainSchema>
|
|
requires(SchemaT::template contains_boundary<utils::domain::StellarSurface>())
|
|
[[nodiscard]] ScalarBoundaryDofMap
|
|
make_stellar_surface_scalar_dof_map(const mfem::ParFiniteElementSpace &finiteElementSpace) {
|
|
return make_scalar_boundary_dof_map<utils::domain::StellarSurface, SchemaT>(finiteElementSpace);
|
|
}
|
|
|
|
/*
|
|
* Boundary rows expressed in a field's reduced solver ordering.
|
|
*
|
|
* This object is deliberately independent of any particular physical
|
|
* surface condition. Its template constructor below combines a field,
|
|
* a semantic boundary, and a domain schema. Consequently the same
|
|
* topology machinery can be used by any compiled surface formulation;
|
|
* it is not tied to enthalpy or pressure.
|
|
*/
|
|
class FieldBoundaryDofMap final {
|
|
public:
|
|
FieldBoundaryDofMap() = default;
|
|
|
|
FieldBoundaryDofMap(
|
|
const int fieldReducedSize,
|
|
const mfem::Array<int> &boundaryReducedDofs
|
|
)
|
|
: m_fieldReducedSize(fieldReducedSize),
|
|
m_boundaryReducedDofs(boundaryReducedDofs) {
|
|
if (m_fieldReducedSize < 0) {
|
|
throw std::invalid_argument("FieldBoundaryDofMap requires a non-negative field size.");
|
|
}
|
|
|
|
m_boundaryReducedDofMarker.SetSize(m_fieldReducedSize);
|
|
m_boundaryReducedDofMarker = 0;
|
|
|
|
int previousReducedDof = -1;
|
|
for (const int reducedDof : m_boundaryReducedDofs) {
|
|
if (reducedDof < 0 || reducedDof >= m_fieldReducedSize) {
|
|
throw std::invalid_argument("FieldBoundaryDofMap contains a DOF outside the reduced field vector.");
|
|
}
|
|
if (reducedDof <= previousReducedDof) {
|
|
throw std::invalid_argument("FieldBoundaryDofMap indices must be strictly increasing and unique.");
|
|
}
|
|
|
|
m_boundaryReducedDofMarker[reducedDof] = 1;
|
|
previousReducedDof = reducedDof;
|
|
}
|
|
}
|
|
|
|
[[nodiscard]] int field_size() const noexcept {
|
|
return m_fieldReducedSize;
|
|
}
|
|
|
|
[[nodiscard]] int size() const noexcept {
|
|
return m_boundaryReducedDofs.Size();
|
|
}
|
|
|
|
[[nodiscard]] bool empty() const noexcept {
|
|
return size() == 0;
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Array<int> &reduced_dofs() const noexcept {
|
|
return m_boundaryReducedDofs;
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Array<int> &reduced_dof_marker() const noexcept {
|
|
return m_boundaryReducedDofMarker;
|
|
}
|
|
|
|
[[nodiscard]] bool contains(const int reducedDof) const {
|
|
if (reducedDof < 0 || reducedDof >= m_fieldReducedSize) {
|
|
throw std::out_of_range("Reduced DOF index is outside FieldBoundaryDofMap.");
|
|
}
|
|
return m_boundaryReducedDofMarker[reducedDof] != 0;
|
|
}
|
|
|
|
private:
|
|
int m_fieldReducedSize{0};
|
|
mfem::Array<int> m_boundaryReducedDofs;
|
|
mfem::Array<int> m_boundaryReducedDofMarker;
|
|
};
|
|
|
|
/* Point-supported rows in a field's reduced solver ordering. */
|
|
class FieldPointDofMap final {
|
|
public:
|
|
FieldPointDofMap() = default;
|
|
|
|
FieldPointDofMap(
|
|
const int fieldReducedSize,
|
|
const mfem::Array<int> &pointReducedDofs
|
|
)
|
|
: m_selectedDofs(
|
|
fieldReducedSize,
|
|
pointReducedDofs
|
|
) {
|
|
}
|
|
|
|
[[nodiscard]] int field_size() const noexcept {
|
|
return m_selectedDofs.field_size();
|
|
}
|
|
|
|
[[nodiscard]] int size() const noexcept {
|
|
return m_selectedDofs.size();
|
|
}
|
|
|
|
[[nodiscard]] bool empty() const noexcept {
|
|
return m_selectedDofs.empty();
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Array<int> &reduced_dofs() const noexcept {
|
|
return m_selectedDofs.reduced_dofs();
|
|
}
|
|
|
|
[[nodiscard]] const mfem::Array<int> &reduced_dof_marker() const noexcept {
|
|
return m_selectedDofs.reduced_dof_marker();
|
|
}
|
|
|
|
[[nodiscard]] bool contains(const int reducedDof) const {
|
|
return m_selectedDofs.contains(reducedDof);
|
|
}
|
|
|
|
private:
|
|
FieldBoundaryDofMap m_selectedDofs;
|
|
};
|
|
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsBoundary BoundaryT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]] FieldBoundaryDofMap make_field_boundary_dof_map(
|
|
const mfem::ParFiniteElementSpace &finiteElementSpace,
|
|
const FieldDofMap &fieldDofMap
|
|
) {
|
|
static_assert(
|
|
SchemaT::template contains_boundary<BoundaryT>(),
|
|
"The requested boundary is not registered in the supplied DomainSchema."
|
|
);
|
|
|
|
MFEM_VERIFY(
|
|
!finiteElementSpace.Nonconforming(),
|
|
"Field boundary true-DOF resolution currently requires a conforming mfem::ParFiniteElementSpace."
|
|
);
|
|
MFEM_VERIFY(
|
|
fieldDofMap.full_size() == finiteElementSpace.GetTrueVSize(),
|
|
"The field map and finite-element space have incompatible true-DOF sizes."
|
|
);
|
|
|
|
const mfem::Mesh *mesh = finiteElementSpace.GetMesh();
|
|
MFEM_VERIFY(mesh != nullptr, "Field boundary DOF resolution requires an MFEM mesh.");
|
|
|
|
mfem::Array<int> boundaryVDofMarker(finiteElementSpace.GetVSize());
|
|
boundaryVDofMarker = 0;
|
|
|
|
mfem::Array<int> boundaryElementVDofs;
|
|
for (int boundaryElement = 0; boundaryElement < mesh->GetNBE(); ++boundaryElement) {
|
|
if (!SchemaT::template boundary_attribute_matches<BoundaryT>(mesh->GetBdrAttribute(boundaryElement))) {
|
|
continue;
|
|
}
|
|
|
|
finiteElementSpace.GetBdrElementVDofs(boundaryElement, boundaryElementVDofs);
|
|
for (const int encodedVDof : boundaryElementVDofs) {
|
|
const int vdof = mfem::FiniteElementSpace::DecodeDof(encodedVDof);
|
|
MFEM_VERIFY(
|
|
vdof >= 0 && vdof < finiteElementSpace.GetVSize(), "MFEM returned an invalid boundary vector DOF."
|
|
);
|
|
boundaryVDofMarker[vdof] = 1;
|
|
}
|
|
}
|
|
|
|
finiteElementSpace.Synchronize(boundaryVDofMarker);
|
|
|
|
mfem::Array<int> boundaryReducedDofMarker(fieldDofMap.reduced_size());
|
|
boundaryReducedDofMarker = 0;
|
|
|
|
for (int vdof = 0; vdof < boundaryVDofMarker.Size(); ++vdof) {
|
|
if (boundaryVDofMarker[vdof] == 0) {
|
|
continue;
|
|
}
|
|
|
|
const int trueDof = finiteElementSpace.GetLocalTDofNumber(vdof);
|
|
if (trueDof < 0) {
|
|
continue;
|
|
}
|
|
|
|
const std::optional<int> reducedDof = fieldDofMap.reduced_dof(trueDof);
|
|
MFEM_VERIFY(
|
|
reducedDof.has_value(),
|
|
"A boundary DOF selected for the field is absent from that field's reduced solver map."
|
|
);
|
|
boundaryReducedDofMarker[*reducedDof] = 1;
|
|
}
|
|
|
|
mfem::Array<int> boundaryReducedDofs;
|
|
mfem::FiniteElementSpace::MarkerToList(boundaryReducedDofMarker, boundaryReducedDofs);
|
|
return FieldBoundaryDofMap(fieldDofMap.reduced_size(), boundaryReducedDofs);
|
|
}
|
|
|
|
template <MfemDomainField FieldT>
|
|
[[nodiscard]] FieldPointDofMap make_field_point_dof_map(
|
|
const mfem::ParFiniteElementSpace &finiteElementSpace,
|
|
const FieldDofMap &fieldDofMap,
|
|
const mfem::Vector &point,
|
|
const double tolerance
|
|
) {
|
|
MFEM_VERIFY(
|
|
!finiteElementSpace.Nonconforming(),
|
|
"Field point true-DOF resolution currently requires a conforming mfem::ParFiniteElementSpace."
|
|
);
|
|
MFEM_VERIFY(
|
|
fieldDofMap.full_size() == finiteElementSpace.GetTrueVSize(),
|
|
"The field map and finite-element space have incompatible true-DOF sizes."
|
|
);
|
|
MFEM_VERIFY(
|
|
std::isfinite(tolerance) && tolerance >= 0.0, "The field point tolerance must be finite and non-negative."
|
|
);
|
|
|
|
const mfem::Mesh *mesh = finiteElementSpace.GetMesh();
|
|
MFEM_VERIFY(mesh != nullptr, "Field point DOF resolution requires an MFEM mesh.");
|
|
MFEM_VERIFY(
|
|
point.Size() == mesh->SpaceDimension(), "The requested field point has the wrong coordinate dimension."
|
|
);
|
|
|
|
mfem::Array<int> pointVDofMarker(finiteElementSpace.GetVSize());
|
|
pointVDofMarker = 0;
|
|
|
|
mfem::Array<int> vertexVDofs;
|
|
for (int vertex = 0; vertex < mesh->GetNV(); ++vertex) {
|
|
const mfem::real_t *coordinates = mesh->GetVertex(vertex);
|
|
double distanceSquared = 0.0;
|
|
for (int component = 0; component < point.Size(); ++component) {
|
|
const double difference = coordinates[component] - point(component);
|
|
distanceSquared += difference * difference;
|
|
}
|
|
if (std::sqrt(distanceSquared) > tolerance) {
|
|
continue;
|
|
}
|
|
|
|
finiteElementSpace.GetVertexVDofs(vertex, vertexVDofs);
|
|
for (const int encodedVDof : vertexVDofs) {
|
|
const int vdof = mfem::FiniteElementSpace::DecodeDof(encodedVDof);
|
|
MFEM_VERIFY(
|
|
vdof >= 0 && vdof < finiteElementSpace.GetVSize(), "MFEM returned an invalid point vector DOF."
|
|
);
|
|
pointVDofMarker[vdof] = 1;
|
|
}
|
|
}
|
|
|
|
finiteElementSpace.Synchronize(pointVDofMarker);
|
|
|
|
mfem::Array<int> pointReducedDofMarker(fieldDofMap.reduced_size());
|
|
pointReducedDofMarker = 0;
|
|
for (int vdof = 0; vdof < pointVDofMarker.Size(); ++vdof) {
|
|
if (pointVDofMarker[vdof] == 0) {
|
|
continue;
|
|
}
|
|
|
|
const int trueDof = finiteElementSpace.GetLocalTDofNumber(vdof);
|
|
if (trueDof < 0) {
|
|
continue;
|
|
}
|
|
|
|
const std::optional<int> reducedDof = fieldDofMap.reduced_dof(trueDof);
|
|
MFEM_VERIFY(
|
|
reducedDof.has_value(),
|
|
"A point DOF selected for the field is absent from that field's reduced solver map."
|
|
);
|
|
pointReducedDofMarker[*reducedDof] = 1;
|
|
}
|
|
|
|
mfem::Array<int> pointReducedDofs;
|
|
mfem::FiniteElementSpace::MarkerToList(pointReducedDofMarker, pointReducedDofs);
|
|
|
|
const long long localPointDofCount = pointReducedDofs.Size();
|
|
long long globalPointDofCount = 0;
|
|
MPI_Allreduce(
|
|
&localPointDofCount, &globalPointDofCount, 1, MPI_LONG_LONG, MPI_SUM, finiteElementSpace.GetComm()
|
|
);
|
|
MFEM_VERIFY(
|
|
globalPointDofCount == finiteElementSpace.GetVDim(),
|
|
"The requested geometric point must identify exactly one field vertex globally."
|
|
);
|
|
|
|
return FieldPointDofMap(fieldDofMap.reduced_size(), pointReducedDofs);
|
|
}
|
|
|
|
/*
|
|
* Canonical adapter between an MFEM GridFunction and a reduced field
|
|
* vector.
|
|
*
|
|
* FieldDofMap deliberately contains only indexing information. This
|
|
* adapter binds that indexing to the exact finite-element space whose true
|
|
* DOFs the map describes. Consequently, a grid function from another
|
|
* finite-element space is rejected even when it happens to have the same
|
|
* vector size.
|
|
*
|
|
* The finite-element space must outlive the adapter.
|
|
*/
|
|
class FieldDofGridFunctionAdapter {
|
|
public:
|
|
FieldDofGridFunctionAdapter(
|
|
FieldDofMap dofMap,
|
|
const mfem::FiniteElementSpace &finiteElementSpace
|
|
)
|
|
: m_dofMap(std::move(dofMap)),
|
|
m_finiteElementSpace(&finiteElementSpace) {
|
|
if (m_dofMap.full_size() != finiteElementSpace.GetTrueVSize()) {
|
|
throw std::invalid_argument(
|
|
"FieldDofGridFunctionAdapter map and finite-element "
|
|
"space have incompatible true-DOF sizes."
|
|
);
|
|
}
|
|
}
|
|
|
|
[[nodiscard]]
|
|
const FieldDofMap &dof_map() const noexcept {
|
|
return m_dofMap;
|
|
}
|
|
|
|
[[nodiscard]]
|
|
const mfem::FiniteElementSpace &finite_element_space() const noexcept {
|
|
return *m_finiteElementSpace;
|
|
}
|
|
|
|
/*
|
|
* Gather the grid function's true DOFs into reduced field ordering.
|
|
* The output vector is not resized so MFEM vector views remain valid.
|
|
*/
|
|
void gather(
|
|
const mfem::GridFunction &gridFunction,
|
|
mfem::Vector &reduced
|
|
) const {
|
|
validate_grid_function(gridFunction);
|
|
|
|
mfem::Vector full;
|
|
gridFunction.GetTrueDofs(full);
|
|
m_dofMap.gather(full, reduced);
|
|
}
|
|
|
|
[[nodiscard]]
|
|
mfem::Vector gather(const mfem::GridFunction &gridFunction) const {
|
|
mfem::Vector reduced(m_dofMap.reduced_size());
|
|
gather(gridFunction, reduced);
|
|
return reduced;
|
|
}
|
|
|
|
/*
|
|
* Scatter with projection semantics. Unsupported true DOFs are zeroed
|
|
* before the complete true vector is distributed to the grid function.
|
|
*/
|
|
void scatter(
|
|
const mfem::Vector &reduced,
|
|
mfem::GridFunction &gridFunction
|
|
) const {
|
|
validate_grid_function(gridFunction);
|
|
|
|
const mfem::Vector full = m_dofMap.scatter(reduced);
|
|
gridFunction.SetFromTrueDofs(full);
|
|
}
|
|
|
|
/*
|
|
* Scatter while preserving the grid function's existing unsupported
|
|
* true DOFs.
|
|
*/
|
|
void scatter_into(
|
|
const mfem::Vector &reduced,
|
|
mfem::GridFunction &gridFunction
|
|
) const {
|
|
validate_grid_function(gridFunction);
|
|
|
|
mfem::Vector full;
|
|
gridFunction.GetTrueDofs(full);
|
|
m_dofMap.scatter_into(reduced, full);
|
|
gridFunction.SetFromTrueDofs(full);
|
|
}
|
|
|
|
private:
|
|
void validate_grid_function(const mfem::GridFunction &gridFunction) const {
|
|
if (gridFunction.FESpace() != m_finiteElementSpace) {
|
|
throw std::invalid_argument(
|
|
"FieldDofGridFunctionAdapter received a grid function "
|
|
"from a different finite-element space."
|
|
);
|
|
}
|
|
}
|
|
|
|
FieldDofMap m_dofMap;
|
|
const mfem::FiniteElementSpace *m_finiteElementSpace;
|
|
};
|
|
|
|
/*
|
|
* Construct the canonical solver map for a registered spatial field.
|
|
*
|
|
* Field and domain semantics are used only while constructing the map.
|
|
* Consumers receive a plain FieldDofMap and therefore do not need to
|
|
* understand DomainSchema or field-support types.
|
|
*/
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
FieldDofMap make_field_dof_map(const mfem::ParFiniteElementSpace &finiteElementSpace) {
|
|
const FieldDofSupport support = resolve_field_dof_support<FieldT, SchemaT>(finiteElementSpace);
|
|
|
|
return FieldDofMap(support);
|
|
}
|
|
|
|
template <
|
|
MfemDomainField FieldT,
|
|
utils::domain::IsSchema SchemaT>
|
|
[[nodiscard]]
|
|
FieldDofGridFunctionAdapter
|
|
make_field_dof_grid_function_adapter(const mfem::ParFiniteElementSpace &finiteElementSpace) {
|
|
return FieldDofGridFunctionAdapter(make_field_dof_map<FieldT, SchemaT>(finiteElementSpace), finiteElementSpace);
|
|
}
|
|
} // namespace mean_field::field
|