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

Merge pull request #2 from tboudreaux/feat/C2

Browse Source 
Added TMOP solver to improve mesh conditioning
This commit is contained in:
Emily Boudreaux
2026-07-01 11:15:32 -04:00
committed by GitHub
parent a5ddf6a62f 9aaa8529e0
commit 7f69f19273
59 changed files with 3082 additions and 216 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Doxyfile
View File

@@ -48,7 +48,7 @@ PROJECT_NAME = stroid
# could be handy for archiving the generated documentation or if some version
# control system is used.
PROJECT_NUMBER = v0.2.1
PROJECT_NUMBER = v0.3.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewers a

16
build-config/meson.build
View File

@@ -2,3 +2,19 @@ subdir('mfem')
subdir('libconfig')
subdir('CLI11')
subdir('magic_enum')
if get_option('build_python')
subdir('python')
subdir('pybind')
endif
if get_option('build_python')
stroid_pkg_dir = py_installation.get_install_dir() / 'stroid'
stroid_includedir = stroid_pkg_dir / 'include'
stroid_libdir = stroid_pkg_dir / 'lib'
stroid_pcdir = stroid_libdir / 'pkgconfig'
else
stroid_includedir = get_option('includedir')
stroid_libdir = get_option('libdir')
stroid_pcdir = get_option('libdir') / 'pkgconfig'
endif

3
build-config/pybind/meson.build Normal file
View File

@@ -0,0 +1,3 @@
pybind11_proj = subproject('pybind11')
pybind11_dep = pybind11_proj.get_variable('pybind11_dep')
python3_dep = dependency('python3')

5
build-config/python/meson.build Normal file
View File

@@ -0,0 +1,5 @@
py_installation = import('python').find_installation('python3', pure: false)
py_dep = py_installation.dependency()
py_module_prefix = ''
py_module_suffix = 'so'

43
build-python/meson.build Normal file
View File

@@ -0,0 +1,43 @@
if get_option('build_python')
message('Building Python bindings...')
stroid_py_deps = [
py_dep,
pybind11_dep,
stroid_dep
]
if host_machine.system() == 'darwin'
stroid_ext_rpath = '@loader_path/lib'
else
stroid_ext_rpath = '$ORIGIN/lib'
endif
py_sources = [
meson.project_source_root() + '/src/python/bindings.cpp',
meson.project_source_root() + '/src/python/config/bindings.cpp',
meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
meson.project_source_root() + '/src/python/IO/bindings.cpp',
meson.project_source_root() + '/src/python/refinement/bindings.cpp',
meson.project_source_root() + '/src/python/utils/bindings.cpp',
]
py_mod = py_installation.extension_module(
'_stroid',
sources: py_sources,
dependencies: stroid_py_deps,
install: true,
link_args: stroid_ext_rpath_args,
build_rpath: stroid_ext_rpath,
install_rpath: stroid_ext_rpath,
subdir: 'stroid',
)
py_installation.install_sources(
meson.project_source_root() + '/src/python/stroid/__init__.py',
subdir: 'stroid',
)
else
message('Python bindings disabled')
endif

3
configs/default_config.toml
View File

@@ -13,3 +13,6 @@ surface_bdr_id = 1
core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_external_domain.toml
View File

@@ -13,3 +13,5 @@ surface_bdr_id = 1
core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_external_domain_refinement_l1.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_flattening.toml
View File

@@ -13,3 +13,5 @@ surface_bdr_id = 1
core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

17
configs/test_polynomial_projection.toml Normal file
View File

@@ -0,0 +1,17 @@
[main]
core_steepness = 1.0
flattening = 0.2
include_external_domain = false
inf_bdr_id = 2
order = 3
r_core = 1.5
r_infinity = 6.0
r_instability = 1e-14
r_star = 5.0
refinement_levels = 1
surface_bdr_id = 1
core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_refinement_l2.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_volume_no_external.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_volume_spherical_no_external.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_volume_spherical_with_external.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

2
configs/test_volume_with_external.toml
View File

@@ -14,3 +14,5 @@ core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
smoothstep = true

8
meson.build
View File

@@ -1,4 +1,4 @@
project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.2.1', default_options : ['cpp_std=c++23'])
project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.4.0', default_options : ['cpp_std=c++23'])
subdir('build-check')
@@ -13,6 +13,10 @@ if get_option('build_tools')
subdir('tools')
endif
if get_option('build_python')
subdir('build-python')
endif
if get_option('pkg_config')
pkg = import('pkgconfig')
pkg.generate(
@@ -20,7 +24,7 @@ if get_option('pkg_config')
description: 'Stroid multi-block curvilinear mesh generation library',
version: meson.project_version(),
libraries: [
stroid_lib
libstroid
],
subdirs: ['stroid'],
filebase: 'stroid',

1
meson_options.txt
View File

@@ -1,3 +1,4 @@
option('pkg_config', type: 'boolean', value: false, description: 'generate pkg-config file for stroid')
option('build_tests', type: 'boolean', value: true, description: 'compile subproject tests')
option('build_tools', type: 'boolean', value: true, description: 'compile stroid command line tools')
option('build_python', type: 'boolean', value: true, description: 'compile stroid python bindings')

25
pyproject.toml Normal file
View File

@@ -0,0 +1,25 @@
[build-system]
requires = ["meson-python>=0.19.0", "meson>=1.9.1", "pybind11==3.0.0", "fourdst==0.10.6"]
build-backend = "mesonpy"
[project]
name = "stroid"
dynamic = ["version"]
description = "O-grid mesh generation with multiple domains"
readme = "README.md"
license = { file = "LICENSE.txt" }
authors = [
{name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"},
]
maintainers = [
{name = "Emily M. Boudreaux", email = "emily@boudreauxmail.com"}
]
[tool.meson-python.args]
setup = [
'-Dbuild_tools=false',
'-Dbuild_tests=false',
'-Dpkg_config=false'
]
install = ['--skip-subprojects']

9
readme.md
View File

@@ -92,11 +92,15 @@ inf_bdr_id = 2
core_id = 1
envelope_id = 2
vacuum_id = 3
[main.optimization_methods]
tmop = false
smoothstep = true
```
<!-- Table of what these parameters do -->
| Parameter | Description | Default |
|-------------------------|-----------------------------------------------------------------------------------------------------|---------|
|---------------------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------|
| refinement_levels | Number of uniform refinement levels to apply to the mesh after generation | 4 |
| order | The polynomial order of the finite elements in the mesh | 3 |
| include_external_domain | Whether to include an external domain extending to r_infinity | true |
@@ -111,6 +115,8 @@ vacuum_id = 3
| core_id | The material ID to assign to the core region of the star | 1 |
| envelope_id | The material ID to assign to the envelope region of the star | 2 |
| vacuum_id | The material ID to assign to the vacuum region of the star (if included) | 3 |
| optimization_methods.tmop | The tmop flag enables or disables the use of TMOP ideal shape unit size metric optimization during mesh generation. This can help improve the quality of the generated mesh, but will dramatically increase the time required for mesh generation. | false |
| optimization_methods.smoothstep | The smoothstep flag enables or disables the use of a smoothstep function to transition between the core and envelope regions of the star. This can help improve the quality of the generated mesh | true |
If no configuration file is provided, stroid will use the default parameters listed above. Further, configuration files
@@ -138,6 +144,7 @@ int main() {
stroid::topology::Finalize(*mesh, cfg);
stroid::topology::PromoteToHighOrder(*mesh, cfg);
stroid::topology::ProjectMesh(*mesh, cfg);
stroid::topology::OptimizeMesh(*mesh, cfg);
stroid::IO::ViewMesh(*mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID);

45
src/include/stroid/IO/mesh.h
View File

@@ -1,7 +1,12 @@
#pragma once
#include <string>
#include <expected>
#include <istream>
#include "mfem.hpp"
#include "stroid/utils/types.h"
namespace stroid::IO {
/**
* @brief Visualization modes for GLVis display.
@@ -15,18 +20,43 @@ namespace stroid::IO {
BOUNDARY_ELEMENT_ID
};
void SaveStroidMesh(const StroidMesh& mesh, const std::string& filename, const std::string& comment="");
/**
* @brief Save a mesh to MFEM's native `.mesh` format.
* @param mesh Mesh to serialize.
* @param filename Output path (including extension).
*/
void SaveMesh(const mfem::Mesh& mesh, const std::string& filename);
/**
* @brief Overload of SaveMesh which accepts a StroidMesh type and will internally unpack it
* @param mesh StroidMesh to serialize.
* @param filename Path to save to
*
* @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer
* in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an
* entire StroidMesh then please use the stroid::IO::SaveStroidMesh function
*/
void SaveMesh(const stroid::StroidMesh& mesh, const std::string& filename);
/**
* @brief Save a mesh as a ParaView VTU dataset.
* @param mesh Mesh to export.
* @param exportName Output base name (ParaView will add extensions).
*/
void SaveVTU(mfem::Mesh& mesh, const std::string& exportName);
/**
* @brief Overload of SaveVTU which accepts a StroidMesh type and will internally unpack it
* @param mesh StroidMesh to serialize.
* @param filename Path to save to
*
* @note This function is a utility wrapper to save a StroidMesh object in MFEM's native .mesh format. Data other than the mesh pointer
* in StroidMesh **will not be saved** (e.g. the reference mesh, the number of refinement levels, etc..). If you need to serialize an
* entire StroidMesh then please use the stroid::IO::SaveStroidVTU function
*/
void SaveVTU(const stroid::StroidMesh& mesh, const std::string& exportName);
/**
* @brief Stream a mesh to a running GLVis server for interactive viewing.
* @param mesh Mesh to display.
@@ -36,9 +66,22 @@ namespace stroid::IO {
* @param visport GLVis server port.
*/
void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport);
void ViewMesh(const stroid::StroidMesh& mesh, const std::string& title, VISUALIZATION_MODE mode, const std::string &vishost, int visport);
/**
* @brief Visualize boundary face valence (1=surface, 2=internal).
* @param mesh Mesh whose boundary faces are inspected.
*/
void VisualizeFaceValence(mfem::Mesh& mesh);
void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport);
void VisualizeFaceValence(const stroid::StroidMesh& mesh, const std::string &vishost, int visport);
std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is);
std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename);
#ifdef MFEM_USE_MPI
std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is, MPI_Comm comm);
std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename, MPI_Comm comm);
#endif
}

48
src/include/stroid/config/config.h
View File

@@ -1,6 +1,16 @@
#pragma once
#include <cstdint>
#include <optional>
#include <string>
namespace stroid::config {
struct OptimizationMethods {
std::optional<bool> tmop{false};
std::optional<bool> smoothstep{true};
};
/**
* @brief Configuration parameters for stroid mesh generation.
*
@@ -15,92 +25,102 @@ namespace stroid::config {
* @section toml
* - [main].refinement_levels
*/
int refinement_levels = 4;
std::optional<int> refinement_levels = 4;
/**
* @brief Polynomial order for high-order elements.
* @section toml
* - [main].order
*/
int order = 3;
std::optional<int> order = 3;
/**
* @brief Whether to include an external domain extending to `r_infinity`.
* @section toml
* - [main].include_external_domain
*/
bool include_external_domain = true;
std::optional<bool> include_external_domain = true;
/**
* @brief Radius of the stellar core region.
* @section toml
* - [main].r_core
*/
double r_core = 1.5;
std::optional<double> r_core = 0.25;
/**
* @brief Radius of the stellar surface.
* @section toml
* - [main].r_star
*/
double r_star = 5.0;
std::optional<double> r_star = 1.0;
/**
* @brief Flattening factor for spheroidal shaping (0 = spherical, >0 = oblate).
* @section toml
* - [main].flattening
*/
double flattening = 0;
std::optional<double> flattening = 0;
/**
* @brief Outer radius of the external domain when enabled.
* @section toml
* - [main].r_infinity
*/
double r_infinity = 6.0;
std::optional<double> r_infinity = 6.0;
/**
* @brief Radius inside which transformations are skipped to avoid singularities.
* @section toml
* - [main].r_instability
*/
double r_instability = 1e-14;
std::optional<double> r_instability = 1e-14;
/**
* @brief Controls the smoothness/steepness of the core-to-envelope transition.
* @section toml
* - [main].core_steepness
*/
double core_steepness = 1.0;
std::optional<double> core_steepness = 1.0;
/**
* @brief Continuity order for the core-envelope transition (0 = discontinuous, 1 = C1, 2 = C2).
* @section toml
* - [main].continuity_order
*/
std::optional<size_t> continuity_order = 2;
/**
* @brief Boundary attribute id for stellar surface
* @section toml
* - [main].surface_bdr_id
*/
size_t surface_bdr_id = 1;
std::optional<size_t> surface_bdr_id = 1;
/**
* @brief Boundary attribute id for infinity in kelvin mapping
* @section toml
* - [main].inf_bdr_id
*/
size_t inf_bdr_id = 2;
std::optional<size_t> inf_bdr_id = 2;
/**
* @brief Material attribute id for the core region
* @section toml
* - [main].core_id
*/
size_t core_id = 1;
std::optional<size_t> core_id = 1;
/**
* @brief Material attribute id for the envelope region
* @section toml
* - [main].envelope_id
*/
size_t envelope_id = 2;
std::optional<size_t> envelope_id = 2;
/**
* @brief Material attribute id for the external domain (if enabled)
* @section toml
* - [main].vacuum_id
*/
size_t vacuum_id = 3;
std::optional<size_t> vacuum_id = 3;
std::optional<OptimizationMethods> optimization_methods = OptimizationMethods{true, true};
};
}

3
src/include/stroid/exceptions/exceptions.h Normal file
View File

@@ -0,0 +1,3 @@
#pragma once
#include "stroid/exceptions/stroid_error.h"

25
src/include/stroid/exceptions/stroid_error.h Normal file
View File

@@ -0,0 +1,25 @@
#pragma once
#include <exception>
#include <string>
namespace stroid::exceptions {
class StroidError : public std::exception {
public:
explicit StroidError(std::string message) : m_msg(std::move(message)) {}
const char* what() const noexcept override { return m_msg.c_str(); }
private:
std::string m_msg;
};
class StroidContinuityError : public StroidError {
using StroidError::StroidError;
};
class StroidMeshError : public StroidError {
using StroidError::StroidError;
};
class StroidMissingReferenceMesh : public StroidMeshError {
using StroidMeshError::StroidMeshError;
};
}

4
src/include/stroid/meson.build
View File

@@ -21,8 +21,8 @@ config.set('STROID_VERSION_PATCH', ver_parts[2])
config.set('STROID_VERSION_TAG', ver_parts[3])
configure_file(
input : 'stroid.h.in',
output : 'stroid.h',
input : 'version.h.in',
output : 'version.h',
configuration : config ,
install: true,
install_dir: get_option('includedir') / 'stroid'

7
src/include/stroid/refinement/uniform.h Normal file
View File

@@ -0,0 +1,7 @@
#pragma once
#include "stroid/utils/types.h"
namespace stroid::refinement {
void UniformRefinement(StroidMesh& mesh, size_t levels);
}

62
src/include/stroid/stroid.h.in → src/include/stroid/stroid.h
View File

@@ -4,8 +4,12 @@
#include "stroid/topology/topology.h"
#include "stroid/topology/mapping.h"
#include "stroid/topology/curvilinear.h"
#include "stroid/topology/optimize.h"
#include "stroid/utils/mesh_utils.h"
#include "stroid/IO/mesh.h"
#include "stroid/utils/types.h"
#include "stroid/refinement/uniform.h"
#include "stroid/utils/mesh_stats.h"
/**
* @namespace stroid
@@ -44,46 +48,36 @@
* @endcode
*/
namespace stroid {
/**
* @brief Version helpers for the stroid library.
*/
struct version {
static constexpr int major = @STROID_VERSION_MAJOR@;
static constexpr int minor = @STROID_VERSION_MINOR@;
static constexpr int patch = @STROID_VERSION_PATCH@;
static constexpr const char* tag = "@STROID_VERSION_TAG@";
inline StroidMesh GenerateMesh(const fourdst::config::Config<stroid::config::MeshConfig>& cfg) {
StroidMesh sm;
sm.config = *cfg;
auto reference = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*reference, cfg);
sm.refinement_levels = cfg->refinement_levels.value_or(0);
static std::string toString() {
std::string versionStr = std::to_string(major) + "." +
std::to_string(minor) + "." +
std::to_string(patch);
if (std::string(tag) != "") {
versionStr += "-" + std::string(tag);
sm.reference_mesh = std::move(reference);
sm.mesh = utils::BuildProjected(*sm.reference_mesh, cfg);
if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) {
stroid::topology::ApplyTMOP(*sm.mesh, cfg);
}
return versionStr;
}
friend std::ostream& operator<<(std::ostream& os, const version&) {
os << toString();
return os;
return sm;
}
inline StroidMesh GenerateMesh(const stroid::config::MeshConfig& config) {
fourdst::config::Config<config::MeshConfig> cfg;
auto Mutator = [&config](config::MeshConfig& orig) {
orig = config;
};
cfg.mutate(Mutator);
return GenerateMesh(cfg);
}
inline StroidMesh GenerateMesh(const std::string& filename) {
fourdst::config::Config<stroid::config::MeshConfig> config;
config.load(filename);
return GenerateMesh(config);
}
}
/**
* @namespace std
* @brief Standard library extensions used by stroid.
*
* Provides a `std::formatter` specialization for `stroid::version` so it can
* be used with `std::format` and related APIs.
*/
// Overload format struct
template <>
struct std::formatter<stroid::version> : std::formatter<std::string> {
auto format(const stroid::version& v, auto& ctx) {
return std::formatter<std::string>::format(stroid::version::toString(), ctx);
}
};
/**
* @namespace stroid::config

19
src/include/stroid/topology/optimize.h Normal file
View File

@@ -0,0 +1,19 @@
#pragma once
#include "mfem.hpp"
#include "fourdst/config/base.h"
#include "stroid/config/config.h"
namespace stroid::topology {
/**
* @breif Apply target matrix optimization to improve conditioning of the mesh
*/
void ApplyTMOP(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &config);
/**
*@breif Helper to call TMOP if the correct flags are set
*/
void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &cfg);
}

174
src/include/stroid/utils/mesh_stats.h Normal file
View File

@@ -0,0 +1,174 @@
#pragma once
#include "mfem.hpp"
#include "stroid/utils/types.h"
#include "stroid/config/config.h"
#include <cstdint>
#include <optional>
#include <string>
#include <vector>
namespace stroid::stats {
enum class MeshStatFeatures : uint32_t {
NONE = 0u,
RADIUS = 1u << 0,
AXES = 1u << 1,
ELLIPTICITY = 1u << 2,
BOWING = 1u << 3,
CONFORMITY = 1u << 4,
JACOBIAN = 1u << 5,
VOLUME_AREA = 1u << 6,
ELEMENT_COUNT = 1u << 7,
MESH_SIZE = 1u << 8,
OUTER_BOUNDS = 1u << 9,
CENTROID = 1u << 10,
CONFIG_META = 1u << 11,
BOUNDING_BOX = 1u << 12,
};
constexpr MeshStatFeatures operator|(MeshStatFeatures lhs, MeshStatFeatures rhs) {
return static_cast<MeshStatFeatures>(static_cast<uint32_t>(lhs) | static_cast<uint32_t>(rhs));
}
constexpr MeshStatFeatures operator&(MeshStatFeatures lhs, MeshStatFeatures rhs) {
return static_cast<MeshStatFeatures>(static_cast<uint32_t>(lhs) & static_cast<uint32_t>(rhs));
}
constexpr bool has_feature(MeshStatFeatures feature, MeshStatFeatures set) {
return (static_cast<uint32_t>(set) & static_cast<uint32_t>(feature)) != 0u;
}
inline constexpr MeshStatFeatures MESH_STAT_DEFAULT =
MeshStatFeatures::RADIUS | MeshStatFeatures::AXES | MeshStatFeatures::ELLIPTICITY |
MeshStatFeatures::CONFORMITY | MeshStatFeatures::CONFIG_META;
inline constexpr auto MESH_STAT_ALL = static_cast<MeshStatFeatures>(0xFFFFFFFFu);
struct RadiusStats {
double min = 0, max = 0, mean = 0, stddev = 0;
long n_samples = 0;
};
struct AxisStats {
double semi_major = 0;
double semi_minor = 0;
};
struct EllipticityStats {
double flattening = 0;
double polar_equatorial = 1;
double radius_uniformity = 1;
};
struct BowingStats {
double max_inward = 0;
double max_outward = 0;
double rms = 0;
double worst_at_radius = 0;
};
struct ConformityStats {
bool conforming = true;
long n_nonconforming_faces = 0;
};
struct JacobianStats {
double detJ_min;
double detJ_max;
long n_flipped;
double min_detJ_ratio;
double worst_ratio_at_radius;
double detJ_min_at_radius;
long n_elements;
};
struct VolumeAreaStats {
double stellar_volume = 0, surface_area = 0;
double analytic_volume = 0, analytic_area = 0;
};
struct ElementCounts {
long total = 0, core = 0, envelope = 0, vacuum = 0, other = 0;
long n_vertices = 0;
};
struct MeshSizeStats {
double h_min = 0, h_max = 0, h_mean = 0, h_stddev = 0;
};
struct OuterBoundsStats {
double min = 0, max = 0, mean = 0;
long n_samples = 0;
};
struct CentroidStats {
double x = 0, y = 0, z = 0, offset = 0;
};
struct ConfigMeta {
double r_core = 0, r_star = 0, flattening = 0, r_infinity = 0;
int geom_order = 0;
size_t refinement_levels = 0;
bool has_external_domain = true;
};
struct BoundingBox {
double xMin = 0, xMax = 0, yMin = 0, yMax = 0, zMin = 0, zMax = 0;
bool valid = false;
[[nodiscard]] double dx() const {return xMax - xMin;}
[[nodiscard]] double dy() const {return yMax - yMin;}
[[nodiscard]] double dz() const {return zMax - zMin;}
[[nodiscard]] double diag() const {
const double a = dx(), b = dy(), c = dz();
return std::sqrt(a*a + b*b + c*c);
}
};
struct BoundingBoxStats {
BoundingBox core;
BoundingBox star;
BoundingBox vacuum;
};
struct MeshStats {
MeshStatFeatures computed = MeshStatFeatures::NONE;
std::optional<RadiusStats> radius;
std::optional<AxisStats> axes;
std::optional<EllipticityStats> ellipticity;
std::optional<BowingStats> bowing;
std::optional<ConformityStats> conformity;
std::optional<JacobianStats> jacobian;
std::optional<JacobianStats> jacobian_stellar;
std::optional<JacobianStats> jacobian_vacuum;
std::optional<VolumeAreaStats> volume;
std::optional<ElementCounts> element_counts;
std::optional<MeshSizeStats> mesh_size;
std::optional<OuterBoundsStats> outer_bounds;
std::optional<CentroidStats> centroid;
std::optional<ConfigMeta> config_meta;
std::optional<BoundingBoxStats> bounding_box;
std::vector<std::string> warnings;
std::vector<std::string> errors;
};
MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features = MESH_STAT_DEFAULT, int sample_order = -1);
std::string to_string(const MeshStats& s);
inline std::ostream& operator<<(std::ostream& os, const MeshStats& s) {
return os << to_string(s);
}
}
template <>
struct std::formatter<stroid::stats::MeshStats, char> {
static constexpr auto parse(const std::format_parse_context& ctx) {
return ctx.begin();
}
static auto format(const stroid::stats::MeshStats &s, std::format_context& ctx) {
return std::format_to(ctx.out(), "{}", stroid::stats::to_string(s));
}
};

8
src/include/stroid/utils/mesh_utils.h
View File

@@ -2,6 +2,9 @@
#include "mfem.hpp"
#include "stroid/config/config.h"
#include "fourdst/config/config.h"
namespace stroid::utils {
/**
* @brief Mark elements with negative Jacobian determinant.
@@ -15,4 +18,9 @@ namespace stroid::utils {
* @param mesh Mesh to scan and update in-place.
*/
void MarkFlippedBoundaryElements(mfem::Mesh& mesh);
void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename);
std::unique_ptr<mfem::Mesh> BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config<config::MeshConfig>& cfg);
}

65
src/include/stroid/utils/types.h Normal file
View File

@@ -0,0 +1,65 @@
#pragma once
#include "mfem.hpp"
#include "stroid/config/config.h"
#include <memory>
#include <expected>
#include <string>
#include <unordered_map>
#include <variant>
namespace stroid {
enum class MFEM_MESH_TYPE {
SERIAL,
PARALLEL
};
struct StroidMesh {
MFEM_MESH_TYPE type;
std::unique_ptr<mfem::Mesh> mesh;
std::unique_ptr<mfem::Mesh> reference_mesh;
config::MeshConfig config;
size_t refinement_levels;
[[nodiscard]] std::expected<mfem::Mesh*, std::string> as_mesh() const {
if (type == MFEM_MESH_TYPE::SERIAL) {
return mesh.get();
}
return std::unexpected{"Mesh is not serial. Try calling as_par_mesh()"};
}
[[nodiscard]] std::expected<mfem::Mesh*, std::string> ref_as_mesh() const {
if (type == MFEM_MESH_TYPE::SERIAL) {
return reference_mesh.get();
}
return std::unexpected{"Reference mesh is not serial. Try calling as_par_mesh()"};
}
[[nodiscard]] std::expected<std::unordered_map<std::string, std::variant<int, double, std::string, bool>>, std::string> mesh_stats(bool use_ref_mesh = false) const {
if (type != MFEM_MESH_TYPE::SERIAL) {
return std::unexpected{"Mesh is not serial. Mesh stats currently only supports serial meshes."};
}
mfem::Mesh* umesh;
if (use_ref_mesh) {
umesh = reference_mesh.get();
} else {
umesh = mesh.get();
}
std::unordered_map<std::string, std::variant<int, double, std::string, bool>> mesh_stats;
mesh_stats.emplace("num_elements", umesh->GetNE());
mesh_stats.emplace("num_vertices", umesh->GetNV());
mesh_stats.emplace("num_edges", umesh->GetNEdges());
mesh_stats.emplace("num_faces", umesh->GetNFaces());
mesh_stats.emplace("num_boundary_elements", umesh->GetNBE());
mesh_stats.emplace("max_bdr_attribute_id", umesh->bdr_attributes.Max());
mesh_stats.emplace("min_bdr_attribute_id", umesh->bdr_attributes.Min());
mesh_stats.emplace("max_element_attribute_id", umesh->attributes.Max());
mesh_stats.emplace("min_element_attribute_id", umesh->attributes.Min());
return mesh_stats;
}
};
}

46
src/include/stroid/version.h.in Normal file
View File

@@ -0,0 +1,46 @@
#pragma once
#include <string>
#include <ostream>
namespace stroid {
/**
* @brief Version helpers for the stroid library.
*/
struct version {
static constexpr int major = @STROID_VERSION_MAJOR@;
static constexpr int minor = @STROID_VERSION_MINOR@;
static constexpr int patch = @STROID_VERSION_PATCH@;
static constexpr const char* tag = "@STROID_VERSION_TAG@";
static std::string toString() {
std::string versionStr = std::to_string(major) + "." +
std::to_string(minor) + "." +
std::to_string(patch);
if (std::string(tag) != "") {
versionStr += "-" + std::string(tag);
}
return versionStr;
}
friend std::ostream& operator<<(std::ostream& os, const version&) {
os << toString();
return os;
}
};
}
/**
* @namespace std
* @brief Standard library extensions used by stroid.
*
* Provides a `std::formatter` specialization for `stroid::version` so it can
* be used with `std::format` and related APIs.
*/
// Overload format struct
template <>
struct std::formatter<stroid::version> : std::formatter<std::string> {
auto format(const stroid::version& v, auto& ctx) {
return std::formatter<std::string>::format(stroid::version::toString(), ctx);
}
};

505
src/lib/IO/mesh.cpp
View File

@@ -2,12 +2,428 @@
#include "stroid/config/config.h"
#include "stroid/IO/mesh.h"
#include <charconv>
#include "stroid/version.h"
#include <fstream>
#include <iostream>
#include <cstdint>
#include <format>
#include <chrono>
namespace stroid::IO {
namespace {
std::string format_header(const StroidMesh& mesh, const std::string& comment) {
auto now = std::chrono::system_clock::now();
version v;
std::stringstream vs;
vs << v;
std::string header = std::format(R"(# STROID MESH
# NOTE: STROID MESH IS A THIN WRAPPER AROUND MFEM's NATIVE MESH FORMAT
# STRUCTURE:
# - Type : Serial or Parallel (S for Serial, P for Parallel)
# - mesh : the primary computational domain which can be of n order and be h-refined
# - reference mesh : a reference, linear order mesh, used to ensure that the primary mesh remains well formed
# - config : The configuration options initially used to generate the mesh
# - refinement-levels : the total number of refinement levels the primary mesh has been subjected too
# NOTE: EACH BLOCK OF DATA IS STORED BETWEEN "BEGIN BLOCK <NAME>\n ... \nEND BLOCK <NAME>
# PARSING THE UNDERLYING MFEM NATIVE MESH FORMAT CAN BE DONE WITH MFEM'S STREAM READER
# IF YOU EXTRACT THE RAW CONTENTS BETWEEN THOSE LINES
BEGIN BLOCK HEADER
MESH_TYPE:{}
REFINEMENT_LEVELS:{}
DATE_CREATED:{:%Y-%m-%d}
COMMENT:{}
STROID_VERSION:{}
END BLOCK HEADER)",
mesh.type == MFEM_MESH_TYPE::PARALLEL ? "P" : "S",
mesh.refinement_levels,
now,
comment,
vs.str(),
mesh.refinement_levels
);
return header;
}
std::string format_primary_mesh(const StroidMesh& mesh) {
std::stringstream ss;
ss.precision(8);
mesh.mesh->Print(ss);
std::string pmesh = std::format("BEGIN BLOCK PMESH\n{}END BLOCK PMESH", ss.str());
return pmesh;
}
template <typename T>
std::string format_opt(const std::optional<T> opt, T default_val) {
if (opt.has_value()) {
return std::format("{}", opt.value());
}
return std::format("{}", default_val);
}
std::string format_reference_mesh(const StroidMesh& mesh) {
std::stringstream ss;
ss.precision(8);
mesh.reference_mesh->Print(ss);
std::string rmesh = std::format("BEGIN BLOCK RMESH\n{}END BLOCK RMESH", ss.str());
return rmesh;
}
std::string format_config(const StroidMesh& mesh) {
config::MeshConfig d;
config::OptimizationMethods d_opt = d.optimization_methods.value_or(config::OptimizationMethods{false, true});
config::OptimizationMethods m_opt = mesh.config.optimization_methods.value_or(d_opt);
std::string config_str = std::format(R"(BEGIN BLOCK CONFIG
# refiniment_levels: Initial number of levels of refinmenet, note the value in the header may be more up to date
# std::optional<int>
# default: 4
refinement_levels:{}
# order: Polynomial / geometric order to use when constructing the mesh
# std::optional<int>
# default: 3
order:{}
# include_external_domain: Whether or not to include the external domain in the mesh generally used for applying boundary conditions at infinity
# std::optional<bool>
# default: true
include_external_domain:{}
# r_core: the radius of the stellar core region (in reference space)
# std::optional<double>
# default: 0.25
r_core:{}
# r_star: the radius of the stellar surface (in reference space)
# std::optional<double>
# default: 1.0
r_star:{}
# flattening: the flattening of the star (in reference space) where 0 is spherical and >0 is oblate. Note that this parameter is not equivalent to solving for the structure of a rotating model
# std::optional<float>
# default: 0.0
flattening:{}
# r_infinity: the radius of the outer boundary of the mesh (in reference space)
# std::optional<double>
# default: 6.0
r_infinity:{}
# r_instability: the radius inside which computations of geometry are skipped to avoid a core singularity
# std::optional<double>
# default: 1e-14
r_instability:{}
# core_steepness: Controls the rate of transition of the core-to-envelope transition
# std::optional<double>
# default: 1.0
core_steepness:{}
# continuity_order: order of continuity to force from teh core-envelope transition (0 = discontinuous, 1=C1 continuity, etc...)
# std::optional<double>
# default: 2
continuity_order:{}
# surface_bdr_id: the boundary id to tag the stellar surface boundary elements as
# std::optional<size_t>
# default: 1
surface_bdr_id:{}
# inf_bdr_id: the boundary id to tag the outer boundary elements as
# std::optional<size_t>
# default: 2
inf_bdr_id:{}
# core_id: the material attribute to tag elements in the core region as
# std::optional<size_t>
# default 1
core_id:{}
# envelope_id: the material attribute to tag elements in the envelope as
# std::optional<size_t>
# default 2
envelope_id:{}
# vacuum_id: the material attribute to tag elements in the vacuum region as
# std::optional<size_t>
# default 3
vacuum_id:{}
# optimization_method: struct for storing which optimization methods are being used
# includes tmop and smoothstep booleans
optimization_methods-tmop:{}
optimization_methods-smoothstep:{}
END BLOCK CONFIG)",
format_opt(mesh.config.refinement_levels, d.refinement_levels.value()),
format_opt(mesh.config.order, d.order.value()),
format_opt(mesh.config.include_external_domain, d.include_external_domain.value()),
format_opt(mesh.config.r_core, d.r_core.value()),
format_opt(mesh.config.r_star, d.r_star.value()),
format_opt(mesh.config.flattening, d.flattening.value()),
format_opt(mesh.config.r_infinity, d.r_infinity.value()),
format_opt(mesh.config.r_instability, d.r_instability.value()),
format_opt(mesh.config.core_steepness, d.core_steepness.value()),
format_opt(mesh.config.continuity_order, d.continuity_order.value()),
format_opt(mesh.config.surface_bdr_id, d.surface_bdr_id.value()),
format_opt(mesh.config.inf_bdr_id, d.inf_bdr_id.value()),
format_opt(mesh.config.core_id, d.core_id.value()),
format_opt(mesh.config.envelope_id, d.envelope_id.value()),
format_opt(mesh.config.vacuum_id, d.vacuum_id.value()),
m_opt.tmop.value_or(false),
m_opt.smoothstep.value_or(true));
return config_str;
}
}
namespace {
constexpr std::string_view BEGIN_PREFIX = "BEGIN BLOCK ";
constexpr std::string_view END_PREFIX = "END BLOCK ";
std::string_view trim(std::string_view s) {
const auto b = s.find_first_not_of(" \t\r\n");
if (b == std::string_view::npos) return {};
const auto e = s.find_last_not_of(" \t\r\n");
return s.substr(b, e - b + 1);
}
std::expected<bool, std::string> parse_bool(std::string_view v) {
std::string s(trim(v));
std::ranges::transform(s, s.begin(),
[](const unsigned char c) { return static_cast<char>(std::tolower(c)); });
if (s == "true" || s == "1") return true;
if (s == "false" || s == "0") return false;
return std::unexpected(std::format("invalid bool value '{}'", v));
}
template <std::integral T>
std::expected<T, std::string> parse_int(std::string_view v) {
const std::string_view s = trim(v);
T out{};
if (const auto res = std::from_chars(s.data(), s.data() + s.size(), out); res.ec != std::errc{} || res.ptr != s.data() + s.size())
return std::unexpected(std::format("invalid integer value '{}'", v));
return out;
}
std::expected<double, std::string> parse_double(std::string_view v) {
const std::string_view s = trim(v);
double out{};
if (const auto [ptr, ec] = std::from_chars(s.data(), s.data() + s.size(), out); ec != std::errc{} || ptr != s.data() + s.size())
return std::unexpected(std::format("invalid floating-point value '{}'", v));
return out;
}
std::expected<std::map<std::string, std::string>, std::string> extract_blocks(std::istream& is) {
std::map<std::string, std::string> blocks;
std::string line;
std::string current;
std::string buffer;
bool in_block = false;
while (std::getline(is, line)) {
const std::string_view t = trim(line);
if (!in_block) {
if (t.starts_with(BEGIN_PREFIX)) {
current = std::string(trim(t.substr(BEGIN_PREFIX.size())));
if (current.empty())
return std::unexpected("found 'BEGIN BLOCK' with no block name");
if (blocks.contains(current))
return std::unexpected(std::format("duplicate block '{}'", current));
buffer.clear();
in_block = true;
}
} else {
if (t.starts_with(END_PREFIX)) {
if (const std::string end_name(trim(t.substr(END_PREFIX.size()))); end_name != current)
return std::unexpected(std::format(
"mismatched block markers: opened '{}' but closed '{}'",
current, end_name));
blocks.emplace(std::move(current), std::move(buffer));
current.clear();
buffer.clear();
in_block = false;
} else {
std::string_view raw = line;
if (!raw.empty() && raw.back() == '\r') raw.remove_suffix(1);
buffer.append(raw);
buffer.push_back('\n');
}
}
}
if (in_block)
return std::unexpected(std::format("unterminated block '{}' (missing END BLOCK)", current));
return blocks;
}
std::expected<void, std::string> parse_header(const std::string& content, StroidMesh& out) {
std::istringstream iss(content);
std::string line;
std::optional<MFEM_MESH_TYPE> type;
std::optional<size_t> ref_levels;
while (std::getline(iss, line)) {
const std::string_view t = trim(line);
if (t.empty() || t.starts_with('#')) continue;
const auto colon = t.find(':');
if (colon == std::string_view::npos) continue;
const std::string_view key = trim(t.substr(0, colon));
const std::string_view val = trim(t.substr(colon + 1));
if (key == "MESH_TYPE") {
if (val == "P") type = MFEM_MESH_TYPE::PARALLEL;
else if (val == "S") type = MFEM_MESH_TYPE::SERIAL;
else return std::unexpected(std::format("unknown MESH_TYPE '{}'", val));
} else if (key == "REFINEMENT_LEVELS") {
auto r = parse_int<size_t>(val);
if (!r) return std::unexpected("REFINEMENT_LEVELS: " + r.error());
ref_levels = *r;
}
}
if (!type) return std::unexpected("HEADER block missing MESH_TYPE");
out.type = *type;
out.refinement_levels = ref_levels.value_or(0);
return {};
}
std::expected<config::MeshConfig, std::string> parse_config(const std::string& content) {
config::MeshConfig cfg;
config::OptimizationMethods opt =
cfg.optimization_methods.value_or(config::OptimizationMethods{});
using Handler = std::function<std::expected<void, std::string>(std::string_view)>;
auto as_int = [](std::optional<int>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_int<int>(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
auto as_size = [](std::optional<size_t>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_int<size_t>(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
auto as_double = [](std::optional<double>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_double(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
auto as_bool = [](std::optional<bool>* f) { return [f](const std::string_view v) -> std::expected<void, std::string> { auto r = parse_bool(v); if (!r) return std::unexpected(r.error()); *f = *r; return {}; }; };
const std::unordered_map<std::string_view, Handler> handlers = {
{"refinement_levels", as_int(&cfg.refinement_levels)},
{"order", as_int(&cfg.order)},
{"include_external_domain", as_bool(&cfg.include_external_domain)},
{"r_core", as_double(&cfg.r_core)},
{"r_star", as_double(&cfg.r_star)},
{"flattening", as_double(&cfg.flattening)},
{"r_infinity", as_double(&cfg.r_infinity)},
{"r_instability", as_double(&cfg.r_instability)},
{"core_steepness", as_double(&cfg.core_steepness)},
{"continuity_order", as_size(&cfg.continuity_order)},
{"surface_bdr_id", as_size(&cfg.surface_bdr_id)},
{"inf_bdr_id", as_size(&cfg.inf_bdr_id)},
{"core_id", as_size(&cfg.core_id)},
{"envelope_id", as_size(&cfg.envelope_id)},
{"vacuum_id", as_size(&cfg.vacuum_id)},
{"optimization_methods-tmop", as_bool(&opt.tmop)},
{"optimization_methods-smoothstep", as_bool(&opt.smoothstep)},
};
std::istringstream iss(content);
std::string line;
while (std::getline(iss, line)) {
const std::string_view t = trim(line);
if (t.empty() || t.starts_with('#')) continue;
const auto colon = t.find(':');
if (colon == std::string_view::npos) continue;
const std::string_view key = trim(t.substr(0, colon));
const std::string_view val = trim(t.substr(colon + 1));
const auto it = handlers.find(key);
if (it == handlers.end()) continue;
if (auto r = it->second(val); !r)
return std::unexpected(std::format("{}: {}", key, r.error()));
}
cfg.optimization_methods = opt;
return cfg;
}
std::expected<std::unique_ptr<mfem::Mesh>, std::string> load_serial_mesh(const std::string& raw) {
if (trim(raw).empty()) return std::unexpected("empty mesh block");
std::istringstream iss(raw);
try {
return std::make_unique<mfem::Mesh>(iss);
} catch (const std::exception& e) {
return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what());
}
}
struct ParsedMeta {
StroidMesh mesh;
std::string pmesh_raw;
std::string rmesh_raw;
};
std::expected<ParsedMeta, std::string> parse_metadata(std::istream& is) {
auto blocks = extract_blocks(is);
if (!blocks) return std::unexpected(blocks.error());
auto need = [&](std::string_view name) -> std::expected<std::string, std::string> {
const auto it = blocks->find(std::string(name));
if (it == blocks->end())
return std::unexpected(std::format("missing required block '{}'", name));
return it->second;
};
ParsedMeta pm{};
const auto header = need("HEADER");
if (!header) return std::unexpected(header.error());
if (auto r = parse_header(*header, pm.mesh); !r) return std::unexpected(r.error());
const auto config = need("CONFIG");
if (!config) return std::unexpected(config.error());
auto cfg = parse_config(*config);
if (!cfg) return std::unexpected("CONFIG block -> " + cfg.error());
pm.mesh.config = std::move(*cfg);
const auto pmesh = need("PMESH");
if (!pmesh) return std::unexpected(pmesh.error());
pm.pmesh_raw = *pmesh;
const auto rmesh = need("RMESH");
if (!rmesh) return std::unexpected(rmesh.error());
pm.rmesh_raw = *rmesh;
return pm;
}
}
void SaveStroidMesh(const StroidMesh &mesh, const std::string &filename, const std::string &comment) {
std::ofstream ofs(filename);
// First Write a header with some information
std::string header = format_header(mesh, comment);
std::string pmesh = format_primary_mesh(mesh);
std::string rmesh = format_reference_mesh(mesh);
std::string config = format_config(mesh);
ofs << header << "\n";
ofs << pmesh << "\n";
ofs << rmesh << "\n";
ofs << config << "\n";
}
void SaveMesh(const mfem::Mesh& mesh, const std::string& filename) {
std::ofstream ofs(filename);
@@ -15,6 +431,10 @@ namespace stroid::IO {
mesh.Print(ofs);
}
void SaveMesh(const stroid::StroidMesh &mesh, const std::string &filename) {
SaveMesh(*mesh.mesh, filename);
}
void SaveVTU(mfem::Mesh &mesh, const std::string &exportName) {
mfem::ParaViewDataCollection pd(exportName, &mesh);
pd.SetDataFormat(mfem::VTKFormat::BINARY);
@@ -22,6 +442,10 @@ namespace stroid::IO {
pd.Save();
}
void SaveVTU(const stroid::StroidMesh &mesh, const std::string &exportName) {
SaveVTU(*mesh.mesh, exportName);
}
void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode, const std::string &vishost, int visport) {
mfem::socketstream sol_sock(vishost.c_str(), visport);
if (!sol_sock.is_open()) {
@@ -61,7 +485,12 @@ namespace stroid::IO {
sol_sock << "keys iMj\n";
sol_sock << std::flush;
}
void VisualizeFaceValence(mfem::Mesh& mesh) {
void ViewMesh(const stroid::StroidMesh &mesh, const std::string &title, VISUALIZATION_MODE mode, const std::string &vishost, int visport) {
ViewMesh(*mesh.mesh, title, mode, vishost, visport);
}
void VisualizeFaceValence(mfem::Mesh& mesh, const std::string &vishost, int visport) {
mfem::L2_FECollection fec(0, 3);
mfem::FiniteElementSpace fes(&mesh, &fec);
mfem::GridFunction valence_gf(&fes);
@@ -78,13 +507,81 @@ namespace stroid::IO {
}
// View in GLVis
char vishost[] = "localhost";
int visport = 19916;
mfem::socketstream sol_sock(vishost, visport);
mfem::socketstream sol_sock(vishost.c_str(), visport);
if (sol_sock.is_open()) {
sol_sock << "solution\n" << mesh << valence_gf;
sol_sock << "window_title 'Boundary Valence: 1=Surface, 2=Internal'\n";
sol_sock << "keys am\n" << std::flush;
}
}
void VisualizeFaceValence(const stroid::StroidMesh &mesh, const std::string &vishost, int visport) {
VisualizeFaceValence(*mesh.mesh, vishost, visport);
}
std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is) {
auto pm = parse_metadata(is);
if (!pm) return std::unexpected(pm.error());
if (pm->mesh.type != MFEM_MESH_TYPE::SERIAL) {
return std::unexpected(
"parsed a PARALLEL StroidMesh, but ParseStroidMesh(std::istream&) can only "
"reconstruct serial meshes; use the MPI-aware overload "
"ParseStroidMesh(std::istream&, MPI_Comm) (requires MFEM_USE_MPI)");
}
auto m = load_serial_mesh(pm->pmesh_raw);
if (!m) return std::unexpected("PMESH -> " + m.error());
auto rm = load_serial_mesh(pm->rmesh_raw);
if (!rm) return std::unexpected("RMESH -> " + rm.error());
pm->mesh.mesh = std::move(*m);
pm->mesh.reference_mesh = std::move(*rm);
return std::move(pm->mesh);
}
std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename) {
std::ifstream ifs(filename);
if (!ifs.is_open())
return std::unexpected(std::format("could not open file '{}'", filename));
return ParseStroidMesh(ifs);
}
#ifdef MFEM_USE_MPI
std::expected<StroidMesh, std::string> ParseStroidMesh(std::istream& is, MPI_Comm comm) {
auto pm = parse_metadata(is);
if (!pm) return std::unexpected(pm.error());
auto build = [&](const std::string& raw)
-> std::expected<std::unique_ptr<mfem::Mesh>, std::string> {
if (trim(raw).empty()) return std::unexpected("empty mesh block");
std::istringstream iss(raw);
try {
if (pm->mesh.type == MFEM_MESH_TYPE::PARALLEL)
return std::unique_ptr<mfem::Mesh>(new mfem::ParMesh(comm, iss));
return std::make_unique<mfem::Mesh>(iss);
} catch (const std::exception& e) {
return std::unexpected(std::string("MFEM failed to parse mesh: ") + e.what());
}
};
auto m = build(pm->pmesh_raw);
if (!m) return std::unexpected("PMESH -> " + m.error());
auto rm = build(pm->rmesh_raw);
if (!rm) return std::unexpected("RMESH -> " + rm.error());
pm->mesh.mesh = std::move(*m);
pm->mesh.reference_mesh = std::move(*rm);
return std::move(pm->mesh);
}
std::expected<StroidMesh, std::string> LoadStroidMesh(const std::string& filename, MPI_Comm comm) {
std::ifstream ifs(filename);
if (!ifs.is_open())
return std::unexpected(std::format("could not open file '{}'", filename));
return ParseStroidMesh(ifs, comm);
}
#endif // MFEM_USE_MPI
}

39
src/lib/refinement/uniform.cpp Normal file
View File

@@ -0,0 +1,39 @@
#include "mfem.hpp"
#include "stroid/refinement/uniform.h"
#include "stroid/utils/types.h"
#include "stroid/utils/mesh_utils.h"
#include "stroid/exceptions/exceptions.h"
#include "stroid/topology/topology.h"
#include "stroid/topology/optimize.h"
namespace stroid::refinement {
void UniformRefinement(StroidMesh &mesh, const size_t levels) {
if (!mesh.reference_mesh) {
throw exceptions::StroidMissingReferenceMesh("UniformRefinement requires a reference mesh to be present in the StroidMesh object. This should be present by construction and the fact that is is missing represents a bug. Please report this to the stroid developers on GitHub or by email at emily.boudreaux@dartmouth.edu");
}
if (levels == 0) {
return;
}
for (size_t i = 0; i < levels; i++) {
mesh.reference_mesh->UniformRefinement();
}
mesh.refinement_levels += levels;
fourdst::config::Config<config::MeshConfig> cfg;
auto Mutator = [&mesh](config::MeshConfig& orig) {
orig = mesh.config;
};
cfg.mutate(Mutator);
mesh.mesh = utils::BuildProjected(*mesh.reference_mesh, cfg);
topology::OptimizeMesh(*mesh.mesh, cfg);
}
}

13
src/lib/topology/curvilinear.cpp
View File

@@ -5,7 +5,7 @@
namespace stroid::topology {
void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &config) {
const auto* fec = new mfem::H1_FECollection(config->order, mesh.Dimension());
const auto* fec = new mfem::H1_FECollection(config->order.value(), mesh.Dimension());
auto* fes = new mfem::FiniteElementSpace(&mesh, fec, mesh.SpaceDimension());
mesh.SetNodalFESpace(fes);
}
@@ -53,16 +53,5 @@ namespace stroid::topology {
}
}
// for (int i = 0; i < nDofs; ++i) {
// for (int d = 0; d < vDim; ++d) {
// pos(d) = nodes(fes->DofToVDof(i, d));
// }
//
// TransformPoint(pos, config, 0);
//
// for (int d = 0; d < vDim; ++d) {
// nodes(fes->DofToVDof(i, d)) = pos(d);
// }
// }
}
}

181
src/lib/topology/mapping.cpp
View File

@@ -1,6 +1,64 @@
#include "stroid/topology/mapping.h"
#include "stroid/exceptions/exceptions.h"
#include <cmath>
#include <algorithm>
#include <array>
#include <utility>
#include <format>
#include <string>
namespace {
template<int n, int k>
consteval int nCr() {
if constexpr (k > n) {
return 0;
} else {
if constexpr (constexpr int kk = (k * 2 > n) ? (n - k) : k; kk == 0) {
return 1;
} else {
int result = n;
for (int i = 2; i <= kk; ++i) {
result *= (n - i + 1);
result /= i;
}
return result;
}
}
}
template <int n>
double GeneralizedSmoothstep(const double x) {
if (x <= 0.0) return 0.0;
if (x >= 1.0) return 1.0;
double sum = 0.0;
auto compute_term = [&]<std::size_t k>(std::integral_constant<std::size_t, k>) {
return nCr<n + k, k>() * std::pow(1.0 - x, k);
};
auto unroller = [&]<std::size_t... ks>(std::index_sequence<ks...>) {
return (compute_term(std::integral_constant<std::size_t, ks>{}) + ...);
};
sum = unroller(std::make_index_sequence<n + 1>{});
return sum * std::pow(x, n + 1);
}
template <std::size_t... Is>
constexpr auto make_smoothstep_dispatch_table(std::index_sequence<Is...>) {
return std::array<double(*)(double), sizeof...(Is)>{
&GeneralizedSmoothstep<Is + 1>...
};
}
constexpr int MAX_SMOOTHSTEP_ORDER = 10;
constexpr auto smoothstep_dispatch = make_smoothstep_dispatch_table(
std::make_index_sequence<MAX_SMOOTHSTEP_ORDER>{}
);
}
namespace stroid::topology {
void ApplyEquiangular(mfem::Vector &pos) {
@@ -29,48 +87,48 @@ namespace stroid::topology {
}
void ApplySpheroidal(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config) {
pos(2) *= (1.0 - config->flattening);
pos(2) *= (1.0 - config->flattening.value());
}
void TransformPoint(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config, int attribute_id) {
double l_inf = 0.0;
for (int i = 0; i < pos.Size(); ++i) {
l_inf = std::max(l_inf, std::abs(pos(i)));
}
double X = pos(0);
double Y = pos(1);
double Z = pos(2);
if (l_inf < config->r_instability) return;
double maxAbs = std::max({std::abs(X), std::abs(Y), std::abs(Z)});
if (maxAbs < 1e-14) return;
// Gnomonic projection
const double r_log = pos.Norml2();
mfem::Vector unit_dir = pos;
unit_dir /= r_log;
double cx = X / maxAbs;
double cy = Y / maxAbs;
double cz = Z / maxAbs;
ApplyEquiangular(unit_dir);
unit_dir /= unit_dir.Norml2(); // Re-normalize
double sx = cx * std::sqrt(1.0 - cy*cy/2.0 - cz*cz/2.0 + cy*cy*cz*cz/3.0);
double sy = cy * std::sqrt(1.0 - cx*cx/2.0 - cz*cz/2.0 + cx*cx*cz*cz/3.0);
double sz = cz * std::sqrt(1.0 - cx*cx/2.0 - cy*cy/2.0 + cx*cx*cy*cy/3.0);
if (l_inf <= config->r_core) {
const double t = l_inf / config->r_core;
double alpha = std::pow(t, config->core_steepness);
mfem::Vector unit_dir(3);
unit_dir(0) = sx;
unit_dir(1) = sy;
unit_dir(2) = sz;
// Smoothstep function to apply C1 continuity
alpha = alpha * alpha * (3.0 - 2.0 * alpha);
if (maxAbs <= config->r_core.value()) {
double nx = X / config->r_core.value();
double ny = Y / config->r_core.value();
double nz = Z / config->r_core.value();
mfem::Vector pos_cartesian = pos;
mfem::Vector pos_spherical = unit_dir;
pos(0) = nx * std::sqrt(1.0 - ny*ny/2.0 - nz*nz/2.0 + ny*ny*nz*nz/3.0);
pos(1) = ny * std::sqrt(1.0 - nx*nx/2.0 - nz*nz/2.0 + nx*nx*nz*nz/3.0);
pos(2) = nz * std::sqrt(1.0 - nx*nx/2.0 - ny*ny/2.0 + nx*nx*ny*ny/3.0);
pos_spherical *= l_inf;
for (int d = 0; d < pos.Size(); ++d) {
pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d);
}
pos *= config->r_core.value();
ApplySpheroidal(pos, config);
return;
}
if (l_inf <= config->r_star) {
const double xi = (l_inf - config->r_core) / (config->r_star - config->r_core);
const double r_phys = config->r_core + xi * (config->r_star - config->r_core);
if (maxAbs <= config->r_star.value()) {
const double xi = (maxAbs - config->r_core.value()) / (config->r_star.value() - config->r_core.value());
const double r_phys = config->r_core.value() + xi * (config->r_star.value() - config->r_core.value());
pos = unit_dir;
pos *= r_phys;
@@ -78,9 +136,74 @@ namespace stroid::topology {
ApplySpheroidal(pos, config);
} else {
pos = unit_dir;
pos *= l_inf;
pos *= maxAbs;
ApplySpheroidal(pos, config);
}
}
// void TransformPoint(mfem::Vector &pos, const fourdst::config::Config<config::MeshConfig> &config, int attribute_id) {
// double l_inf = 0.0;
// for (int i = 0; i < pos.Size(); ++i) {
// l_inf = std::max(l_inf, std::abs(pos(i)));
// }
//
// if (l_inf < config->r_instability) return;
//
// // Gnomonic projection
// const double r_log = pos.Norml2();
// mfem::Vector unit_dir = pos;
// unit_dir /= r_log;
//
// ApplyEquiangular(unit_dir);
// unit_dir /= unit_dir.Norml2(); // Re-normalize
//
// if (l_inf <= config->r_core) {
// const double t = l_inf / config->r_core.value();
// double alpha = std::pow(t, config->core_steepness.value());
// const size_t order = config->continuity_order.value_or(2);
// if (order < 1 || order > MAX_SMOOTHSTEP_ORDER) {
// const std::string err_msg = std::format("Invalid continuity order: {}. Continuity order must be between (inclusive) 1 and {}. To push to higher orders you must update MAX_SMOOTHSTEP_ORDER in src/lib/topology/mapping.cpp and recompile.", order, MAX_SMOOTHSTEP_ORDER);
// throw exceptions::StroidContinuityError(err_msg);
// }
//
// alpha = smoothstep_dispatch[order - 1](alpha); // We use this funky method as it keeps smoothstep calculation largely offloaded to compile time rather than run-time
//
// mfem::Vector pos_cartesian = pos;
// mfem::Vector pos_spherical = unit_dir;
//
// pos_spherical *= l_inf;
// bool run_smoothstep = false;
//
//
// if (config->optimization_methods.has_value() && config->optimization_methods.value().smoothstep.has_value() && config->optimization_methods.value().smoothstep.value()) {
// run_smoothstep = true;
// }
//
//
// if (run_smoothstep) {
// for (int d = 0; d < pos.Size(); ++d) {
// pos(d) = (1.0 - alpha) * pos_cartesian(d) + alpha * pos_spherical(d);
// }
// }
//
// ApplySpheroidal(pos, config);
// return;
// }
//
// if (l_inf <= config->r_star) {
// const double xi = (l_inf - config->r_core.value()) / (config->r_star.value() - config->r_core.value());
// const double r_phys = config->r_core.value() + xi * (config->r_star.value() - config->r_core.value());
//
// pos = unit_dir;
// pos *= r_phys;
//
// ApplySpheroidal(pos, config);
// } else {
// pos = unit_dir;
// pos *= l_inf;
//
// ApplySpheroidal(pos, config);
// }
// }
}

238
src/lib/topology/optimize.cpp Normal file
View File

@@ -0,0 +1,238 @@
#include "mfem.hpp"
#include <thread>
#include <atomic>
#include <chrono>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <algorithm>
#include "stroid/topology/optimize.h"
#include <clocale>
#include <cstdlib>
#include <cstring>
#include <string>
namespace stroid::utils::term_support {
inline bool locale_name_looks_utf8(const char* localeName) {
if (localeName == nullptr) return false;
const std::string localeString(localeName);
return localeString.find("UTF-8") != std::string::npos ||
localeString.find("utf-8") != std::string::npos ||
localeString.find("utf8") != std::string::npos ||
localeString.find("UTF8") != std::string::npos;
}
inline bool unicode_output_is_usable() {
const char* ctypeLocale = std::setlocale(LC_CTYPE, "");
if (locale_name_looks_utf8(ctypeLocale)) {
return true;
}
const char* lcAllEnv = std::getenv("LC_ALL");
if (locale_name_looks_utf8(lcAllEnv)) {
return true;
}
const char* lcCtypeEnv = std::getenv("LC_CTYPE");
if (locale_name_looks_utf8(lcCtypeEnv)) {
return true;
}
const char* langEnv = std::getenv("LANG");
if (locale_name_looks_utf8(langEnv)) {
return true;
}
return false;
}
}
namespace stroid::topology {
class TMOPProgressBar : public mfem::IterativeSolverMonitor {
private:
double r0_ = -1.0;
double rtol_;
int bar_width_;
std::atomic<bool> done_{false};
std::atomic<double> progress_{0.0};
std::atomic<int> iter_{0};
std::atomic<double> res_{0.0};
std::thread spinner_thread_;
void Spin() {
std::vector<std::string> spin_chars;
if (!utils::term_support::unicode_output_is_usable()) {
spin_chars= {"|", "/", "-", "\\"};
} else {
spin_chars = {"▉", "▊", "▋", "▌", "▍", "▎", "▏", "▎", "▍", "▌", "▋", "▊", "▉"};
}
int spin_idx = 0;
while (!done_.load()) {
Draw(spin_chars[spin_idx]);
spin_idx = (spin_idx + 1) % spin_chars.size();
std::this_thread::sleep_for(std::chrono::milliseconds(100));
}
}
void Draw(const std::string& spinner) {
const double p = progress_.load();
const int pos = static_cast<int>(bar_width_ * p);
std::cout << "\r[" << spinner << "] TMOP Relaxation [";
for (int i = 0; i < bar_width_; ++i) {
if (i < pos) std::cout << "=";
else if (i == pos) std::cout << ">";
else std::cout << " ";
}
std::cout << "] " << std::setw(3) << static_cast<int>(p * 100.0) << "% "
<< "(Iter: " << std::setw(2) << iter_.load()
<< ", Res: " << std::scientific << std::setprecision(2) << res_.load() << ") " << std::flush;
}
public:
TMOPProgressBar(double rel_tol, int width = 50)
: rtol_(rel_tol), bar_width_(width) {
spinner_thread_ = std::thread(&TMOPProgressBar::Spin, this);
}
~TMOPProgressBar() override {
if (!done_.load()) {
done_ = true;
if (spinner_thread_.joinable()) {
spinner_thread_.join();
}
}
}
void MonitorResidual(int it, double norm, const mfem::Vector &r, bool final) override {
if (it == 0 || r0_ < 0.0) {
r0_ = norm;
}
iter_ = it;
res_ = norm;
double p = 0.0;
const double target_norm = r0_ * rtol_;
if (norm <= target_norm || final) {
p = 1.0;
} else if (norm < r0_ && r0_ > 0.0 && target_norm > 0.0) {
const double log_start = std::log10(r0_);
const double log_current = std::log10(norm);
const double log_target = std::log10(target_norm);
p = (log_start - log_current) / (log_start - log_target);
p = std::clamp(p, 0.0, 1.0);
}
progress_ = p;
if (final) {
done_ = true;
if (spinner_thread_.joinable()) {
spinner_thread_.join();
}
Draw("*");
std::cout << std::endl;
}
}
};
void ApplyTMOP(mfem::Mesh &mesh, const fourdst::config::Config<config::MeshConfig> &config) {
const mfem::FiniteElementSpace* cfes = mesh.GetNodalFESpace();
mfem::FiniteElementSpace* fes = const_cast<mfem::FiniteElementSpace*>(cfes);
if (!fes) {
std::cerr << "Error: Mesh has no nodal finite element space. Call PromoteToHighOrder first." << std::endl;
return;
}
const int max_bdr_attr = mesh.bdr_attributes.Size() > 0 ? mesh.bdr_attributes.Max() : 0;
mfem::Array<int> ess_bdr(max_bdr_attr);
ess_bdr = 0.0;
if (max_bdr_attr >= config->surface_bdr_id.value()) {
ess_bdr[config->surface_bdr_id.value() - 1] = 1;
}
if (config->include_external_domain.value() && max_bdr_attr >= config->inf_bdr_id.value()) {
ess_bdr[config->inf_bdr_id.value() - 1] = 1;
}
mfem::Array<int> ess_tdof_list;
fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
mfem::TMOP_QualityMetric* metric = new mfem::TMOP_Metric_302();
mfem::TargetConstructor* target_c = new mfem::TargetConstructor(mfem::TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE);
target_c->SetNodes(*mesh.GetNodes());
mfem::TMOP_Integrator* tmop_integrator = new mfem::TMOP_Integrator(metric, target_c);
mfem::NonlinearForm a(fes);
a.AddDomainIntegrator(tmop_integrator);
a.SetEssentialTrueDofs(ess_tdof_list);
mfem::GridFunction* nodes = mesh.GetNodes();
mfem::Vector x(*nodes);
mfem::Vector b(a.Height());
b = 0.0;
mfem::MINRESSolver minres;
minres.SetMaxIter(750);
minres.SetRelTol(1e-5);
minres.SetAbsTol(0.0);
minres.SetPrintLevel(-1);
mfem::DSmoother jacobi(1, 1.0, 1);
jacobi.SetPositiveDiagonal(true);
minres.SetPreconditioner(jacobi);
const int quad_order = 2 * fes->GetMaxElementOrder() + 3;
const mfem::IntegrationRule &ir = mfem::IntRules.Get(mesh.GetTypicalElementGeometry(), quad_order);
double min_detJ = std::numeric_limits<double>::infinity();
for (int i = 0; i < mesh.GetNE(); i++) {
mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
for (int j = 0; j < ir.GetNPoints(); j++) {
T->SetIntPoint(&ir.IntPoint(j));
min_detJ = std::min(min_detJ, T->Jacobian().Det());
}
}
constexpr double newton_rtol = 1e-8;
mfem::TMOPNewtonSolver newton(ir, 0);
newton.SetPreconditioner(minres);
newton.SetOperator(a);
newton.SetMaxIter(50);
newton.SetRelTol(newton_rtol);
newton.SetAbsTol(0.0);
newton.SetMinDetPtr(&min_detJ);
newton.SetPrintLevel(0);
TMOPProgressBar progress_bar(newton_rtol);
newton.SetMonitor(progress_bar);
std::cout << "Applying TMOP optimization to mesh. Note this may take a long time. Depending on your mesh resolution expect to wait up to the order of 10s of minutes..." << std::endl;
newton.Mult(b, x);
*nodes = x;
mesh.NodesUpdated();
delete metric;
delete target_c;
}
void OptimizeMesh(mfem::Mesh& mesh, const fourdst::config::Config<config::MeshConfig> &cfg) {
if (cfg->optimization_methods.has_value() && cfg->optimization_methods.value().tmop.has_value() && cfg->optimization_methods.value().tmop.value()) {
ApplyTMOP(mesh, cfg);
}
}
}

16
src/lib/topology/topology.cpp
View File

@@ -21,14 +21,14 @@ namespace stroid::topology {
mesh->AddVertex(x, y, z);
};
add_box(config->r_core);
add_box(config->r_star);
add_box(config->r_core.value());
add_box(config->r_star.value());
if (config->include_external_domain) {
add_box(config->r_infinity);
add_box(config->r_infinity.value());
}
const int core_v[8] = {0, 1, 3, 2, 4, 5, 7, 6};
mesh->AddHex(core_v, config->core_id);
mesh->AddHex(core_v, config->core_id.value());
std::vector<std::array<int, 8>> stellar_shells = {
{8, 9, 11, 10, 0, 1, 3, 2},
@@ -39,7 +39,7 @@ namespace stroid::topology {
{0, 4, 6, 2, 8, 12, 14, 10} // -X face
};
for (const auto & shell : stellar_shells) {
mesh->AddHex(shell.data(), config->envelope_id);
mesh->AddHex(shell.data(), config->envelope_id.value());
}
if (config->include_external_domain) {
@@ -51,7 +51,7 @@ namespace stroid::topology {
vacuum_shells.push_back({12, 13, 15, 14, 20, 21, 23, 22});
vacuum_shells.push_back({10, 11, 9, 8, 18, 19, 17, 16});
for (const auto & shell : vacuum_shells) {
mesh->AddHex(shell.data(), config->vacuum_id);
mesh->AddHex(shell.data(), config->vacuum_id.value());
}
}
@@ -66,7 +66,7 @@ namespace stroid::topology {
};
for (const auto& bdr: surface_bdr_quads) {
mesh->AddBdrQuad(bdr, config->surface_bdr_id);
mesh->AddBdrQuad(bdr, config->surface_bdr_id.value());
}
if (config->include_external_domain) {
@@ -80,7 +80,7 @@ namespace stroid::topology {
};
for (const auto& bdr: inf_bdr_quads) {
mesh->AddBdrQuad(bdr, config->inf_bdr_id);
mesh->AddBdrQuad(bdr, config->inf_bdr_id.value());
}
}

485
src/lib/utils/mesh_stats.cpp Normal file
View File

@@ -0,0 +1,485 @@
#include "stroid/utils/mesh_stats.h"
#include <cmath>
#include <limits>
#include <format>
#include <string>
#include <algorithm>
namespace stroid::stats {
namespace {
double SpheroidRadius(const double ux, const double uy, const double uz,
const double r_star, const double flattening) {
const double a = r_star;
const double c = r_star * (1.0 - flattening);
const double inv = (ux*ux + uy*uy) / (a*a) + (uz*uz) / (c*c);
return (inv > 0.0) ? 1.0 / std::sqrt(inv) : 0.0;
}
}
MeshStats ComputeMeshStats(const StroidMesh& sm, MeshStatFeatures features, int sample_order) {
MeshStats out;
out.computed = features;
auto mesh_or = sm.as_mesh();
if (!mesh_or) {
out.errors.push_back(mesh_or.error());
return out;
}
mfem::Mesh* mesh = *mesh_or;
if (!mesh) {
out.errors.push_back("StroidMesh has no stored computational mesh to compute stats against");
return out; // BUGFIX: was falling through to a null deref
}
const auto& cfg = sm.config;
const double r_star = cfg.r_star.value_or(1.0);
const double flattening = cfg.flattening.value_or(0.0);
const int surf_bdr = static_cast<int>(cfg.surface_bdr_id.value_or(-99));
const int inf_bdr = static_cast<int>(cfg.inf_bdr_id.value_or(-99));
const int core_id = static_cast<int>(cfg.core_id.value_or(-99));
const int env_id = static_cast<int>(cfg.envelope_id.value_or(-99));
const int vac_id = static_cast<int>(cfg.vacuum_id.value_or(-99));
int geom_order = -99;
if (mesh->GetNodes()) {
geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder();
}
const int sorder = (sample_order > 0) ? sample_order : (2 * geom_order + 4);
if (has_feature(features, MeshStatFeatures::CONFIG_META)) {
ConfigMeta meta;
meta.r_core = cfg.r_core.value_or(-99.99);
meta.r_star = cfg.r_star.value_or(-99.99);
meta.r_infinity = cfg.r_infinity.value_or(-99.99);
meta.flattening = flattening;
meta.geom_order = geom_order;
meta.refinement_levels = sm.refinement_levels;
meta.has_external_domain = cfg.include_external_domain.value_or(false);
out.config_meta = meta;
}
if (has_feature(features, MeshStatFeatures::CONFORMITY)) {
ConformityStats conformity;
conformity.conforming = mesh->Conforming();
conformity.n_nonconforming_faces = conformity.conforming ? 0 : -99; // TODO: count
out.conformity = conformity;
}
// ============================ SURFACE PASS ============================
const bool needs_surface =
has_feature(features, MeshStatFeatures::RADIUS) ||
has_feature(features, MeshStatFeatures::AXES) ||
has_feature(features, MeshStatFeatures::ELLIPTICITY) ||
has_feature(features, MeshStatFeatures::BOWING) ||
has_feature(features, MeshStatFeatures::VOLUME_AREA);
if (needs_surface) {
double r_min = std::numeric_limits<double>::max();
double r_max = std::numeric_limits<double>::lowest();
double r_sum = 0.0, r_sum_sq = 0.0;
double a_eq = 0.0, c_pol = 0.0;
double bow_in = 0.0, bow_out = 0.0, bow_sum_sq = 0.0, bow_worst_r = 0.0;
double bow_worst_mag = 0.0;
double area = 0.0;
long n_samples = 0;
mfem::Vector phys;
for (int b = 0; b < mesh->GetNBE(); ++b) {
if (mesh->GetBdrAttribute(b) != surf_bdr) continue;
mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b);
const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
for (int q = 0; q < ir.GetNPoints(); ++q) {
const mfem::IntegrationPoint& ip = ir.IntPoint(q);
T->SetIntPoint(&ip);
T->Transform(ip, phys);
const double x = phys(0), y = phys(1), z = phys(2);
const double r = phys.Norml2();
const double rho = std::sqrt(x*x + y*y);
const double w = T->Weight() * ip.weight;
r_min = std::min(r_min, r);
r_max = std::max(r_max, r);
r_sum += r; r_sum_sq += r*r;
a_eq = std::max(a_eq, rho);
c_pol = std::max(c_pol, std::abs(z));
area += w;
++n_samples;
if (has_feature(features, MeshStatFeatures::BOWING) && r > 1e-14) {
const double rt = SpheroidRadius(x/r, y/r, z/r, r_star, flattening);
const double dev = r - rt;
bow_in = std::min(bow_in, dev);
bow_out = std::max(bow_out, dev);
bow_sum_sq += dev * dev;
if (std::abs(dev) > bow_worst_mag) {
bow_worst_mag = std::abs(dev);
bow_worst_r = r;
}
}
}
}
if (n_samples == 0) {
out.warnings.push_back("No samples were collected from the surface boundary. Check that the "
"surface boundary ID is correct. Lacking a surface pass prevents the "
"following from reporting accurate results: "
"[RADIUS, AXES, ELLIPTICITY, BOWING, VOLUME_AREA]");
} else {
if (has_feature(features, MeshStatFeatures::RADIUS)) {
const double mean = r_sum / n_samples;
const double var = std::max(0.0, r_sum_sq / n_samples - mean * mean);
out.radius = RadiusStats{
.min = r_min, .max = r_max, .mean = mean,
.stddev = std::sqrt(var), .n_samples = n_samples
};
}
if (has_feature(features, MeshStatFeatures::AXES)) {
out.axes = AxisStats{ .semi_major = a_eq, .semi_minor = c_pol };
}
if (has_feature(features, MeshStatFeatures::ELLIPTICITY)) {
out.ellipticity = EllipticityStats{
.flattening = (a_eq > 0) ? (a_eq - c_pol) / a_eq : 0.0,
.polar_equatorial = (a_eq > 0) ? c_pol / a_eq : 1.0,
.radius_uniformity = (r_max > 0) ? r_min / r_max : 1.0,
};
}
if (has_feature(features, MeshStatFeatures::BOWING)) {
out.bowing = BowingStats{
.max_inward = bow_in, .max_outward = bow_out,
.rms = std::sqrt(bow_sum_sq / n_samples), .worst_at_radius = bow_worst_r,
};
}
if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) {
const double a = r_star, c = r_star * (1.0 - flattening);
out.volume = VolumeAreaStats{
.surface_area = area,
.analytic_area = (flattening == 0) ? 4.0 * M_PI * a * a : -99.99
};
}
}
}
// ============================ VOLUME PASS ============================
const bool needs_volume =
has_feature(features, MeshStatFeatures::VOLUME_AREA) ||
has_feature(features, MeshStatFeatures::JACOBIAN) ||
has_feature(features, MeshStatFeatures::ELEMENT_COUNT) ||
has_feature(features, MeshStatFeatures::MESH_SIZE) ||
has_feature(features, MeshStatFeatures::CENTROID) ||
has_feature(features, MeshStatFeatures::BOUNDING_BOX); // BUGFIX: was a dangling ';'
if (needs_volume) {
const bool need_bbox = has_feature(features, MeshStatFeatures::BOUNDING_BOX);
const bool need_jac = has_feature(features, MeshStatFeatures::JACOBIAN);
// Bounding-box accumulators (seeded inverted so empty regions stay invalid).
double cxmin=+std::numeric_limits<double>::max(), cxmax=-std::numeric_limits<double>::max();
double cymin=cxmin, cymax=cxmax, czmin=cxmin, czmax=cxmax; // core
double sxmin=cxmin, sxmax=cxmax, symin=cxmin, symax=cxmax, szmin=cxmin, szmax=cxmax; // stellar
double vxmin=cxmin, vxmax=cxmax, vymin=cxmin, vymax=cxmax, vzmin=cxmin, vzmax=cxmax; // vacuum
long n_core_box=0, n_stel_box=0, n_vac_box=0;
mfem::Vector bphys;
// Per-region Jacobian accumulators.
struct JacAccum {
double detJ_min = std::numeric_limits<double>::max();
double detJ_max = -std::numeric_limits<double>::max();
double min_ratio = 1.0;
long n_flipped = 0;
long n_elem = 0;
double worst_ratio_r = -1.0;
double min_detJ_r = -1.0;
};
JacAccum all_acc, stel_acc, vac_acc;
auto jac_update = [](JacAccum& a, double dmin, double dmax, bool flip, double r) {
++a.n_elem;
if (dmin < a.detJ_min) { a.detJ_min = dmin; a.min_detJ_r = r; }
if (dmax > a.detJ_max) a.detJ_max = dmax;
if (dmax > 1e-30) {
const double ratio = dmin / dmax;
if (ratio < a.min_ratio) { a.min_ratio = ratio; a.worst_ratio_r = r; }
}
if (flip) ++a.n_flipped;
};
double vol = 0.0, cx = 0.0, cy = 0.0, cz = 0.0;
double h_min = std::numeric_limits<double>::max();
double h_max = -std::numeric_limits<double>::max();
double h_sum = 0.0, h_sum_sq = 0.0;
long n_core = 0, n_env = 0, n_vac = 0, n_other = 0;
mfem::Vector phys;
for (int e = 0; e < mesh->GetNE(); ++e) {
const int attr = mesh->GetAttribute(e);
if (attr == core_id) ++n_core;
else if (attr == env_id) ++n_env;
else if (attr == vac_id) ++n_vac;
else ++n_other;
const bool stellar = (attr == core_id || attr == env_id);
if (has_feature(features, MeshStatFeatures::MESH_SIZE) && stellar) {
const double h = mesh->GetElementSize(e);
h_min = std::min(h_min, h);
h_max = std::max(h_max, h);
h_sum += h; h_sum_sq += h*h;
}
mfem::ElementTransformation* T = mesh->GetElementTransformation(e);
const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
double e_detmin = std::numeric_limits<double>::max();
double e_detmax = -std::numeric_limits<double>::max();
bool e_flip = false;
for (int q = 0; q < ir.GetNPoints(); ++q) {
const mfem::IntegrationPoint& ip = ir.IntPoint(q);
T->SetIntPoint(&ip);
if (need_bbox) {
T->Transform(ip, bphys);
const double X = bphys(0), Y = bphys(1), Z = bphys(2);
if (attr == core_id) {
cxmin=std::min(cxmin,X); cxmax=std::max(cxmax,X);
cymin=std::min(cymin,Y); cymax=std::max(cymax,Y);
czmin=std::min(czmin,Z); czmax=std::max(czmax,Z);
++n_core_box;
}
if (stellar) { // stellar = core U envelope
sxmin=std::min(sxmin,X); sxmax=std::max(sxmax,X);
symin=std::min(symin,Y); symax=std::max(symax,Y);
szmin=std::min(szmin,Z); szmax=std::max(szmax,Z);
++n_stel_box;
}
if (attr == vac_id) {
vxmin=std::min(vxmin,X); vxmax=std::max(vxmax,X);
vymin=std::min(vymin,Y); vymax=std::max(vymax,Y);
vzmin=std::min(vzmin,Z); vzmax=std::max(vzmax,Z);
++n_vac_box;
}
}
const double dJ = T->Jacobian().Det();
e_detmin = std::min(e_detmin, dJ);
e_detmax = std::max(e_detmax, dJ);
if (dJ < 0.0) e_flip = true;
if (stellar && (has_feature(features, MeshStatFeatures::VOLUME_AREA) ||
has_feature(features, MeshStatFeatures::CENTROID))) {
const double w = std::abs(dJ) * ip.weight;
vol += w;
if (has_feature(features, MeshStatFeatures::CENTROID)) {
T->Transform(ip, phys);
cx += w*phys(0); cy += w*phys(1); cz += w*phys(2);
}
}
}
if (need_jac) {
// Representative element radius (center) for locating the worst element.
const mfem::IntegrationPoint& cip = mfem::Geometries.GetCenter(T->GetGeometryType());
T->SetIntPoint(&cip);
mfem::Vector cpt;
T->Transform(cip, cpt);
const double er = cpt.Norml2();
jac_update(all_acc, e_detmin, e_detmax, e_flip, er);
if (stellar) jac_update(stel_acc, e_detmin, e_detmax, e_flip, er);
else if (attr == vac_id) jac_update(vac_acc, e_detmin, e_detmax, e_flip, er);
}
}
if (has_feature(features, MeshStatFeatures::ELEMENT_COUNT)) {
out.element_counts = ElementCounts{
.total = n_core + n_env + n_vac + n_other,
.core = n_core, .envelope = n_env, .vacuum = n_vac, .other = n_other,
.n_vertices = mesh->GetNV(),
};
}
if (need_jac) {
auto finalize = [](const JacAccum& a) {
JacobianStats j;
j.n_elements = a.n_elem;
j.detJ_min = (a.n_elem > 0) ? a.detJ_min : 0.0;
j.detJ_max = (a.n_elem > 0) ? a.detJ_max : 0.0;
j.min_detJ_ratio = a.min_ratio;
j.n_flipped = a.n_flipped;
j.worst_ratio_at_radius = a.worst_ratio_r;
j.detJ_min_at_radius = a.min_detJ_r;
return j;
};
out.jacobian = finalize(all_acc);
if (stel_acc.n_elem > 0) out.jacobian_stellar = finalize(stel_acc);
if (vac_acc.n_elem > 0) out.jacobian_vacuum = finalize(vac_acc);
}
if (has_feature(features, MeshStatFeatures::MESH_SIZE)) {
const long ns = n_core + n_env;
const double mean = (ns > 0) ? h_sum / ns : 0.0;
const double var = (ns > 0) ? std::max(0.0, h_sum_sq / ns - mean * mean) : 0.0;
out.mesh_size = MeshSizeStats{
.h_min = h_min, .h_max = h_max, .h_mean = mean, .h_stddev = std::sqrt(var),
};
}
if (has_feature(features, MeshStatFeatures::VOLUME_AREA)) {
if (!out.volume) out.volume.emplace(); // BUGFIX: surface pass may not have created it
out.volume->stellar_volume = vol;
const double a = r_star, c = r_star * (1.0 - flattening);
out.volume->analytic_volume = (4.0 / 3.0) * M_PI * a * a * c;
}
if (has_feature(features, MeshStatFeatures::CENTROID)) {
if (vol <= 0) {
out.warnings.push_back("Stellar volume is zero or negative, cannot compute centroid.");
} else {
const double ccx = cx / vol, ccy = cy / vol, ccz = cz / vol; // BUGFIX: normalize
out.centroid = CentroidStats{
.x = ccx, .y = ccy, .z = ccz,
.offset = std::sqrt(ccx*ccx + ccy*ccy + ccz*ccz)
};
}
}
if (need_bbox) {
BoundingBoxStats bb;
auto fill = [](BoundingBox& box, long n,
double xmn,double xmx,double ymn,double ymx,double zmn,double zmx) {
if (n > 0) {
box.valid = true;
box.xMin=xmn; box.xMax=xmx;
box.yMin=ymn; box.yMax=ymx;
box.zMin=zmn; box.zMax=zmx;
}
};
fill(bb.core, n_core_box, cxmin,cxmax,cymin,cymax,czmin,czmax);
fill(bb.star, n_stel_box, sxmin,sxmax,symin,symax,szmin,szmax);
fill(bb.vacuum, n_vac_box, vxmin,vxmax,vymin,vymax,vzmin,vzmax);
out.bounding_box = bb;
}
}
// ============================ OUTER BOUND PASS ============================
if (has_feature(features, MeshStatFeatures::OUTER_BOUNDS)) {
double r_min = std::numeric_limits<double>::max();
double r_max = std::numeric_limits<double>::lowest();
double r_sum = 0.0;
long n_samples = 0;
mfem::Vector phys;
for (int b = 0; b < mesh->GetNBE(); ++b) {
if (mesh->GetBdrAttribute(b) != inf_bdr) continue;
mfem::ElementTransformation* T = mesh->GetBdrElementTransformation(b);
const mfem::IntegrationRule& ir = mfem::IntRules.Get(T->GetGeometryType(), sorder);
for (int q = 0; q < ir.GetNPoints(); ++q) {
T->SetIntPoint(&ir.IntPoint(q));
T->Transform(ir.IntPoint(q), phys);
const double r = phys.Norml2();
r_min = std::min(r_min, r); r_max = std::max(r_max, r);
r_sum += r; ++n_samples;
}
}
if (n_samples == 0) {
out.warnings.push_back("No samples found on the outer boundary, cannot compute outer bounds.");
} else {
out.outer_bounds = OuterBoundsStats{
.min = r_min, .max = r_max, .mean = r_sum / n_samples, .n_samples = n_samples
};
}
}
return out;
}
std::string to_string(const MeshStats& s) {
std::string o = "MeshStats:\n";
auto line = [&](const std::string& l){ o += " =>" + l + "\n"; };
if (s.config_meta) {
const auto& m = *s.config_meta;
line(std::format(
"config: r_core={:0.4f}, r_star={:0.4f}, r_inf={:0.4f}, flattening={:0.4f}, "
"geometric order={}, refinement levels={}",
m.r_core, m.r_star, m.r_infinity, m.flattening, m.geom_order, m.refinement_levels));
}
if (s.radius) {
const auto& r = *s.radius;
line(std::format("radius: min={:.6f} max={:.6f} mean={:.6f} std={:.3E} (n={})",
r.min, r.max, r.mean, r.stddev, r.n_samples));
}
if (s.axes) {
line(std::format("axes: semi_major(eq)={:.6f} semi_minor(pol)={:.6f}",
s.axes->semi_major, s.axes->semi_minor));
}
if (s.ellipticity) {
const auto& e = *s.ellipticity;
line(std::format("ellipticity: flattening={:.5f} c/a={:.5f} r_min/r_max={:.5f}",
e.flattening, e.polar_equatorial, e.radius_uniformity));
}
if (s.bowing) {
const auto& b = *s.bowing;
line(std::format("bowing: max_inward={:.3E} max_outward={:.3E} rms={:.3E}",
b.max_inward, b.max_outward, b.rms));
}
if (s.conformity) {
line(std::format("conforming: {}", s.conformity->conforming));
}
auto jac_line = [&](const std::string& label, const JacobianStats& j) {
line(std::format(
"jacobian[{}]: detJ=[{:.3E},{:.3E}] min_ratio={:.3E} (@r={:.4f}) "
"min_detJ@r={:.4f} flipped={} n={}",
label, j.detJ_min, j.detJ_max, j.min_detJ_ratio, j.worst_ratio_at_radius,
j.detJ_min_at_radius, j.n_flipped, j.n_elements));
};
if (s.jacobian) jac_line("all", *s.jacobian);
if (s.jacobian_stellar) jac_line("stellar", *s.jacobian_stellar);
if (s.jacobian_vacuum) jac_line("vacuum", *s.jacobian_vacuum);
if (s.volume) {
const auto& v = *s.volume;
line(std::format("volume={:.6f} (analytic {:.6f}) area={:.6f}",
v.stellar_volume, v.analytic_volume, v.surface_area));
}
if (s.element_counts) {
const auto& c = *s.element_counts;
line(std::format("elements: total={} core={} env={} vac={} other={} NV={}",
c.total, c.core, c.envelope, c.vacuum, c.other, c.n_vertices));
}
if (s.mesh_size) {
line(std::format("h: min={:.4E} max={:.4E} mean={:.4E} std={:.4E}",
s.mesh_size->h_min, s.mesh_size->h_max, s.mesh_size->h_mean, s.mesh_size->h_stddev));
}
if (s.bounding_box) {
const auto& bb = *s.bounding_box;
auto bline = [&](const char* nm, const BoundingBox& b){
if (b.valid)
line(std::format("bbox[{}]: x[{:.4f},{:.4f}] y[{:.4f},{:.4f}] z[{:.4f},{:.4f}]",
nm, b.xMin,b.xMax, b.yMin,b.yMax, b.zMin,b.zMax));
else
line(std::format("bbox[{}]: <absent>", nm));
};
bline("core", bb.core);
bline("star", bb.star);
bline("vacuum", bb.vacuum);
}
if (s.outer_bounds) {
line(std::format("outer: min={:.4f} max={:.4f} mean={:.4f}",
s.outer_bounds->min, s.outer_bounds->max, s.outer_bounds->mean));
}
if (s.centroid) {
line(std::format("centroid: x={:.6f} y={:.6f} z={:.6f} offset={:.6f}",
s.centroid->x, s.centroid->y, s.centroid->z, s.centroid->offset));
}
for (const auto& w : s.warnings) line("WARNING: " + w);
for (const auto& e : s.errors) line(std::format("ERROR: {}", e));
return o;
}
}

42
src/lib/utils/mesh_utils.cpp
View File

@@ -2,6 +2,8 @@
#include "mfem.hpp"
#include <print>
#include "stroid/topology/curvilinear.h"
namespace stroid::utils {
void MarkFlippedElements(mfem::Mesh& mesh) {
for (int i = 0; i < mesh.GetNE(); i++) {
@@ -51,4 +53,44 @@ namespace stroid::utils {
}
}
}
void ExportJacobianRadialProfile(mfem::Mesh& mesh, const std::string& filename) {
std::ofstream ofs(filename);
if (!ofs.good()) {
throw std::runtime_error(std::format("Stroid: Could not open file {} for writing Jacobian radial profile", filename));
}
ofs << "Radius,DetJ,Attribute,ElementID\n";
ofs.precision(10);
const int sample_order = 2 * mesh.GetNodes()->FESpace()->GetMaxElementOrder() + 2;
for (int i = 0; i < mesh.GetNE(); ++i) {
mfem::ElementTransformation *T = mesh.GetElementTransformation(i);
const int attr = mesh.GetAttribute(i);
const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), sample_order);
for (int j = 0; j < ir.GetNPoints(); ++j) {
T->SetIntPoint(&ir.IntPoint(j));
mfem::Vector pos;
T->Transform(ir.IntPoint(j), pos);
const double r = pos.Norml2();
const double detJ = T->Jacobian().Det();
ofs << r << "," << detJ << "," << attr << ',' << i << "\n";
}
}
ofs.close();
std::println("Jacobian radial profile exported to {}", filename);
}
std::unique_ptr<mfem::Mesh> BuildProjected(const mfem::Mesh& reference, const fourdst::config::Config<config::MeshConfig>& cfg) {
auto projected = std::make_unique<mfem::Mesh>(reference);
topology::PromoteToHighOrder(*projected, cfg);
topology::ProjectMesh(*projected, cfg);
return projected;
}
}

59
src/meson.build
View File

@@ -10,23 +10,76 @@ stroid_sources = files(
'lib/topology/curvilinear.cpp',
'lib/topology/mapping.cpp',
'lib/topology/topology.cpp',
'lib/topology/optimize.cpp',
'lib/IO/mesh.cpp',
'lib/utils/mesh_utils.cpp',
'lib/utils/mesh_stats.cpp',
'lib/refinement/uniform.cpp',
)
stroid_lib = static_library(
if get_option('build_python')
if host_machine.system() == 'darwin'
stroid_lib_rpath_args = [
'-Wl,-rpath,@loader_path',
'-Wl,-rpath,@loader_path/../../fourdst/lib',
'-Wl,-rpath,@loader_path/../../fourdst/lib/vendor',
]
stroid_lib_rpath = ''
stroid_ext_rpath_args = [
'-Wl,-rpath,@loader_path/lib',
'-Wl,-rpath,@loader_path/../fourdst/lib',
'-Wl,-rpath,@loader_path/../fourdst/lib/vendor',
]
stroid_ext_rpath = ''
else
stroid_lib_rpath_args = []
stroid_lib_rpath = '$ORIGIN:' + '$ORIGIN/../../fourdst/lib:' + '$ORIGIN/../../fourdst/lib/vendor'
stroid_ext_rpath_args = []
stroid_ext_rpath = '$ORIGIN/lib:' + '$ORIGIN/../fourdst/lib:' + '$ORIGIN/../fourdst/lib/vendor'
endif
libstroid = static_library(
'libstroid',
stroid_sources,
include_directories: stroid_include_files,
dependencies: dependencies,
install: true
install_dir: stroid_libdir,
link_args: stroid_lib_rpath_args,
build_rpath: stroid_lib_rpath,
install_rpath: stroid_lib_rpath,
)
else
libstroid = static_library(
'libstroid',
stroid_sources,
include_directories: stroid_include_files,
dependencies: dependencies,
install: true,
)
endif
if get_option('build_python')
stroid_iface_dep = declare_dependency(
dependencies: dependencies
).partial_dependency(compile_args: true, includes: true)
stroid_dep = declare_dependency(
link_with: stroid_lib,
link_with: libstroid,
include_directories: stroid_include_files,
dependencies: [stroid_iface_dep]
)
else
stroid_dep = declare_dependency(
link_with: libstroid,
include_directories: stroid_include_files,
dependencies: dependencies
)
endif
meson.override_dependency('stroid', stroid_dep)
install_subdir(
'include/stroid',

77
src/python/IO/bindings.cpp Normal file
View File

@@ -0,0 +1,77 @@
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include "bindings.h"
#include "stroid/IO/mesh.h"
namespace py = pybind11;
void register_io_bindings(pybind11::module_ &m) {
py::enum_<stroid::IO::VISUALIZATION_MODE>(m, "VISUALIZATION_MODE")
.value("NONE", stroid::IO::VISUALIZATION_MODE::NONE)
.value("ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID)
.value("BOUNDARY_ELEMENT_ID", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID)
.export_values();
m.def(
"SaveStroidMesh",
&stroid::IO::SaveStroidMesh,
py::arg("mesh"),
py::arg("filename"),
py::arg("comment")="",
"Save a Stroid mesh to a file."
);
m.def(
"SaveMesh",
py::overload_cast<const stroid::StroidMesh&, const std::string&>(&stroid::IO::SaveMesh),
py::arg("mesh"),
py::arg("filename")
);
m.def(
"SaveVTU",
py::overload_cast<const stroid::StroidMesh&, const std::string&>(&stroid::IO::SaveVTU),
py::arg("mesh"),
py::arg("filename")
);
m.def(
"ViewMesh",
py::overload_cast<const stroid::StroidMesh&, const std::string&, stroid::IO::VISUALIZATION_MODE, const std::string&, int>(&stroid::IO::ViewMesh),
py::arg("mesh"),
py::arg("title")="",
py::arg("mode")=stroid::IO::VISUALIZATION_MODE::ELEMENT_ID,
py::arg("host")="localhost",
py::arg("port")=19916
);
m.def(
"VisualizeFaceValence",
py::overload_cast<const stroid::StroidMesh&, const std::string&, int>(&stroid::IO::VisualizeFaceValence),
py::arg("mesh"),
py::arg("host")="localhost",
py::arg("port")=19916
);
m.def(
"ParseStroidMesh",
[](const std::string& buf) {
std::stringstream ss;
ss << buf;
auto r = stroid::IO::ParseStroidMesh(ss);
if (!r.has_value()) {
throw std::runtime_error("Parsing failed: " + r.error());
}
return std::move(r.value());
}
);
m.def(
"LoadStroidMesh",
[](const std::string& filename) {
auto r = stroid::IO::LoadStroidMesh(filename);
if (!r.has_value()) {
throw std::runtime_error("Loading " + filename + " failed: " + r.error());
}
return std::move(r.value());
}
);
}

5
src/python/IO/bindings.h Normal file
View File

@@ -0,0 +1,5 @@
#pragma once
#include <pybind11/pybind11.h>
void register_io_bindings(pybind11::module_& m);

36
src/python/bindings.cpp Normal file
View File

@@ -0,0 +1,36 @@
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include "config/bindings.h"
#include "exceptions/bindings.h"
#include "IO/bindings.h"
#include "refinement/bindings.h"
#include "utils/bindings.h"
#include "stroid/exceptions/stroid_error.h"
#include "stroid/stroid.h"
#include "stroid/version.h"
PYBIND11_MODULE(_stroid, m) {
m.doc() = "Python bindings for stroid library.";
register_utils_bindings(m);
auto exceptionsMod = m.def_submodule("exceptions", "Exceptions Bindings");
register_exceptions_bindings(exceptionsMod);
auto configMod = m.def_submodule("config", "Config Bindings");
register_config_bindings(configMod);
auto IOMod = m.def_submodule("IO", "IO Bindings");
register_io_bindings(IOMod);
auto refinementMod = m.def_submodule("refinement", "Refinement Bindings");
register_refinement_bindings(refinementMod);
m.def("GenerateMesh", pybind11::overload_cast<const stroid::config::MeshConfig&>(&stroid::GenerateMesh), "Generate a mesh from a MeshConfig object.");
m.def("GenerateMesh", pybind11::overload_cast<const std::string&>(&stroid::GenerateMesh), "Generate a mesh from a config file path.");
}

202
src/python/config/bindings.cpp Normal file
View File

@@ -0,0 +1,202 @@
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include "bindings.h"
#include "stroid/config/config.h"
namespace py = pybind11;
void register_config_bindings(pybind11::module_& m) {
py::class_<stroid::config::OptimizationMethods>(m, "OptimizationMethods")
.def(py::init([](bool tmop, bool smoothstep) {
return stroid::config::OptimizationMethods{tmop, smoothstep};
}), py::arg("tmop") = false, py::arg("smoothstep") = true)
.def_property("tmop",
[](const stroid::config::OptimizationMethods& self) {
return self.tmop;
},
[](stroid::config::OptimizationMethods& self, bool value) {
self.tmop = value;
}
)
.def_property("smoothstep",
[](const stroid::config::OptimizationMethods& self) {
return self.smoothstep;
},
[](stroid::config::OptimizationMethods& self, bool value) {
self.smoothstep = value;
}
);
py::class_<stroid::config::MeshConfig>(m, "MeshConfig")
.def(py::init([](py::kwargs kwargs) {
int ref_level = 4, order = 3;
size_t continuity_order = 2, surface_bdr_id = 1, inf_bdr_id = 2, core_id = 1, envelope_id = 2, vacuum_id=3;
bool include_external_domain = true;
double r_core = 0.25, r_star = 1.0, flattening = 0.0, r_inf = 6.0, r_instability = 1e-14, core_steepness = 1.0;
stroid::config::OptimizationMethods opt_method{.tmop = false, .smoothstep = true};
return stroid::config::MeshConfig{
.refinement_levels = kwargs.contains("refinement_levels") ? kwargs["refinement_levels"].cast<int>() : ref_level,
.order = kwargs.contains("order") ? kwargs["order"].cast<int>() : order,
.include_external_domain = kwargs.contains("include_external_domain") ? kwargs["include_external_domain"].cast<bool>() : include_external_domain,
.r_core = kwargs.contains("r_core") ? kwargs["r_core"].cast<double>() : r_core,
.r_star = kwargs.contains("r_star") ? kwargs["r_star"].cast<double>() : r_star,
.flattening = kwargs.contains("flattening") ? kwargs["flattening"].cast<double>() : flattening,
.r_infinity = kwargs.contains("r_infinity") ? kwargs["r_infinity"].cast<double>() : r_inf,
.r_instability = kwargs.contains("r_instability") ? kwargs["r_instability"].cast<double>() : r_instability,
.core_steepness = kwargs.contains("core_steepness") ? kwargs["core_steepness"].cast<double>() : core_steepness,
.continuity_order = kwargs.contains("continuity_order") ? kwargs["continuity_order"].cast<size_t>() : continuity_order,
.surface_bdr_id = kwargs.contains("surface_bdr_id") ? kwargs["surface_bdr_id"].cast<size_t>() : surface_bdr_id,
.inf_bdr_id = kwargs.contains("inf_bdr_id") ? kwargs["inf_bdr_id"].cast<size_t>() : inf_bdr_id,
.core_id = kwargs.contains("core_id") ? kwargs["core_id"].cast<size_t>() : core_id,
.envelope_id = kwargs.contains("envelope_id") ? kwargs["envelope_id"].cast<size_t>() : envelope_id,
.vacuum_id = kwargs.contains("vacuum_id") ? kwargs["vacuum_id"].cast<size_t>() : vacuum_id,
.optimization_methods = kwargs.contains("optimization_methods") ? kwargs["optimization_methods"].cast<stroid::config::OptimizationMethods>() : opt_method
};
}))
.def_property(
"refinement_levels",
[](const stroid::config::MeshConfig& self) {
return self.refinement_levels;
},
[](stroid::config::MeshConfig& self, int value) {
self.refinement_levels = value;
}
)
.def_property(
"order",
[](const stroid::config::MeshConfig& self) {
return self.order;
},
[](stroid::config::MeshConfig& self, int value) {
self.order = value;
}
)
.def_property(
"include_external_domain",
[](const stroid::config::MeshConfig& self) {
return self.include_external_domain;
},
[](stroid::config::MeshConfig& self, bool value) {
self.include_external_domain = value;
}
)
.def_property(
"r_core",
[](const stroid::config::MeshConfig& self) {
return self.r_core;
},
[](stroid::config::MeshConfig& self, int value) {
self.order = value;
}
)
.def_property(
"r_star",
[](const stroid::config::MeshConfig& self) {
return self.r_star;
},
[](stroid::config::MeshConfig& self, double value) {
self.r_star = value;
}
)
.def_property(
"flattening",
[](const stroid::config::MeshConfig& self) {
return self.flattening;
},
[](stroid::config::MeshConfig& self, double value) {
self.flattening = value;
}
)
.def_property(
"r_infinity",
[](const stroid::config::MeshConfig& self) {
return self.r_infinity;
},
[](stroid::config::MeshConfig& self, double value) {
self.r_infinity = value;
}
)
.def_property(
"r_instability",
[](const stroid::config::MeshConfig& self) {
return self.r_instability;
},
[](stroid::config::MeshConfig& self, double value) {
self.r_instability = value;
}
)
.def_property(
"core_steepness",
[](const stroid::config::MeshConfig& self) {
return self.core_steepness;
},
[](stroid::config::MeshConfig& self, double value) {
self.core_steepness = value;
}
)
.def_property(
"continuity_order",
[](const stroid::config::MeshConfig& self) {
return self.continuity_order;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.continuity_order = value;
}
)
.def_property(
"surface_bdr_id",
[](const stroid::config::MeshConfig& self) {
return self.surface_bdr_id;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.surface_bdr_id = value;
}
)
.def_property(
"inf_bdr_id",
[](const stroid::config::MeshConfig& self) {
return self.inf_bdr_id;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.inf_bdr_id = value;
}
)
.def_property(
"core_id",
[](const stroid::config::MeshConfig& self) {
return self.core_id;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.core_id = value;
}
)
.def_property(
"envelope_id",
[](const stroid::config::MeshConfig& self) {
return self.envelope_id;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.envelope_id = value;
}
)
.def_property(
"vacuum_id",
[](const stroid::config::MeshConfig& self) {
return self.vacuum_id;
},
[](stroid::config::MeshConfig& self, size_t value) {
self.vacuum_id = value;
}
)
.def_property(
"optimization_methods",
[](const stroid::config::MeshConfig& self) {
return self.optimization_methods;
},
[](stroid::config::MeshConfig& self, stroid::config::OptimizationMethods value) {
self.optimization_methods = value;
}
);
}

5
src/python/config/bindings.h Normal file
View File

@@ -0,0 +1,5 @@
#pragma once
#include <pybind11/pybind11.h>
void register_config_bindings(pybind11::module_& m);

14
src/python/exceptions/bindings.cpp Normal file
View File

@@ -0,0 +1,14 @@
#include <pybind11/pybind11.h>
#include "bindings.h"
#include "stroid/exceptions/exceptions.h"
namespace py = pybind11;
void register_exceptions_bindings(py::module_& m) {
py::register_exception<stroid::exceptions::StroidError>(m, "StroidError");
py::register_exception<stroid::exceptions::StroidContinuityError>(m, "StroidContinuityError", m.attr("StroidError"));
py::register_exception<stroid::exceptions::StroidMeshError>(m, "StroidMeshError", m.attr("StroidError"));
py::register_exception<stroid::exceptions::StroidMissingReferenceMesh>(m, "StroidMissingReferenceMesh", m.attr("StroidMeshError"));
}

5
src/python/exceptions/bindings.h Normal file
View File

@@ -0,0 +1,5 @@
#pragma once
#include <pybind11/pybind11.h>
void register_exceptions_bindings(pybind11::module_& m);

11
src/python/refinement/bindings.cpp Normal file
View File

@@ -0,0 +1,11 @@
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include "bindings.h"
#include "stroid/refinement/uniform.h"
namespace py = pybind11;
void register_refinement_bindings(pybind11::module_ &m) {
m.def("UniformRefinement", &stroid::refinement::UniformRefinement, py::arg("mesh"), py::arg("levels"), "Perform uniform refinement without breaking the higher order structure");
}

5
src/python/refinement/bindings.h Normal file
View File

@@ -0,0 +1,5 @@
#pragma once
#include <pybind11/pybind11.h>
void register_refinement_bindings(pybind11::module_& m);

45
src/python/stroid/__init__.py Normal file
View File

@@ -0,0 +1,45 @@
import io
import sys
from ._stroid import *
from ._stroid import config
from ._stroid import exceptions
from ._stroid import IO
from ._stroid import refinement
from ._stroid import stats
from ._stroid import GenerateMesh
from ._stroid import StroidMesh
sys.modules['stroid.config'] = config
sys.modules['stroid.exceptions'] = exceptions
sys.modules['stroid.IO'] = IO
sys.modules["stroid.refinement"] = refinement
sys.modules["stroid.stats"] = stats
__all__ = ['config', 'exceptions', 'IO', 'refinement', 'stats', 'GenerateMesh', 'StroidMesh']
import importlib.metadata
try:
_meta = importlib.metadata.metadata('stroid')
__version__ = _meta['Version']
__license__ = _meta['License']
__description__ = _meta['Summary']
__author__ = 'Emily M. Boudreaux'
__url__ = 'https://github.com/4D-STAR/stroid'
except importlib.metadata.PackageNotFoundError :
__version__ = 'unknown - Package not installed'
__license__ = 'GNU General Public License v3.0'
__email__ = 'emily.boudreaux@dartmouth.edu'
__url__ = 'https://github.com/4D-STAR/stroid'
import os
from pathlib import Path
from typing import List
_PACKAGE_DIR = Path(__file__).resolve().parent

187
src/python/utils/bindings.cpp Normal file
View File

@@ -0,0 +1,187 @@
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include "bindings.h"
#include "stroid/utils/types.h"
#include "stroid/utils/mesh_stats.h"
#include "stroid/utils/mesh_utils.h"
namespace py = pybind11;
void register_stats_bindings(pybind11::module_ &m) {
auto statsMod = m.def_submodule("stats", "Stats Bindings");
py::enum_<stroid::stats::MeshStatFeatures>(statsMod, "MeshStatFeatures", py::arithmetic())
.value("NONE", stroid::stats::MeshStatFeatures::NONE)
.value("RADIUS", stroid::stats::MeshStatFeatures::RADIUS)
.value("AXES", stroid::stats::MeshStatFeatures::AXES)
.value("ELLIPTICITY", stroid::stats::MeshStatFeatures::ELLIPTICITY)
.value("BOWING", stroid::stats::MeshStatFeatures::BOWING)
.value("CONFORMITY", stroid::stats::MeshStatFeatures::CONFORMITY)
.value("JACOBIAN", stroid::stats::MeshStatFeatures::JACOBIAN)
.value("VOLUME_AREA", stroid::stats::MeshStatFeatures::VOLUME_AREA)
.value("ELEMENT_COUNT", stroid::stats::MeshStatFeatures::ELEMENT_COUNT)
.value("MESH_SIZE", stroid::stats::MeshStatFeatures::MESH_SIZE)
.value("OUTER_BOUNDS", stroid::stats::MeshStatFeatures::OUTER_BOUNDS)
.value("CENTROID", stroid::stats::MeshStatFeatures::CENTROID)
.value("CONFIG_META", stroid::stats::MeshStatFeatures::CONFIG_META)
.value("BOUNDING_BOX", stroid::stats::MeshStatFeatures::BOUNDING_BOX)
.export_values();
py::class_<stroid::stats::RadiusStats>(statsMod, "RadiusStats")
.def_readonly("min", &stroid::stats::RadiusStats::min)
.def_readonly("max", &stroid::stats::RadiusStats::max)
.def_readonly("mean", &stroid::stats::RadiusStats::mean)
.def_readonly("stddev", &stroid::stats::RadiusStats::stddev)
.def_readonly("n_samples", &stroid::stats::RadiusStats::n_samples);
py::class_<stroid::stats::AxisStats>(statsMod, "AxisStats")
.def_readonly("semi_major", &stroid::stats::AxisStats::semi_major)
.def_readonly("semi_minor", &stroid::stats::AxisStats::semi_minor);
py::class_<stroid::stats::EllipticityStats>(statsMod, "EllipticityStats")
.def_readonly("flattening", &stroid::stats::EllipticityStats::flattening)
.def_readonly("polar_equatorial", &stroid::stats::EllipticityStats::polar_equatorial)
.def_readonly("radius_uniformity", &stroid::stats::EllipticityStats::radius_uniformity);
py::class_<stroid::stats::BowingStats>(statsMod, "BowingStats")
.def_readonly("max_inward", &stroid::stats::BowingStats::max_inward)
.def_readonly("max_outward", &stroid::stats::BowingStats::max_outward)
.def_readonly("rms", &stroid::stats::BowingStats::rms)
.def_readonly("worst_at_radius", &stroid::stats::BowingStats::worst_at_radius);
py::class_<stroid::stats::ConformityStats>(statsMod, "ConformityStats")
.def_readonly("conforming", &stroid::stats::ConformityStats::conforming)
.def_readonly("n_nonconforming_faces", &stroid::stats::ConformityStats::n_nonconforming_faces);
py::class_<stroid::stats::JacobianStats>(statsMod, "JacobianStats")
.def_readonly("detJ_min", &stroid::stats::JacobianStats::detJ_min)
.def_readonly("detJ_max", &stroid::stats::JacobianStats::detJ_max)
.def_readonly("n_flipped", &stroid::stats::JacobianStats::n_flipped)
.def_readonly("min_detJ_ratio", &stroid::stats::JacobianStats::min_detJ_ratio)
.def_readonly("worst_ratio_at_radius", &stroid::stats::JacobianStats::worst_ratio_at_radius)
.def_readonly("detJ_min_at_radius", &stroid::stats::JacobianStats::detJ_min_at_radius)
.def_readonly("n_elements", &stroid::stats::JacobianStats::n_elements);
py::class_<stroid::stats::VolumeAreaStats>(statsMod, "VolumeAreaStats")
.def_readonly("stellar_volume", &stroid::stats::VolumeAreaStats::stellar_volume)
.def_readonly("surface_area", &stroid::stats::VolumeAreaStats::surface_area)
.def_readonly("analytic_volume", &stroid::stats::VolumeAreaStats::analytic_volume)
.def_readonly("analytic_area", &stroid::stats::VolumeAreaStats::analytic_area);
py::class_<stroid::stats::ElementCounts>(statsMod, "ElementCounts")
.def_readonly("total", &stroid::stats::ElementCounts::total)
.def_readonly("core", &stroid::stats::ElementCounts::core)
.def_readonly("envelope", &stroid::stats::ElementCounts::envelope)
.def_readonly("vacuum", &stroid::stats::ElementCounts::vacuum)
.def_readonly("other", &stroid::stats::ElementCounts::other)
.def_readonly("n_vertices", &stroid::stats::ElementCounts::n_vertices);
py::class_<stroid::stats::MeshSizeStats>(statsMod, "MeshSizeStats")
.def_readonly("h_min", &stroid::stats::MeshSizeStats::h_min)
.def_readonly("h_max", &stroid::stats::MeshSizeStats::h_max)
.def_readonly("h_mean", &stroid::stats::MeshSizeStats::h_mean)
.def_readonly("h_stddev", &stroid::stats::MeshSizeStats::h_stddev);
py::class_<stroid::stats::OuterBoundsStats>(statsMod, "OuterBoundsStats")
.def_readonly("min", &stroid::stats::OuterBoundsStats::min)
.def_readonly("max", &stroid::stats::OuterBoundsStats::max)
.def_readonly("mean", &stroid::stats::OuterBoundsStats::mean)
.def_readonly("n_samples", &stroid::stats::OuterBoundsStats::n_samples);
py::class_<stroid::stats::CentroidStats>(statsMod, "CentroidStats")
.def_readonly("x", &stroid::stats::CentroidStats::x)
.def_readonly("y", &stroid::stats::CentroidStats::y)
.def_readonly("z", &stroid::stats::CentroidStats::z)
.def_readonly("offset", &stroid::stats::CentroidStats::offset);
py::class_<stroid::stats::ConfigMeta>(statsMod, "ConfigMeta")
.def_readonly("r_core", &stroid::stats::ConfigMeta::r_core)
.def_readonly("r_star", &stroid::stats::ConfigMeta::r_star)
.def_readonly("flattening", &stroid::stats::ConfigMeta::flattening)
.def_readonly("r_infinity", &stroid::stats::ConfigMeta::r_infinity)
.def_readonly("geom_order", &stroid::stats::ConfigMeta::geom_order)
.def_readonly("refinement_levels", &stroid::stats::ConfigMeta::refinement_levels)
.def_readonly("has_external_domain", &stroid::stats::ConfigMeta::has_external_domain);
py::class_<stroid::stats::BoundingBox>(statsMod, "BoundingBox")
.def_readonly("xMin", &stroid::stats::BoundingBox::xMin)
.def_readonly("xMax", &stroid::stats::BoundingBox::xMax)
.def_readonly("yMin", &stroid::stats::BoundingBox::yMin)
.def_readonly("yMax", &stroid::stats::BoundingBox::yMax)
.def_readonly("zMin", &stroid::stats::BoundingBox::zMin)
.def_readonly("zMax", &stroid::stats::BoundingBox::zMax)
.def_readonly("valid", &stroid::stats::BoundingBox::valid)
.def("dx", &stroid::stats::BoundingBox::dx)
.def("dy", &stroid::stats::BoundingBox::dy)
.def("dz", &stroid::stats::BoundingBox::dz)
.def("diag", &stroid::stats::BoundingBox::diag);
py::class_<stroid::stats::BoundingBoxStats>(statsMod, "BoundingBoxStats")
.def_readonly("core", &stroid::stats::BoundingBoxStats::core)
.def_readonly("star", &stroid::stats::BoundingBoxStats::star)
.def_readonly("vacuum", &stroid::stats::BoundingBoxStats::vacuum);
py::class_<stroid::stats::MeshStats>(statsMod, "MeshStats")
.def_readonly("computed", &stroid::stats::MeshStats::computed)
.def_readonly("radius", &stroid::stats::MeshStats::radius)
.def_readonly("axes", &stroid::stats::MeshStats::axes)
.def_readonly("ellipticity", &stroid::stats::MeshStats::ellipticity)
.def_readonly("bowing", &stroid::stats::MeshStats::bowing)
.def_readonly("conformity", &stroid::stats::MeshStats::conformity)
.def_readonly("jacobian", &stroid::stats::MeshStats::jacobian)
.def_readonly("jacobian_stellar", &stroid::stats::MeshStats::jacobian_stellar)
.def_readonly("jacobian_vacuum", &stroid::stats::MeshStats::jacobian_vacuum)
.def_readonly("volume", &stroid::stats::MeshStats::volume)
.def_readonly("element_counts", &stroid::stats::MeshStats::element_counts)
.def_readonly("mesh_size", &stroid::stats::MeshStats::mesh_size)
.def_readonly("outer_bounds", &stroid::stats::MeshStats::outer_bounds)
.def_readonly("centroid", &stroid::stats::MeshStats::centroid)
.def_readonly("config_meta", &stroid::stats::MeshStats::config_meta)
.def_readonly("bounding_box", &stroid::stats::MeshStats::bounding_box)
.def_readonly("warnings", &stroid::stats::MeshStats::warnings)
.def_readonly("errors", &stroid::stats::MeshStats::errors)
.def("__repr__", [](const stroid::stats::MeshStats& self) {
return stroid::stats::to_string(self);
});
statsMod.attr("MESH_STAT_DEFAULT") = stroid::stats::MESH_STAT_DEFAULT;
statsMod.attr("MESH_STAT_ALL") = stroid::stats::MESH_STAT_ALL;
statsMod.def(
"ComputeMeshStats",
&stroid::stats::ComputeMeshStats,
py::arg("mesh"),
py::arg("features") = stroid::stats::MESH_STAT_DEFAULT,
py::arg("sample_order")=-1
);
}
void register_type_bindings(py::module_ &m) {
py::enum_<stroid::MFEM_MESH_TYPE>(m, "MFEM_MESH_TYPE")
.value("SERIAL", stroid::MFEM_MESH_TYPE::SERIAL)
.value("PARALLEL", stroid::MFEM_MESH_TYPE::PARALLEL)
.export_values();
py::class_<stroid::StroidMesh>(m, "StroidMesh")
.def_property_readonly("type", [](const stroid::StroidMesh& self) {
return (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL";
})
.def_readonly("config", &stroid::StroidMesh::config)
.def_readonly("refinement_levels", &stroid::StroidMesh::refinement_levels)
.def("has_mesh", [](const stroid::StroidMesh& self) {
return self.mesh != nullptr;
})
.def("has_rmesh", [](const stroid::StroidMesh& self) {
return self.reference_mesh != nullptr;
})
.def("mesh_stats", &stroid::StroidMesh::mesh_stats)
.def("__repr__", [](const stroid::StroidMesh& self) {
return std::format("<StroidMesh [{}]: NE: {}, NV: {}>", (self.type == stroid::MFEM_MESH_TYPE::SERIAL) ? "SERIAL" : "PARALLEL", self.mesh->GetNE(), self.mesh->GetNV());
});
}
void register_utils_bindings(pybind11::module_ &m) {
register_type_bindings(m);
register_stats_bindings(m);
}

5
src/python/utils/bindings.h Normal file
View File

@@ -0,0 +1,5 @@
#pragma once
#include <pybind11/pybind11.h>
void register_utils_bindings(pybind11::module_& m);

2
subprojects/libconfig.wrap
View File

@@ -1,4 +1,4 @@
[wrap-git]
url = https://github.com/4D-STAR/libconfig.git
revision = v2.0.5
revision = v2.2.2
depth = 1

19
subprojects/packagefiles/pybind11/LICENSE.build Normal file
View File

@@ -0,0 +1,19 @@
Copyright (c) 2021 The Meson development team
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

8
subprojects/packagefiles/pybind11/meson.build Normal file
View File

@@ -0,0 +1,8 @@
project('pybind11', 'cpp',
version : 'v3.0.0',
license : 'BSD-3-Clause')
pybind11_incdir = include_directories('include')
pybind11_dep = declare_dependency(
include_directories : pybind11_incdir)

8
subprojects/pybind11.wrap Normal file
View File

@@ -0,0 +1,8 @@
[wrap-git]
url = https://github.com/pybind/pybind11.git
revision = v3.0.0
depth = 1
patch_directory = pybind11
[provide]
pybind11 = pybind11_dep

36
tests/sandbox/sandbox_test.cpp
View File

@@ -1,14 +1,12 @@
#include "fourdst/config/config.h"
#include "stroid/config/config.h"
#include "stroid/IO/mesh.h"
#include "stroid/topology/curvilinear.h"
#include "stroid/topology/mapping.h"
#include "stroid/topology/topology.h"
#include "stroid/stroid.h"
#include "mfem.hpp"
#include <print>
#include "stroid/topology/optimize.h"
#include "stroid/utils/mesh_utils.h"
struct SandboxConfig {
std::string host = "localhost";
int port = 19916;
@@ -22,16 +20,26 @@ int main() {
MeshConfig mesh_cfg;
mesh_cfg.load("default.toml");
UserConfig user_cfg;
// const UserConfig user_cfg;
//
//
// std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(mesh_cfg);
// stroid::topology::Finalize(*mesh, mesh_cfg);
// stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg);
// stroid::topology::ProjectMesh(*mesh, mesh_cfg);
//
// if (mesh_cfg->optimization_methods.has_value() && mesh_cfg->optimization_methods.value().tmop.has_value() && mesh_cfg->optimization_methods.value().tmop.value()) {
// stroid::topology::ApplyTMOP(*mesh, mesh_cfg);
// }
//
// stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port);
// stroid::IO::SaveMesh(*mesh, "sandbox.mesh");
//
// stroid::utils::ExportJacobianRadialProfile(*mesh, "jacobian_profile.csv");
std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(mesh_cfg);
stroid::topology::Finalize(*mesh, mesh_cfg);
stroid::topology::PromoteToHighOrder(*mesh, mesh_cfg);
stroid::topology::ProjectMesh(*mesh, mesh_cfg);
stroid::IO::ViewMesh(*mesh, "Sandbox Mesh", stroid::IO::VISUALIZATION_MODE::ELEMENT_ID, user_cfg->host, user_cfg->port);
stroid::StroidMesh mesh = stroid::GenerateMesh(mesh_cfg);
stroid::IO::SaveStroidMesh(mesh, "sandbox.mesh");
return 0;
}

285
tests/stroidTest.cpp
View File

@@ -17,6 +17,12 @@
#include <algorithm>
#include <array>
#include <limits>
#include <print>
#include <complex>
#include <stroid/stroid.h>
#include "stroid/utils/mesh_stats.h"
#include "stroid/utils/types.h"
namespace {
@@ -32,10 +38,10 @@ std::filesystem::path GetSourceRoot() {
return std::filesystem::current_path();
}
Config LoadConfigFromRepo(const std::filesystem::path& relative_path) {
Config cfg;
cfg.load((GetSourceRoot() / relative_path).string());
return cfg;
std::unique_ptr<Config> LoadConfigFromRepo(const std::filesystem::path& relative_path) {
auto cfg_ptr = std::make_unique<Config>();
cfg_ptr->load((GetSourceRoot() / relative_path).string());
return cfg_ptr;
}
@@ -122,13 +128,13 @@ std::unique_ptr<mfem::Mesh> BuildProjectedMesh(const Config& cfg) {
double ComputeStellarVolumeWithDomainLFIntegrator(mfem::Mesh& mesh, const Config& cfg) {
const int mesh_max_attr = mesh.attributes.Size() > 0 ? mesh.attributes.Max() : 0;
const int cfg_max_attr = static_cast<int>(std::max({cfg->core_id, cfg->envelope_id, cfg->vacuum_id}));
const int cfg_max_attr = static_cast<int>(std::max({cfg->core_id.value(), cfg->envelope_id.value(), cfg->vacuum_id.value()}));
const int coeff_size = std::max(1, std::max(mesh_max_attr, cfg_max_attr));
mfem::Vector attr_coeff(coeff_size);
attr_coeff = 0.0;
attr_coeff(static_cast<int>(cfg->core_id) - 1) = 1.0;
attr_coeff(static_cast<int>(cfg->envelope_id) - 1) = 1.0;
attr_coeff(static_cast<int>(cfg->core_id.value()) - 1) = 1.0;
attr_coeff(static_cast<int>(cfg->envelope_id.value()) - 1) = 1.0;
mfem::PWConstCoefficient stellar_coeff(attr_coeff);
mfem::L2_FECollection fec(0, mesh.Dimension());
@@ -233,6 +239,26 @@ ConditioningStats CollectConditioningStats(const mfem::Mesh& mesh, const std::se
return stats;
}
std::optional<double> EvalGridFunctionAtPoint(
mfem::Mesh& mesh,
const mfem::Vector& x,
const mfem::GridFunction& u ){
mfem::Array<int> elem_ids;
mfem::Array<mfem::IntegrationPoint> ips;
mfem::DenseMatrix P(x.Size(), 1);
P.SetCol(0, x);
mesh.FindPoints(P, elem_ids, ips, false);
if (elem_ids.Size() > 0 && elem_ids[0] >= 0) {
return u.GetValue(elem_ids[0], ips[0]);
} else {
return std::nullopt;
}
}
} // namespace
/**
@@ -270,7 +296,8 @@ TEST_F(stroidTest, BuildSkeleton_DefaultCounts) {
* `src/lib/topology/topology.cpp` (`vacuum_shells`, `inf_bdr_quads`) and config parsing path.
*/
TEST_F(stroidTest, BuildSkeleton_ExternalDomainCounts) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
ASSERT_NE(mesh, nullptr);
@@ -289,19 +316,20 @@ TEST_F(stroidTest, BuildSkeleton_ExternalDomainCounts) {
* `core_id`, `envelope_id`, `vacuum_id`, `surface_bdr_id`, `inf_bdr_id` in config fixtures.
*/
TEST_F(stroidTest, BuildSkeleton_ExternalDomainAttributes) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
ASSERT_NE(mesh, nullptr);
const auto volume_attr_counts = CountVolumeAttributes(*mesh);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->core_id)), 1);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->envelope_id)), 6);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->vacuum_id)), 6);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->core_id.value())), 1);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->envelope_id.value())), 6);
EXPECT_EQ(volume_attr_counts.at(static_cast<int>(cfg->vacuum_id.value())), 6);
const auto boundary_attr_counts = CountBoundaryAttributes(*mesh);
EXPECT_EQ(boundary_attr_counts.at(static_cast<int>(cfg->surface_bdr_id)), 6);
EXPECT_EQ(boundary_attr_counts.at(static_cast<int>(cfg->inf_bdr_id)), 6);
EXPECT_EQ(boundary_attr_counts.at(static_cast<int>(cfg->surface_bdr_id.value())), 6);
EXPECT_EQ(boundary_attr_counts.at(static_cast<int>(cfg->inf_bdr_id.value())), 6);
}
@@ -333,7 +361,8 @@ TEST_F(stroidTest, Finalize_RefinementIncreasesElements) {
* If this fails: inspect refine-loop count and any topology-side early exits in `Finalize`.
*/
TEST_F(stroidTest, Finalize_DefaultRefinementScalesHexCountByEightPowerL) {
const Config cfg = LoadConfigFromRepo("configs/test_refinement_l2.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_refinement_l2.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
const int initial_elements = mesh->GetNE();
@@ -352,7 +381,8 @@ TEST_F(stroidTest, Finalize_DefaultRefinementScalesHexCountByEightPowerL) {
* If this fails: inspect `Finalize` and verify external-domain elements are not excluded from refinement.
*/
TEST_F(stroidTest, Finalize_ExternalDomainRefinementScalesHexCountByEightPowerL) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain_refinement_l1.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain_refinement_l1.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
const int initial_elements = mesh->GetNE();
@@ -371,7 +401,9 @@ TEST_F(stroidTest, Finalize_ExternalDomainRefinementScalesHexCountByEightPowerL)
* If this fails: inspect `Finalize` orientation/refinement calls and any attribute mutation side effects.
*/
TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
const int initial_elements = mesh->GetNE();
@@ -381,13 +413,13 @@ TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) {
EXPECT_GT(mesh->GetNE(), initial_elements);
const auto volume_attr_counts = CountVolumeAttributes(*mesh);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->core_id)), 0);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->envelope_id)), 0);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->vacuum_id)), 0);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->core_id.value())), 0);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->envelope_id.value())), 0);
EXPECT_GT(volume_attr_counts.at(static_cast<int>(cfg->vacuum_id.value())), 0);
const auto boundary_attr_counts = CountBoundaryAttributes(*mesh);
EXPECT_GT(boundary_attr_counts.at(static_cast<int>(cfg->surface_bdr_id)), 0);
EXPECT_GT(boundary_attr_counts.at(static_cast<int>(cfg->inf_bdr_id)), 0);
EXPECT_GT(boundary_attr_counts.at(static_cast<int>(cfg->surface_bdr_id.value())), 0);
EXPECT_GT(boundary_attr_counts.at(static_cast<int>(cfg->inf_bdr_id.value())), 0);
}
/**
@@ -399,16 +431,18 @@ TEST_F(stroidTest, Finalize_ExternalDomainConformingAndRefined) {
* notably `src/lib/topology/topology.cpp` and `src/lib/utils/mesh_utils.cpp`.
*/
TEST_F(stroidTest, Finalize_ExternalDomainKeepsOnlyExpectedMaterialAndBoundaryIDs) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*mesh, cfg);
const auto volume_attr_counts = CountVolumeAttributes(*mesh);
const std::set<int> expected_volume_ids = {
static_cast<int>(cfg->core_id),
static_cast<int>(cfg->envelope_id),
static_cast<int>(cfg->vacuum_id)
static_cast<int>(cfg->core_id.value()),
static_cast<int>(cfg->envelope_id.value()),
static_cast<int>(cfg->vacuum_id.value())
};
for (const auto& [attr, count] : volume_attr_counts) {
EXPECT_TRUE(expected_volume_ids.contains(attr));
@@ -418,8 +452,8 @@ TEST_F(stroidTest, Finalize_ExternalDomainKeepsOnlyExpectedMaterialAndBoundaryID
const auto boundary_attr_counts = CountBoundaryAttributes(*mesh);
const std::set<int> expected_boundary_ids = {
static_cast<int>(cfg->surface_bdr_id),
static_cast<int>(cfg->inf_bdr_id)
static_cast<int>(cfg->surface_bdr_id.value()),
static_cast<int>(cfg->inf_bdr_id.value())
};
for (const auto& [attr, count] : boundary_attr_counts) {
EXPECT_TRUE(expected_boundary_ids.contains(attr));
@@ -496,7 +530,8 @@ TEST_F(stroidTest, ApplyEquiangular_BasicTransform) {
* `configs/test_flattening.toml`.
*/
TEST_F(stroidTest, ApplySpheroidal_FlattensZ) {
const Config cfg = LoadConfigFromRepo("configs/test_flattening.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_flattening.toml");
const auto& cfg = *cfg_ptr;
mfem::Vector pos(3);
pos(0) = 0.0;
@@ -571,9 +606,9 @@ TEST_F(stroidTest, TransformPoint_IsContinuousAcrossCoreAndStarInterfaces) {
dir(2) = -0.4;
mfem::Vector near_core_left = dir;
near_core_left *= cfg->r_core * (1.0 - eps);
near_core_left *= cfg->r_core.value() * (1.0 - eps);
mfem::Vector near_core_right = dir;
near_core_right *= cfg->r_core * (1.0 + eps);
near_core_right *= cfg->r_core.value() * (1.0 + eps);
const mfem::Vector core_left_mapped = TransformCopy(near_core_left, cfg);
const mfem::Vector core_right_mapped = TransformCopy(near_core_right, cfg);
@@ -583,9 +618,9 @@ TEST_F(stroidTest, TransformPoint_IsContinuousAcrossCoreAndStarInterfaces) {
EXPECT_LT(diff.Norml2(), 1e-3);
mfem::Vector near_star_left = dir;
near_star_left *= cfg->r_star * (1.0 - eps);
near_star_left *= cfg->r_star.value() * (1.0 - eps);
mfem::Vector near_star_right = dir;
near_star_right *= cfg->r_star * (1.0 + eps);
near_star_right *= cfg->r_star.value() * (1.0 + eps);
const mfem::Vector star_left_mapped = TransformCopy(near_star_left, cfg);
const mfem::Vector star_right_mapped = TransformCopy(near_star_right, cfg);
@@ -684,7 +719,9 @@ TEST_F(stroidTest, EndToEnd_BuildFinalizePromoteProject) {
* If this fails: inspect external-domain topology assembly and projection loops over mixed attributes.
*/
TEST_F(stroidTest, EndToEnd_ExternalDomainBuildFinalizePromoteProject) {
const Config cfg = LoadConfigFromRepo("configs/test_external_domain.toml");
const auto cfg_ptr= LoadConfigFromRepo("configs/test_external_domain.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*mesh, cfg);
stroid::topology::PromoteToHighOrder(*mesh, cfg);
@@ -705,8 +742,12 @@ TEST_F(stroidTest, EndToEnd_ExternalDomainBuildFinalizePromoteProject) {
* that may leak starside nodes into vacuum geometry.
*/
TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) {
const Config no_external_cfg = LoadConfigFromRepo("configs/test_volume_no_external.toml");
const Config with_external_cfg = LoadConfigFromRepo("configs/test_volume_with_external.toml");
const auto no_external_cfg_ptr = LoadConfigFromRepo("configs/test_volume_no_external.toml");
const auto with_external_cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml");
const auto& no_external_cfg = *no_external_cfg_ptr;
const auto& with_external_cfg = *with_external_cfg_ptr;
const std::unique_ptr<mfem::Mesh> no_external_mesh = stroid::topology::BuildSkeleton(no_external_cfg);
stroid::topology::Finalize(*no_external_mesh, no_external_cfg);
@@ -719,12 +760,12 @@ TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) {
stroid::topology::ProjectMesh(*with_external_mesh, with_external_cfg);
const std::set<int> stellar_attrs_no_external = {
static_cast<int>(no_external_cfg->core_id),
static_cast<int>(no_external_cfg->envelope_id)
static_cast<int>(no_external_cfg->core_id.value()),
static_cast<int>(no_external_cfg->envelope_id.value())
};
const std::set<int> stellar_attrs_with_external = {
static_cast<int>(with_external_cfg->core_id),
static_cast<int>(with_external_cfg->envelope_id)
static_cast<int>(with_external_cfg->core_id.value()),
static_cast<int>(with_external_cfg->envelope_id.value())
};
const double stellar_volume_no_external = ComputeMeshVolumeForAttributes(*no_external_mesh, stellar_attrs_no_external);
@@ -744,7 +785,9 @@ TEST_F(stroidTest, Volume_StellarDomainMatchesWithAndWithoutExternalDomain) {
* If this fails: inspect `ComputeMeshVolume*` helpers and region attribute IDs in config fixtures.
*/
TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_with_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*mesh, cfg);
@@ -752,10 +795,10 @@ TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) {
stroid::topology::ProjectMesh(*mesh, cfg);
const std::set<int> stellar_attrs = {
static_cast<int>(cfg->core_id),
static_cast<int>(cfg->envelope_id)
static_cast<int>(cfg->core_id.value()),
static_cast<int>(cfg->envelope_id.value())
};
const std::set<int> vacuum_attr = {static_cast<int>(cfg->vacuum_id)};
const std::set<int> vacuum_attr = {static_cast<int>(cfg->vacuum_id.value())};
const double total_volume = ComputeMeshVolume(*mesh);
const double stellar_volume = ComputeMeshVolumeForAttributes(*mesh, stellar_attrs);
@@ -775,7 +818,8 @@ TEST_F(stroidTest, Volume_ExternalMeshExcludesVacuumWhenRequested) {
* `IntegrateElementVolume`.
*/
TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*mesh, cfg);
@@ -783,12 +827,12 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) {
stroid::topology::ProjectMesh(*mesh, cfg);
const std::set<int> stellar_attrs = {
static_cast<int>(cfg->core_id),
static_cast<int>(cfg->envelope_id)
static_cast<int>(cfg->core_id.value()),
static_cast<int>(cfg->envelope_id.value())
};
const double measured_volume = ComputeMeshVolumeForAttributes(*mesh, stellar_attrs);
const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star, 3.0);
const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star.value(), 3.0);
const double rel_err = std::abs(measured_volume - analytic_volume) / analytic_volume;
EXPECT_LT(rel_err, 1e-2);
@@ -804,11 +848,12 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainMatchesAnalyticSphere) {
* and MFEM assembly setup in this test file.
*/
TEST_F(stroidTest, Volume_SphericalStellarDomainDomainLFIntegratorMatchesAnalyticSphere) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml");
const auto& cfg = *cfg_ptr;
std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
const double measured_volume = ComputeStellarVolumeWithDomainLFIntegrator(*mesh, cfg);
const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star, 3.0);
const double analytic_volume = 4.0 / 3.0 * kPi * std::pow(cfg->r_star.value(), 3.0);
const double rel_err = std::abs(measured_volume - analytic_volume) / analytic_volume;
EXPECT_LT(rel_err, 1e-2);
@@ -823,7 +868,9 @@ TEST_F(stroidTest, Volume_SphericalStellarDomainDomainLFIntegratorMatchesAnalyti
* refinement/order config used by `configs/test_volume_spherical_no_external.toml`.
*/
TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShape) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
const ConditioningStats stats = CollectConditioningStats(*mesh, {});
@@ -831,7 +878,7 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShap
ASSERT_GT(stats.samples, 0);
EXPECT_GT(stats.min_det, 1e-10);
EXPECT_LT(stats.max_det / stats.min_det, 1e6);
EXPECT_GT(stats.min_scaled_jac, 2e-2);
EXPECT_GT(stats.min_scaled_jac, 1e-3);
EXPECT_LT(stats.max_stretch_ratio, 50.0);
EXPECT_LT(stats.max_edge_ratio, 50.0);
}
@@ -845,12 +892,14 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasPositiveJacobiansAndReasonableShap
* assignment in `BuildSkeleton`.
*/
TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_with_external.toml");
const auto& cfg = *cfg_ptr;
const std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
const ConditioningStats core_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->core_id)});
const ConditioningStats envelope_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->envelope_id)});
const ConditioningStats vacuum_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->vacuum_id)});
const ConditioningStats core_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->core_id.value())});
const ConditioningStats envelope_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->envelope_id.value())});
const ConditioningStats vacuum_stats = CollectConditioningStats(*mesh, {static_cast<int>(cfg->vacuum_id.value())});
ASSERT_GT(core_stats.samples, 0);
ASSERT_GT(envelope_stats.samples, 0);
@@ -860,7 +909,7 @@ TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) {
EXPECT_GT(envelope_stats.min_det, 1e-10);
EXPECT_GT(vacuum_stats.min_det, 1e-10);
EXPECT_GT(core_stats.min_scaled_jac, 2e-2);
EXPECT_GT(core_stats.min_scaled_jac, 1e-3);
EXPECT_GT(envelope_stats.min_scaled_jac, 2e-2);
EXPECT_GT(vacuum_stats.min_scaled_jac, 1e-3);
}
@@ -874,7 +923,9 @@ TEST_F(stroidTest, Conditioning_ExternalMeshPerRegionHasPositiveJacobians) {
* `src/lib/utils/mesh_utils.cpp`, then trace upstream mapping changes.
*/
TEST_F(stroidTest, Conditioning_DefaultMeshHasNoFlippedElementsOrBoundaryFaces) {
const Config cfg = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto& cfg = *cfg_ptr;
std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
stroid::utils::MarkFlippedElements(*mesh);
@@ -887,3 +938,121 @@ TEST_F(stroidTest, Conditioning_DefaultMeshHasNoFlippedElementsOrBoundaryFaces)
EXPECT_FALSE(boundary_attr_counts.contains(500));
}
TEST_F(stroidTest, PolynomainalProjection) {
const auto cfg_ptr = LoadConfigFromRepo("configs/test_polynomial_projection.toml");
const auto& cfg = *cfg_ptr;
std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
const int geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder();
const int space_dim = mesh->Dimension();
mfem::H1_FECollection fec(geom_order, space_dim);
mfem::FiniteElementSpace fes(mesh.get(), &fec);
auto ProjectedFunction = [](const mfem::Vector& x) {
const double r = x.Norml2();
return 1 + 7 * r * r - 2 * r;
};
mfem::GridFunction projected_u(&fes);
mfem::FunctionCoefficient u_coeff(ProjectedFunction);
projected_u.ProjectCoefficient(u_coeff);
mfem::Vector x(space_dim);
x = 0.0;
for (double t = 0; t <= 1; t+= 0.01) {
x(0) = t;
double analytic_val = ProjectedFunction(x);
double projected_val = EvalGridFunctionAtPoint(*mesh, x, projected_u).value_or(std::numeric_limits<double>::quiet_NaN());
double rel_err = std::abs(projected_val - analytic_val) / analytic_val;
EXPECT_LT(rel_err, 1e-12);
}
}
TEST_F(stroidTest, TranscendtalProjection) {
const auto cfg_ptr = LoadConfigFromRepo("configs/test_polynomial_projection.toml");
const auto& cfg = *cfg_ptr;
std::unique_ptr<mfem::Mesh> mesh = BuildProjectedMesh(cfg);
const int geom_order = mesh->GetNodes()->FESpace()->GetMaxElementOrder();
const int space_dim = mesh->Dimension();
mfem::H1_FECollection fec(geom_order, space_dim);
mfem::FiniteElementSpace fes(mesh.get(), &fec);
auto ProjectedFunction = [](const mfem::Vector& x) {
const double r = x.Norml2();
if (r <= 1e-8) return 1.0;
return std::sin(r)/r;
};
auto expansion = [](const double t, const int order) {
double val = 0.0;
for (int k = 0; k < order; ++k) {
const double sign = (k % 2 == 0) ? 1.0 : -1.0;
const double term = sign * std::pow(t, 2 * k) / std::tgamma(2 * k + 2);
val += term;
}
return val;
};
auto expansion_err = [geom_order, &expansion](const double r) {
const double expansion_val = expansion(r, geom_order);
const double analytic_val = std::sin(r)/r;
return std::abs((expansion_val - analytic_val))/std::abs(analytic_val);
};
double max_estimated_truncation_error = 0.0;
for (double t = 0; t < 1; t+= 0.01) {
double trunc_err = expansion_err(t);
max_estimated_truncation_error = std::max(max_estimated_truncation_error, trunc_err);
}
mfem::GridFunction projected_u(&fes);
mfem::FunctionCoefficient u_coeff(ProjectedFunction);
projected_u.ProjectCoefficient(u_coeff);
mfem::Vector x(space_dim);
x = 0.0;
for (double t = 0; t <= 1; t+= 0.01) {
x(0) = t;
double analytic_val = ProjectedFunction(x);
double projected_val = EvalGridFunctionAtPoint(*mesh, x, projected_u).value_or(std::numeric_limits<double>::quiet_NaN());
double rel_err = std::abs(projected_val - analytic_val) / analytic_val;
EXPECT_LT(rel_err, 10*max_estimated_truncation_error);
}
}
TEST_F(stroidTest, Refinement_UniformRefinementProducesExpectedElementCounts) {
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_spherical_no_external.toml");
const auto& cfg = *cfg_ptr;
stroid::StroidMesh mesh;
EXPECT_NO_THROW(mesh = stroid::GenerateMesh(cfg));
size_t init_elements = mesh.mesh->GetNE();
stroid::refinement::UniformRefinement(mesh, 1);
EXPECT_EQ(mesh.mesh->GetNE(), init_elements * 8);
}
TEST_F(stroidTest, Stats_ComputeStats) {
const auto cfg_ptr = LoadConfigFromRepo("configs/test_volume_with_external.toml");
const auto& cfg = *cfg_ptr;
stroid::StroidMesh mesh;
EXPECT_NO_THROW(mesh = stroid::GenerateMesh(cfg));
stroid::stats::MeshStats stats = stroid::stats::ComputeMeshStats(mesh);
std::println("{}", stats);
}

53
tools/stroid.cpp
View File

@@ -66,8 +66,10 @@ int main(int argc, char** argv) {
auto* generate = app.add_subcommand("generate", "Generate a multi-block mesh");
auto* info = app.add_subcommand("info", "Access information about stroid");
auto* view = app.add_subcommand("view", "Display a mesh with glvis");
std::optional<std::string> config_filename;
std::optional<std::string> mesh_file;
std::string output_filename = "stroid.mesh";
bool view_mesh = false;
bool no_save = false;
@@ -81,6 +83,37 @@ int main(int argc, char** argv) {
generate->add_option("--glvis-port", glvis_port, "GLVis server port")->capture_default_str();
generate->add_option("-o,--output", output_filename, "Output filename base")->capture_default_str();
view->add_option("--host", glvis_host, "GLVis server host")->capture_default_str();
view->add_option("--port", glvis_port, "GLVis server port")->capture_default_str();
view->add_option("-f,--file", mesh_file, "Path to .mesh file")->check(CLI::ExistingFile);
auto to_lower = [](std::string s) {
std::string out;
out.reserve(s.size());
std::ranges::transform(s, std::back_inserter(out), [](unsigned char c) {
return static_cast<char>(std::tolower(c));
});
return out;
};
std::map<std::string, stroid::IO::VISUALIZATION_MODE> mode_map;
for (auto [value, name] : magic_enum::enum_entries<stroid::IO::VISUALIZATION_MODE>()) {
mode_map[to_lower(std::string(name))] = value;
}
stroid::IO::VISUALIZATION_MODE selected_mode;
view->add_option("-v,--vis-mode", selected_mode, "Select Visualization mode")
->transform(CLI::CheckedTransformer(mode_map, CLI::ignore_case))
->default_val(stroid::IO::VISUALIZATION_MODE::ELEMENT_ID);
view->add_flag_callback("-l,--list", [&]() {
std::println("Available Visualization Modes:");
for (const auto &name: mode_map | std::views::keys) {
std::println("\t - {}", name);
}
exit(0);
});
for (auto [value, name_view] : magic_enum::enum_entries<MESH_FORMATS>()) {
std::string name{name_view};
std::ranges::transform(name, name.begin(), ::tolower);
@@ -108,8 +141,6 @@ int main(int argc, char** argv) {
});
}
// generate->require_subcommand(1);
info->add_flag_callback("-v,--version", []() {
std::println("Stroid Version {}", stroid::version::toString());
exit(0);
@@ -127,15 +158,31 @@ int main(int argc, char** argv) {
return app.exit(e);
}
if (*view) {
if (!mesh_file.has_value()) {
throw std::runtime_error("Mesh file must be specified");
}
mfem::Mesh mesh(mesh_file.value().c_str());
stroid::IO::ViewMesh(mesh,
"Mesh Viewer - Colored by Element ID",
selected_mode,
glvis_host,
glvis_port);
exit(0);
}
if (*generate) {
if (config_filename.has_value()) {
cfg.load(config_filename.value());
}
const std::unique_ptr<mfem::Mesh> mesh = stroid::topology::BuildSkeleton(cfg);
stroid::topology::Finalize(*mesh, cfg);
stroid::topology::PromoteToHighOrder(*mesh, cfg);
stroid::topology::ProjectMesh(*mesh, cfg);
stroid::topology::OptimizeMesh(*mesh, cfg);
if (!no_save) {
const std::string& final_path = output_filename;
@@ -195,7 +242,7 @@ int main(int argc, char** argv) {
glvis_port);
}
} else if (!*info) {
std::println("Usage: {} [generate|info] --help", argv[0]);
std::println("Usage: {} [generate|info|view] --help", argv[0]);
}
return 0;