tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Progress on bringing the DomainMapper over

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-20 16:33:06 -04:00
parent d7e7a3a015
commit 806268ccb8
9 changed files with 467 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

127
src/include/serif/mapping/data.hpp Normal file
View File

@@ -0,0 +1,127 @@
#pragma once
#include "mfem/fem/fe/fe_base.hpp"
#include "mfem/fem/fespace.hpp"
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
#include "mfem/mesh/mesh.hpp"
#include <cstdint>
#include <unordered_map>
#include <tuple>
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/mapping/context.hpp"
#include "serif/mapping/data.hpp"
namespace serif::mapping {
enum class FaceElementSide : uint8_t {
ELEMENT_1,
ELEMENT_2
};
class ElementDisplacementData {
public:
ElementDisplacementData(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs, mfem::Ordering::Type ordering = mfem::Ordering::byNODES);
[[nodiscard]] const mfem::FiniteElement& get_element() const noexcept;
[[nodiscard]] const mfem::DenseMatrix &get_dof_matrix() const noexcept;
[[nodiscard]] int get_dimension() const noexcept;
[[nodiscard]] int get_dof_count() const noexcept;
[[nodiscard]] mfem::Ordering::Type get_ordering() const noexcept;
private:
const mfem::FiniteElement *m_element;
mfem::DenseMatrix m_dof_matrix;
int m_dimension;
mfem::Ordering::Type m_ordering;
};
// Maybe this should be moved into the compactification module
struct CompactificationPointData {
double compactification_coordinate{0.0};
mfem::Vector coordinate_gradient;
};
[[nodiscard]] ElementDisplacementData element_displacement_data_from_element_vdofs(
const mfem::FiniteElement& element,
mfem::Vector &displacement_dofs
);
class ElementCompactificationData {
public:
ElementCompactificationData(const mfem::FiniteElement& element, const mfem::Vector& dofs);
[[nodiscard]] const mfem::FiniteElement& get_element() const noexcept;
[[nodiscard]] const mfem::Vector &get_dofs() const noexcept;
[[nodiscard]] int get_dof_count() const noexcept;
private:
const mfem::FiniteElement *m_element;
mfem::Vector m_dofs;
};
struct ElementMappingData {
const ElementDisplacementData &displacement;
const ElementCompactificationData &compactification;
};
class DMWorkspace {
public:
explicit DMWorkspace (int dimension = 3);
void set_dimension(int dimension);
[[nodiscard]] int get_dimension() const noexcept;
public: // This is a workspace class, these members are kept public so we do not need indirection through getters to access them
mfem::Vector m_shape;
mfem::DenseMatrix m_reference_shape;
mfem::DenseMatrix m_mesh_dshape;
mfem::DenseMatrix m_reference_field_jacobian;
mfem::Vector m_field_value;
mfem::DenseMatrix m_field_jacobian;
mfem::Vector m_compactification_shape;
mfem::DenseMatrix m_compactification_dshape;
CompactificationPointData m_compactification_point;
mfem::Vector m_reference_normal;
mfem::Vector m_mapped_normal;
mfem::DenseMatrix m_full_element_jacobian;
mfem::Vector m_vector_temp;
mfem::DenseMatrix m_matrix_temp_1;
mfem::DenseMatrix m_matrix_temp_2;
compactification::ExteriorMapResult m_exterior_result;
compactification::ExteriorMapVariation m_exterior_variation;
private:
int m_dimension;
};
struct DMFixedMeshGeometry {
mfem::Vector position;
mfem::DenseMatrix jacobian;
mfem::DenseMatrix inverse_jacobian;
double weight{0.0};
};
struct DMGeometryCacheStatistics {
std::uint64_t hits{0};
std::uint64_t misses{0};
std::size_t points{0};
};
struct DMReferenceDisplacementSample {
const mfem::Vector& value;
const mfem::DenseMatrix& gradient;
};
struct DMGeometryCache {
bool enabled{true};
const mfem::Mesh* mesh{nullptr};
int32_t sequence{-1};
int32_t nodesSequence{-1};
std::unordered_map<std::tuple<int, double, double, double>, DMFixedMeshGeometry> points;
uint64_t hits{0};
uint64_t misses{0};
};
}

161
src/include/serif/mapping/domain_mapper.hpp Normal file
View File

@@ -0,0 +1,161 @@
#pragma once
#include <cstdint>
#include <memory>
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
#include "mfem/fem/eltrans.hpp"
#include "mfem/fem/intrules.hpp"
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/mapping/context.hpp"
#include "serif/mapping/data.hpp"
#include "serif/utils/user_options/mapping/domain_mapper.hpp"
namespace serif::mapping {
class DomainMapper {
public:
DomainMapper(
utils::user_options::DomainMapperOptions options,
std::unique_ptr<compactification::ExteriorDomainMap> exterior_map
);
DomainMapper(const DomainMapper&) = delete;
DomainMapper& operator=(const DomainMapper&) = delete;
~DomainMapper();
DomainMapper(DomainMapper&&) = delete;
DomainMapper& operator=(DomainMapper&&) = delete;
void set_geometry_cache_enabled(bool enabled) const;
[[nodiscard]] bool is_geometry_cache_enabled() const noexcept;
void clear_geometry_cache() const;
[[nodiscard]] DMGeometryCacheStatistics get_geometry_cache_statistics() const;
[[nodiscard]] const DMFixedMeshGeometry* get_fixed_mesh_geometry(mfem::ElementTransformation& transformation, const mfem::IntegrationPoint& ip) const;
[[nodiscard]] MappingStatus evaluate_point(
const ElementMappingData& element_data,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
MappingPointContext& context
) const;
[[nodiscard]] MappingStatus evaluate_volume(
const ElementMappingData& element_data,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
VolumeMappingContext& context,
const DMReferenceDisplacementSample *displacement_sample = nullptr
) const;
[[nodiscard]] MappingStatus evaluate_face(
const ElementMappingData& element_data,
mfem::ElementTransformation& transformation,
FaceElementSide side,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
FaceMappingContext& context
) const;
[[nodiscard]] MappingStatus evaluate_point_variation(
const ElementMappingData& element_data,
const ElementDisplacementData& direction,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
const MappingPointContext& context,
DMWorkspace& workspace,
MappingPointVariation& variation
);
[[nodiscard]] MappingStatus evaluate_volume_variation(
const ElementMappingData& element_data,
const ElementDisplacementData& direction,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
const VolumeMappingContext& context,
DMWorkspace& workspace,
VolumeMappingVariation& variation
);
[[nodiscard]] MappingStatus evaluate_face_variation(
const ElementMappingData& element_data,
const ElementDisplacementData& direction,
mfem::ElementTransformation& transformation,
FaceElementSide side,
const mfem::IntegrationPoint& ip,
const FaceMappingContext& context,
DMWorkspace& workspace,
FaceMappingVariation& variation
);
[[nodiscard]] bool is_compactified_element(const mfem::ElementTransformation& transformation) const noexcept;
[[nodiscard]] int get_dimension() const noexcept;
[[nodiscard]] const compactification::ExteriorDomainMap& get_exterior_map() const noexcept;
private:
[[nodiscard]] MappingStatus evaluate_point_impl(
const ElementMappingData& element_data,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
MappingPointContext& context,
const DMFixedMeshGeometry* geometry,
const DMReferenceDisplacementSample* displacement_sample = nullptr
) const;
void validate_element_data(
const ElementMappingData& element_data
) const;
void evaluate_field(
const ElementDisplacementData& field,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
mfem::Vector& value,
mfem::DenseMatrix& jacobian,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const;
[[nodiscard]] MappingStatus evaluate_compactification_coordinate(
const ElementCompactificationData& compactification,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
CompactificationPointData& point_data,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const;
[[nodiscard]] MappingStatus evaluate_point_variation_impl(
const ElementMappingData& element_data,
const ElementDisplacementData& direction,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
const MappingPointContext& context,
DMWorkspace& workspace,
MappingPointVariation& variation,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const;
// These are static since they do not use any class members or workspace. Really these could and
// maybe should be free functions...
[[nodiscard]] static mfem::ElementTransformation& select_face_element_transformation(
mfem::FaceElementTransformations& transformations,
FaceElementSide side
);
[[nodiscard]] static const mfem::IntegrationPoint& select_face_element_integration_point(
mfem::FaceElementTransformations& transformations,
FaceElementSide side
);
private:
mutable std::unique_ptr<DMGeometryCache> m_geometryCache;
utils::user_options::DomainMapperOptions m_options;
std::unique_ptr<const compactification::ExteriorDomainMap> m_exterior_map;
};
}

6
src/include/serif/utils/error/codes.hpp
View File

@@ -51,6 +51,12 @@ namespace serif::utils::errors {
NON_FINITE_MAPPING_RESULT,
NON_POSITIVE_DETERMINATE,
INVALID_DOF_COUNT,
INVALID_DOF_ORDERING,
INVALID_COMPACTIFICATION_RANK,
INVALID_COMPACTIFICATION_TYPE,
INVALID_COMPACTIFICATION_DERIVATIVE_VALUE,
IMPROPERLY_ORDERED_REFERENCE_RADII = 450,
INVALID_COORDINATE_TOLERANCE,

8
src/include/serif/utils/user_options/mapping/domain_mapper.hpp Normal file
View File

@@ -0,0 +1,8 @@
#pragma once
namespace serif::utils::user_options {
struct DomainMapperOptions {
int dimension{3};
int vacuum_element_attribute_ID{3};
};
}

153
src/lib/serif/mapping/data.cpp Normal file
View File

@@ -0,0 +1,153 @@
#include <stdexcept>
#include "mfem/linalg/vector.hpp"
#include "mfem/fem/fe/fe_base.hpp"
#include "mfem/fem/fespace.hpp"
#include "serif/mapping/data.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping {
using namespace serif::utils::errors;
using mfem::Ordering;
ElementDisplacementData::ElementDisplacementData(
const mfem::FiniteElement& element,
mfem::Vector& displacement_dofs,
const Ordering::Type ordering
) : m_element(&element), m_dimension(0), m_ordering(ordering) {
const int dof_count = element.GetDof();
if (dof_count <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("ElementDisplacementData: Displacement elements must have at least one degree of freedom. Got an element with {} dofs", dof_count);
}
if (displacement_dofs.Size() <= 0 || displacement_dofs.Size() % dof_count != 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("ElementDisplacementData: The displacement vector size ({}) must be a positive multiple of the element degree of freedom count ({}).", displacement_dofs.Size(), dof_count);
}
m_dimension = displacement_dofs.Size() / dof_count;
m_dof_matrix.SetSize(dof_count, m_dimension);
switch (ordering) {
case Ordering::byNODES: {
for (int componentID = 0; componentID < m_dimension; ++componentID) {
for (int dofID = 0; dofID < dof_count; ++dofID) {
m_dof_matrix(dofID, componentID) = displacement_dofs(dofID + componentID * dof_count);
}
}
break;
}
case Ordering::byVDIM: {
for (int dofID = 0; dofID < dof_count; ++dofID) {
for (int componentID = 0; componentID < m_dimension; ++componentID) {
m_dof_matrix(dofID, componentID) = displacement_dofs(componentID + dofID * m_dimension);
}
}
break;
}
default:
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_ORDERING>("Unsupported ordering");
}
}
const mfem::FiniteElement& ElementDisplacementData::get_element() const noexcept {
return *m_element;
}
const mfem::DenseMatrix& ElementDisplacementData::get_dof_matrix() const noexcept {
return m_dof_matrix;
}
int ElementDisplacementData::get_dimension() const noexcept {
return m_dimension;
}
int ElementDisplacementData::get_dof_count() const noexcept {
return m_element->GetDof();
}
Ordering::Type ElementDisplacementData::get_ordering() const noexcept {
return m_ordering;
}
ElementDisplacementData element_displacement_data_from_element_vdofs(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs) {
return ElementDisplacementData(element, displacement_dofs, Ordering::byNODES);
}
ElementCompactificationData::ElementCompactificationData(
const mfem::FiniteElement& element,
const mfem::Vector& dofs
) : m_element(&element), m_dofs(dofs) {
if (element.GetRangeType() != mfem::FiniteElement::SCALAR) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>("Compactification coordinate requires a scalar finite element.");
}
if (element.GetMapType() != mfem::FiniteElement::VALUE) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_TYPE>("Compactification coordinate requires a value-mapped scalar finite element");
}
if (element.GetDerivType() != mfem::FiniteElement::GRAD) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_VALUE>("Compactification coordinate finite element must provide a gradient");
}
if (element.GetDof() <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Compactification coordinatge finite element must have a positive non-zero number of degrees of freedom. Recived {}", element.GetDof());
}
if (dofs.Size() != element.GetDof()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Compactification coordinate degree of freedom count ({}) does not match its finite element dof count ({})", dofs.Size(), element.GetDof());
}
}
const mfem::FiniteElement& ElementCompactificationData::get_element() const noexcept {
return *m_element;
}
const mfem::Vector& ElementCompactificationData::get_dofs() const noexcept {
return m_dofs;
}
int ElementCompactificationData::get_dof_count() const noexcept {
return m_dofs.Size();
}
DMWorkspace::DMWorkspace(const int dimension) {
set_dimension(dimension);
}
void DMWorkspace::set_dimension(const int dimension) {
if (dimension <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Domain mapping workspace dimension must be positive and finite. Received {}", dimension);
}
m_dimension = dimension;
m_field_value.SetSize(dimension);
m_field_jacobian.SetSize(dimension, dimension);
m_reference_field_jacobian.SetSize(dimension, dimension);
m_compactification_point.compactification_coordinate = 0.0;
m_compactification_point.coordinate_gradient.SetSize(dimension);
m_reference_normal.SetSize(dimension);
m_mapped_normal.SetSize(dimension);
m_full_element_jacobian.SetSize(dimension, dimension);
m_vector_temp.SetSize(dimension);
m_matrix_temp_1.SetSize(dimension, dimension);
m_matrix_temp_2.SetSize(dimension, dimension);
m_exterior_result.physical_position.SetSize(dimension);
m_exterior_result.mapping_jacobian.SetSize(dimension, dimension);
m_exterior_variation.physical_position_variation.SetSize(dimension);
m_exterior_variation.mapping_jacobian_variation.SetSize(dimension, dimension);
}
int DMWorkspace::get_dimension() const noexcept {
return m_dimension;
}
}

1
src/meson.build
View File

@@ -6,6 +6,7 @@ serif_sources = files(
'lib/serif/discretization/quadrature/policy.cpp',
'lib/serif/eos/models/polytropic.cpp',
'lib/serif/mapping/compactification/kelvin.cpp',
'lib/serif/mapping/data.cpp',
)
serif_lib = library(

7
subprojects/packagefiles/stroid/build-check/meson.build Normal file
View File

@@ -0,0 +1,7 @@
# A subproject must not add global arguments after the parent defines targets.
cpp = meson.get_compiler('cpp')
add_project_arguments(
cpp.get_supported_arguments('-Wno-unused-command-line-argument'),
language: 'cpp',
)
subdir('CPPC')

3
subprojects/packagefiles/stroid/build-config/mfem/meson.build Normal file
View File

@@ -0,0 +1,3 @@
# The parent overrides this dependency with its selected MFEM build. Building
# another copy here can change class layouts and cause crashes across the ABI.
mfem_dep = dependency('mfem', version: '>=4.10')

1
tests/sandbox/mapping_sandbox.cpp
View File

@@ -5,6 +5,7 @@
#include "serif/mapping/compactification/compactification.hpp"
#include "serif/mapping/status.hpp"
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/mapping/domain_mapper.hpp"
#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
int main() {