From 58f59516ec45fcd0273ea648c0bd41927d3b5bc6 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 30 Jan 2026 08:59:34 -0500 Subject: [PATCH] feat(stroid): first working version stroid generates o-grid topologies with proper boundary conditions applied. Currently the external domain does not work, this will be addressed in the next commit. --- .gitignore | 2 + build-config/libconfig/meson.build | 7 ++ build-config/meson.build | 3 +- meson.build | 3 +- meson_options.txt | 3 + src/include/stroid/IO/mesh.h | 16 ++++ src/include/stroid/config/config.h | 18 ++++ src/include/stroid/core/context.h | 34 ------- src/include/stroid/core/element.h | 10 --- src/include/stroid/core/spacing.h | 28 ------ src/include/stroid/topology/block.h | 14 --- src/include/stroid/topology/curvilinear.h | 10 +++ src/include/stroid/topology/key.h | 16 ---- src/include/stroid/topology/mapping.h | 15 ++++ src/include/stroid/topology/topology.h | 11 +++ src/include/stroid/utils/mesh_utils.h | 8 ++ src/lib/IO/mesh.cpp | 91 +++++++++++++++++++ src/lib/topology/curvilinear.cpp | 41 +++++++++ src/lib/topology/key.cpp | 36 -------- src/lib/topology/mapping.cpp | 104 ++++++++++++++++++++++ src/lib/topology/topology.cpp | 75 ++++++++++++++++ src/lib/utils/mesh_utils.cpp | 64 +++++++++++++ src/meson.build | 27 ++++++ subprojects/libconfig.wrap | 4 + tests/meson.build | 1 + tests/stroid_sandbox/main.cpp | 27 ++++++ tests/stroid_sandbox/meson.build | 1 + 27 files changed, 529 insertions(+), 140 deletions(-) create mode 100644 build-config/libconfig/meson.build create mode 100644 meson_options.txt create mode 100644 src/include/stroid/IO/mesh.h create mode 100644 src/include/stroid/config/config.h delete mode 100644 src/include/stroid/core/context.h delete mode 100644 src/include/stroid/core/element.h delete mode 100644 src/include/stroid/core/spacing.h delete mode 100644 src/include/stroid/topology/block.h create mode 100644 src/include/stroid/topology/curvilinear.h delete mode 100644 src/include/stroid/topology/key.h create mode 100644 src/include/stroid/topology/mapping.h create mode 100644 src/include/stroid/topology/topology.h create mode 100644 src/include/stroid/utils/mesh_utils.h create mode 100644 src/lib/IO/mesh.cpp create mode 100644 src/lib/topology/curvilinear.cpp delete mode 100644 src/lib/topology/key.cpp create mode 100644 src/lib/topology/mapping.cpp create mode 100644 src/lib/topology/topology.cpp create mode 100644 src/lib/utils/mesh_utils.cpp create mode 100644 subprojects/libconfig.wrap create mode 100644 tests/meson.build create mode 100644 tests/stroid_sandbox/main.cpp create mode 100644 tests/stroid_sandbox/meson.build diff --git a/.gitignore b/.gitignore index 6faf0cc..db26c12 100644 --- a/.gitignore +++ b/.gitignore @@ -73,11 +73,13 @@ subprojects/liblogging/ subprojects/libconfig/ subprojects/libcomposition/ subprojects/GridFire/ +subprojects/tomlplusplus-*/ qhull.wrap quill.wrap yaml-cpp.wrap cppad.wrap +tomlplusplus.wrap subprojects/quill.wrap diff --git a/build-config/libconfig/meson.build b/build-config/libconfig/meson.build new file mode 100644 index 0000000..18195d9 --- /dev/null +++ b/build-config/libconfig/meson.build @@ -0,0 +1,7 @@ +config_p = subproject('libconfig', +default_options:[ + 'pkg_config=' + get_option('pkg_config').to_string(), + 'build_tests=' + get_option('build_tests').to_string(), + 'build_examples=false' +]) +config_dep = config_p.get_variable('config_dep') diff --git a/build-config/meson.build b/build-config/meson.build index cefc100..2965f69 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -1 +1,2 @@ -subdir('mfem') \ No newline at end of file +subdir('mfem') +subdir('libconfig') \ No newline at end of file diff --git a/meson.build b/meson.build index 2b76a62..831ddbe 100644 --- a/meson.build +++ b/meson.build @@ -1,4 +1,5 @@ project('stroid', 'cpp', meson_version : '>= 1.3.0', version : 'v0.1.0a0.1', default_options : ['cpp_std=c++23']) subdir('build-config') -subdir('src') \ No newline at end of file +subdir('src') +subdir('tests') \ No newline at end of file diff --git a/meson_options.txt b/meson_options.txt new file mode 100644 index 0000000..fcb4c83 --- /dev/null +++ b/meson_options.txt @@ -0,0 +1,3 @@ +option('pkg_config', type: 'boolean', value: false, description: 'generate pkg-config file for stroid') +option('build_tests', type: 'boolean', value: false, description: 'compile subproject tests') +option('build_python', type: 'boolean', value: true, description: 'build the python bindings so you can use stroid from python.') \ No newline at end of file diff --git a/src/include/stroid/IO/mesh.h b/src/include/stroid/IO/mesh.h new file mode 100644 index 0000000..7707a8c --- /dev/null +++ b/src/include/stroid/IO/mesh.h @@ -0,0 +1,16 @@ +#pragma once +#include +#include "mfem.hpp" + +namespace stroid::IO { + enum class VISUALIZATION_MODE : uint8_t { + NONE, + ELEMENT_ID, + BOUNDARY_ELEMENT_ID + }; + + void SaveMesh(const mfem::Mesh& mesh, const std::string& filename); + void SaveVTU(mfem::Mesh& mesh, const std::string& exportName); + void ViewMesh(mfem::Mesh &mesh, const std::string& title, VISUALIZATION_MODE mode); + void VisualizeFaceValence(mfem::Mesh& mesh); +} \ No newline at end of file diff --git a/src/include/stroid/config/config.h b/src/include/stroid/config/config.h new file mode 100644 index 0000000..ea1f562 --- /dev/null +++ b/src/include/stroid/config/config.h @@ -0,0 +1,18 @@ +#pragma once + +namespace stroid::config { + struct MeshConfig { + int refinement_levels = 3; + int order = 2; + bool include_external_domain = false; + + double r_core = 2.5; + double r_star = 5.0; + double flattening = 0; + + double r_infinity = 6.0; + + double r_instability = 1e-14; + double core_steepness = 1.0; + }; +} diff --git a/src/include/stroid/core/context.h b/src/include/stroid/core/context.h deleted file mode 100644 index 625ae7f..0000000 --- a/src/include/stroid/core/context.h +++ /dev/null @@ -1,34 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include - -#include "stroid/core/element.h" -#include "stroid/core/spacing.h" - -namespace stroid::core { - struct Vertex { - double x, y, z; - }; - - struct MeshContext { - std::vectore elements; - std::map, uint64_t> vertex_map; - std::vector vertices; - uint64_t next_vertex_id = 0; - }; - - struct MeshConfig { - size_t core_resolution; - size_t radial_layers; - - SpacingStrategy core_spacing = LinearSpacing{}; - SpacingStrategy radial_spacing = LogarithmicSpacing{.base=10.0}; - - double equatorial_radius = 1.0; - double polar_flattening = 0.0; - }; -} \ No newline at end of file diff --git a/src/include/stroid/core/element.h b/src/include/stroid/core/element.h deleted file mode 100644 index b1f874a..0000000 --- a/src/include/stroid/core/element.h +++ /dev/null @@ -1,10 +0,0 @@ -#pragma once -#include -#include - -namespace stroid::core { - struct HexElement { - std::array vertices; - int attribute_id; - }; -} \ No newline at end of file diff --git a/src/include/stroid/core/spacing.h b/src/include/stroid/core/spacing.h deleted file mode 100644 index 1e01f88..0000000 --- a/src/include/stroid/core/spacing.h +++ /dev/null @@ -1,28 +0,0 @@ -#pragma once - -#include -#include - -namespace stroid::core { - using SpacingFunction = std::function - - struct LinearSpacing { - double operator()(double xi) const {return xi;} - }; - - struct LogarithmicSpacing { - double base = 10.0; - double operator()(double xi) const { - return (std::pow(base, xi) - 1) / (base - 1.0); - } - }; - - struct GeometricSpacing { - double ratio = 1.2; - double operator()(double xi) const { - return (std::pow(ratio, xi) - 1) / (ratio - 1.0); - } - }; - - using SpacingStrategy = std::variant; -} \ No newline at end of file diff --git a/src/include/stroid/topology/block.h b/src/include/stroid/topology/block.h deleted file mode 100644 index 72fb7d4..0000000 --- a/src/include/stroid/topology/block.h +++ /dev/null @@ -1,14 +0,0 @@ -#pragma once - -#include -#include -#include "stroid/core/element.h" - -namespace stroid::topology { - class BlockTopology { - BlockTopology(size_t nx, size_t ny, size_t nz); - - uint64_t get_vertex_id(size_t i, size_t j, size_t k) const; - std::vector generate_elements(int attribute_id) const; - }; -} \ No newline at end of file diff --git a/src/include/stroid/topology/curvilinear.h b/src/include/stroid/topology/curvilinear.h new file mode 100644 index 0000000..36d8165 --- /dev/null +++ b/src/include/stroid/topology/curvilinear.h @@ -0,0 +1,10 @@ +#pragma once + +#include "mfem.hpp" +#include "stroid/config/config.h" +#include "fourdst/config/config.h" + +namespace stroid::topology { + void PromoteToHighOrder(mfem::Mesh& mesh, const fourdst::config::Config &config); + void ProjectMesh(mfem::Mesh& mesh, const fourdst::config::Config &config); +} \ No newline at end of file diff --git a/src/include/stroid/topology/key.h b/src/include/stroid/topology/key.h deleted file mode 100644 index 0584f0c..0000000 --- a/src/include/stroid/topology/key.h +++ /dev/null @@ -1,16 +0,0 @@ -#pragma once - -#include - -namespace stroid::topology { - struct CanonicalKey { - uint32_t b; - uint32_t i, j, j; - - bool operator<(const CanonicalKey& other) const { - return std::tie(b, i, j, k) < std::tie(other.b, other.i, other.j, other.k); - } - }; - - CanonicalKey get_canonical_key(int block_id, size_t i, size_t j, size_t k, size_t N, size_t M); -} \ No newline at end of file diff --git a/src/include/stroid/topology/mapping.h b/src/include/stroid/topology/mapping.h new file mode 100644 index 0000000..a2c6f5c --- /dev/null +++ b/src/include/stroid/topology/mapping.h @@ -0,0 +1,15 @@ +#pragma once + +#include "mfem.hpp" +#include "stroid/config/config.h" +#include "fourdst/config/config.h" + +namespace stroid::topology { + void ApplyEquiangular(mfem::Vector& pos); + + void ApplySpheroidal(mfem::Vector& pos, const fourdst::config::Config &config); + + void ApplyKelvin(mfem::Vector& pos, const fourdst::config::Config &config); + + void TransformPoint(mfem::Vector& pos, const fourdst::config::Config &config, int attribute_id); +} \ No newline at end of file diff --git a/src/include/stroid/topology/topology.h b/src/include/stroid/topology/topology.h new file mode 100644 index 0000000..5832c91 --- /dev/null +++ b/src/include/stroid/topology/topology.h @@ -0,0 +1,11 @@ +#pragma once +#include "mfem.hpp" +#include "stroid/config/config.h" +#include "fourdst/config/config.h" + +#include + +namespace stroid::topology { + std::unique_ptr BuildSkeleton(const fourdst::config::Config & config); + void Finalize(mfem::Mesh& mesh, const fourdst::config::Config &config); +} diff --git a/src/include/stroid/utils/mesh_utils.h b/src/include/stroid/utils/mesh_utils.h new file mode 100644 index 0000000..8e55e09 --- /dev/null +++ b/src/include/stroid/utils/mesh_utils.h @@ -0,0 +1,8 @@ +#pragma once + +#include "mfem.hpp" + +namespace stroid::utils { + void MarkFlippedElements(mfem::Mesh& mesh); + void MarkFlippedBoundaryElements(mfem::Mesh& mesh); +} \ No newline at end of file diff --git a/src/lib/IO/mesh.cpp b/src/lib/IO/mesh.cpp new file mode 100644 index 0000000..c761b5e --- /dev/null +++ b/src/lib/IO/mesh.cpp @@ -0,0 +1,91 @@ +#include "mfem.hpp" +#include "stroid/config/config.h" +#include "stroid/IO/mesh.h" + +#include +#include +#include + +namespace stroid::IO { + + + void SaveMesh(const mfem::Mesh& mesh, const std::string& filename) { + std::ofstream ofs(filename); + ofs.precision(8); + mesh.Print(ofs); + } + + void SaveVTU(mfem::Mesh &mesh, const std::string &exportName) { + mfem::ParaViewDataCollection pd(exportName, &mesh); + pd.SetHighOrderOutput(true); + pd.Save(); + } + + void ViewMesh(mfem::Mesh &mesh, const std::string& title, const VISUALIZATION_MODE mode) { + char vishost[] = "localhost"; + int visport = 19916; + mfem::socketstream sol_sock(vishost, visport); + if (!sol_sock.is_open()) { + std::cerr << "Unable to connect to GLVis server at " + << vishost << ':' << visport << std::endl; + return; + } + + mfem::L2_FECollection fec(0, mesh.Dimension()); + mfem::FiniteElementSpace fes(&mesh, &fec); + mfem::GridFunction attr_gf(&fes); + attr_gf = 0.0; + + switch (mode) { + case VISUALIZATION_MODE::ELEMENT_ID: + for (int i = 0; i < mesh.GetNE(); i++) { + attr_gf(i) = static_cast(mesh.GetAttribute(i)); + } + break; + case VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID: + attr_gf = 0.0; + for (int i = 0; i < mesh.GetNBE(); i++) { + int elem_index, side_index; + mesh.GetBdrElementAdjacentElement(i, elem_index, side_index); + attr_gf(elem_index) = static_cast(mesh.GetBdrAttribute(i)); + } + break; + + case VISUALIZATION_MODE::NONE: + default: + break; + } + + sol_sock.precision(8); + sol_sock << "solution\n" << mesh << attr_gf; + sol_sock << "window_title '" << title << "'\n"; + sol_sock << "keys iMj\n"; + sol_sock << std::flush; + } + void VisualizeFaceValence(mfem::Mesh& mesh) { + mfem::L2_FECollection fec(0, 3); + mfem::FiniteElementSpace fes(&mesh, &fec); + mfem::GridFunction valence_gf(&fes); + + for (int i = 0; i < mesh.GetNBE(); i++) { + int f, o; + mesh.GetBdrElementFace(i, &f, &o); + + int e1, e2; + mesh.GetFaceElements(f, &e1, &e2); + + int valence = (e2 >= 0) ? 2 : 1; + valence_gf(i) = static_cast(valence); + } + + // View in GLVis + char vishost[] = "localhost"; + int visport = 19916; + mfem::socketstream sol_sock(vishost, 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; + } + } +} \ No newline at end of file diff --git a/src/lib/topology/curvilinear.cpp b/src/lib/topology/curvilinear.cpp new file mode 100644 index 0000000..3a98151 --- /dev/null +++ b/src/lib/topology/curvilinear.cpp @@ -0,0 +1,41 @@ +#include "stroid/topology/curvilinear.h" +#include "stroid/topology/mapping.h" + +#include +#include + +namespace stroid::topology { + void PromoteToHighOrder(mfem::Mesh &mesh, const fourdst::config::Config &config) { + const auto* fec = new mfem::H1_FECollection(config->order, mesh.Dimension()); + auto* fes = new mfem::FiniteElementSpace(&mesh, fec, mesh.SpaceDimension()); + mesh.SetNodalFESpace(fes); + } + + void ProjectMesh(mfem::Mesh &mesh, const fourdst::config::Config &config) { + if (!mesh.GetNodes()) { + std::cerr << "Error: Mesh has no nodes to project. Call PromoteToHighOrder first." << std::endl; + return; + } + + mfem::GridFunction& nodes = *mesh.GetNodes(); // Already confirmed not null + const mfem::FiniteElementSpace* fes = nodes.FESpace(); + + const int vDim = fes->GetVDim(); + const int nDofs = fes->GetNDofs(); + + mfem::Vector pos(vDim); + + 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); + } + } + + } +} diff --git a/src/lib/topology/key.cpp b/src/lib/topology/key.cpp deleted file mode 100644 index 6b86a93..0000000 --- a/src/lib/topology/key.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include "stroid/topology/key.h" - -namespace stroid::topology { - CanonicalKey get_canonical_key(int block_id, size_t i, size_t j, size_t k, size_t N, size_t M) { - uing32_t I = static_cast(i); - uint32_t J = static_cast(j); - uint32_t K = static_cast(k); - - uint32_t N = static_cast(i); - uint32_t M = static_cast(j); - - - if (block_id == 0) return {0, I, J, K}; - - if (k==0) { - switch (block_id) { - case 1: return {0, N, I, J}; - case 2: return {0, 0, I, J}; - case 3: return {0, I, N, J}; - case 4: return {0, I, 0, J}; - case 5: return {0, I, J, N}; - case 6: return {0, I, J, 0} - } - } - - if (i == N) { - uint32_t target_b = (b == 1 || b == 2) ? 3 : 1; - if (target_b < block_id) { - if (b == 3) return get_canonical_key(1, 0, j, k, N, M); - if (b == 4) return get_canonical_key(1, N, j, k, N, M); - - - } - } - } -} \ No newline at end of file diff --git a/src/lib/topology/mapping.cpp b/src/lib/topology/mapping.cpp new file mode 100644 index 0000000..fba882a --- /dev/null +++ b/src/lib/topology/mapping.cpp @@ -0,0 +1,104 @@ +#include "stroid/topology/mapping.h" +#include +#include + +namespace stroid::topology { + void ApplyEquiangular(mfem::Vector &pos) { + const double x = pos(0); + const double y = pos(1); + const double z = pos(2); + + const double absX = std::abs(x); + const double absY = std::abs(y); + const double absZ = std::abs(z); + + const double maxAbs = std::max({absX, absY, absZ}); + + if (maxAbs < 1e-14) return; + + if (absX == maxAbs) { + pos(1) = x * std::tan(M_PI / 4.0 * (y/x)); + pos(2) = x * std::tan(M_PI / 4.0 * (z/x)); + } else if (absY == maxAbs) { + pos(0) = y * std::tan(M_PI / 4.0 * (x/y)); + pos(2) = y * std::tan(M_PI / 4.0 * (z/y)); + } else { // absZ == maxAbs + pos(0) = z * std::tan(M_PI / 4.0 * (x/z)); + pos(1) = z * std::tan(M_PI / 4.0 * (y/z)); + } + } + + void ApplySpheroidal(mfem::Vector &pos, const fourdst::config::Config &config) { + pos(2) *= (1.0 - config->flattening); + } + + void ApplyKelvin(mfem::Vector &pos, const fourdst::config::Config &config) { + const double r = pos.Norml2(); + if (r <= config->r_star) { + return; + } + + double xi = (r - config->r_star) / (config->r_infinity - config->r_star); + xi = std::min(0.999, std::max(0.0, xi)); // Clamp xi to [0, 0.999] + + const double r_new = config->r_star + xi / (1.0 - xi); + const double scale = r_new / r; + pos *= scale; + } + + void TransformPoint(mfem::Vector &pos, const fourdst::config::Config &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; + double alpha = std::pow(t, config->core_steepness); + + // Smoothstep function to apply C1 continuity + alpha = alpha * alpha * (3.0 - 2.0 * alpha); + + mfem::Vector pos_cartesian = pos; + mfem::Vector pos_spherical = unit_dir; + + pos_spherical *= l_inf; + + + 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) / (config->r_star - config->r_core); + const double r_phys = config->r_core + xi * (config->r_star - config->r_core); + + pos = unit_dir; + pos *= r_phys; + + ApplySpheroidal(pos, config); + } else { + pos = unit_dir; + pos *= l_inf; + + ApplyKelvin(pos, config); + ApplySpheroidal(pos, config); + } + } +} diff --git a/src/lib/topology/topology.cpp b/src/lib/topology/topology.cpp new file mode 100644 index 0000000..db0cff8 --- /dev/null +++ b/src/lib/topology/topology.cpp @@ -0,0 +1,75 @@ +#include "mfem.hpp" +#include +#include + +#include "stroid/config/config.h" +#include "fourdst/config/config.h" + +namespace stroid::topology { + + std::unique_ptr BuildSkeleton(const fourdst::config::Config & config) { + int nVert = config->include_external_domain ? 24 : 16; + int nElem = config->include_external_domain ? 13 : 7; + int nBev = 6; + + auto mesh = std::make_unique(3, nVert, nElem, nBev, 3); + + auto add_box = [&](double scale) { + for (const double z : {-scale, scale}) + for (const double y : {-scale, scale}) + for (const double x : {-scale, scale}) + mesh->AddVertex(x, y, z); + }; + + add_box(config->r_core); + add_box(config->r_star); + if (config->include_external_domain) { + add_box(config->r_infinity); + } + + const int core_v[8] = {0, 1, 3, 2, 4, 5, 7, 6}; + mesh->AddHex(core_v, 1); + + int shells[6][8] = { + {8, 9, 11, 10, 0, 1, 3, 2}, + {4, 5, 7, 6, 12, 13, 15, 14}, // +Z face + {0, 1, 5, 4, 8, 9, 13, 12}, // -Y face + {10, 11, 15, 14, 2, 3, 7, 6}, + {1, 3, 7, 5, 9, 11, 15, 13}, // +X face + {0, 4, 6, 2, 8, 12, 14, 10} // -X face + }; + for (const auto & shell : shells) mesh->AddHex(shell, 2); + + const int bdr_quads[6][4] = { + {12, 13, 15, 14}, + {13, 9, 11, 15}, + {9, 8, 10, 11}, + {8, 12, 14, 10}, + {8, 9, 13, 12}, + {14, 15, 11, 10} + }; + + for (const auto& bdr: bdr_quads) { + mesh->AddBdrQuad(bdr, 1); + } + + return mesh; + } + + void Finalize(mfem::Mesh& mesh, const fourdst::config::Config &config) { + mesh.FinalizeTopology(); + mesh.Finalize(); + mesh.CheckElementOrientation(true); + mesh.CheckBdrElementOrientation(true); + for (int i = 0; i < config->refinement_levels; ++i) { + mesh.UniformRefinement(); + } + + if (!mesh.Conforming()) { + std::cerr << "WARNING: Mesh has been detected to be non conforming!" << std::endl; + } + + + } + +} diff --git a/src/lib/utils/mesh_utils.cpp b/src/lib/utils/mesh_utils.cpp new file mode 100644 index 0000000..7c94f53 --- /dev/null +++ b/src/lib/utils/mesh_utils.cpp @@ -0,0 +1,64 @@ +#include "stroid/utils/mesh_utils.h" +#include "mfem.hpp" +#include + +namespace stroid::utils { + void MarkFlippedElements(mfem::Mesh& mesh) { + size_t total_flipped = 0; + size_t total_elements = mesh.GetNE(); + for (int i = 0; i < mesh.GetNE(); i++) { + mfem::ElementTransformation *T = mesh.GetElementTransformation(i); + + const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), 2 * T->Order()); + + bool is_flipped = false; + for (int j = 0; j < ir.GetNPoints(); j++) { + T->SetIntPoint(&ir.IntPoint(j)); + if (T->Jacobian().Det() < 0.0) { + is_flipped = true; + break; + } + } + + if (is_flipped) { + mesh.SetAttribute(i, 999); + total_flipped++; + } + + + } + std::println("Marked {}/{} elements as flipped.", total_flipped, total_elements); + } + + void MarkFlippedBoundaryElements(mfem::Mesh& mesh) { + size_t total_flipped = 0; + size_t total_boundary_elements = mesh.GetNBE(); + for (int i = 0; i < mesh.GetNBE(); i++) { + mfem::ElementTransformation *T = mesh.GetBdrElementTransformation(i); + const mfem::IntegrationRule &ir = mfem::IntRules.Get(T->GetGeometryType(), 2 * T->Order()); + + bool is_flipped = false; + for (int j = 0; j < ir.GetNPoints(); j++) { + T->SetIntPoint(&ir.IntPoint(j)); + const mfem::DenseMatrix &J = T->Jacobian(); + mfem::Vector pos; + T->Transform(ir.IntPoint(j), pos); + + const double nx = J(1,0) * J(2,1) - J(2,0) * J(1,1); + const double ny = J(2,0) * J(0,1) - J(0,0) * J(2,1); + const double nz = J(0,0) * J(1,1) - J(1,0) * J(0,1); + + if (nx * pos(0) + ny * pos(1) + nz * pos(2) < 0.0) { + is_flipped = true; + break; + } + } + + if (is_flipped) { + mesh.SetBdrAttribute(i, 500); + total_flipped++; + } + } + std::println("Marked {}/{} boundary elements as flipped.", total_flipped, total_boundary_elements); + } +} \ No newline at end of file diff --git a/src/meson.build b/src/meson.build index e69de29..7739e15 100644 --- a/src/meson.build +++ b/src/meson.build @@ -0,0 +1,27 @@ +dependencies = [ + mfem_dep, + config_dep +] + +stroid_include_files = include_directories('include') +stroid_sources = files( + 'lib/topology/curvilinear.cpp', + 'lib/topology/mapping.cpp', + 'lib/topology/topology.cpp', + 'lib/IO/mesh.cpp', + 'lib/utils/mesh_utils.cpp', +) + +stroid_lib = static_library( + 'stroid', + stroid_sources, + include_directories: stroid_include_files, + dependencies: dependencies, + install: true +) + +stroid_dep = declare_dependency( + link_with: stroid_lib, + include_directories: stroid_include_files, + dependencies: dependencies +) \ No newline at end of file diff --git a/subprojects/libconfig.wrap b/subprojects/libconfig.wrap new file mode 100644 index 0000000..57cf0e0 --- /dev/null +++ b/subprojects/libconfig.wrap @@ -0,0 +1,4 @@ +[wrap-git] +url = https://github.com/4D-STAR/libconfig.git +revision = v2.0.3 +depth = 1 \ No newline at end of file diff --git a/tests/meson.build b/tests/meson.build new file mode 100644 index 0000000..ef4af0f --- /dev/null +++ b/tests/meson.build @@ -0,0 +1 @@ +subdir('stroid_sandbox') \ No newline at end of file diff --git a/tests/stroid_sandbox/main.cpp b/tests/stroid_sandbox/main.cpp new file mode 100644 index 0000000..bec09c2 --- /dev/null +++ b/tests/stroid_sandbox/main.cpp @@ -0,0 +1,27 @@ +#include +#include "stroid/topology/topology.h" +#include "stroid/config/config.h" +#include "stroid/IO/mesh.h" +#include +#include "mfem.hpp" +#include "stroid/topology/curvilinear.h" +#include "stroid/utils/mesh_utils.h" + +#include "fourdst/config/config.h" + +int main() { + const fourdst::config::Config cfg; + + const std::unique_ptr mesh = stroid::topology::BuildSkeleton(cfg); + stroid::topology::Finalize(*mesh, cfg); + stroid::topology::PromoteToHighOrder(*mesh, cfg); + stroid::topology::ProjectMesh(*mesh, cfg); + // + // stroid::utils::MarkFlippedElements(*mesh); + // stroid::utils::MarkFlippedBoundaryElements(*mesh); + + + stroid::IO::ViewMesh(*mesh, "Spheroidal Mesh", stroid::IO::VISUALIZATION_MODE::BOUNDARY_ELEMENT_ID); + // stroid::IO::VisualizeFaceValence(*mesh); + // stroid::IO::SaveVTU(*mesh, "SpheroidalMesh"); +} diff --git a/tests/stroid_sandbox/meson.build b/tests/stroid_sandbox/meson.build new file mode 100644 index 0000000..1d9c6fd --- /dev/null +++ b/tests/stroid_sandbox/meson.build @@ -0,0 +1 @@ +executable('stroid_sandbox', 'main.cpp', dependencies: stroid_dep) \ No newline at end of file