diff --git a/.gitignore b/.gitignore index 07227a2..3d33190 100644 --- a/.gitignore +++ b/.gitignore @@ -56,4 +56,9 @@ tags *.log *.cache *.private -*.private/ \ No newline at end of file +*.private/ + +subprojects/mfem/ +subprojects/tetgen/ + +.vscode/ diff --git a/Readme.md b/Readme.md index 93b9b91..e1beda1 100644 --- a/Readme.md +++ b/Readme.md @@ -1,8 +1,16 @@ # New implimentation of 3+1D SSE +New (as yet unnamed) 4DSSE code. We need an exciting name. -# Building -In order to build run +This code is very early in development and should not be used for scientific purposes yet. + + +## Building +In order to build you will need meson installed on your system. The easiest way to do this is to use the python package manager (pip) +```bash +pip install meson +``` +You can then either use the mk script or meson commands automatically. When running either the script or meson commands manually `MFEM` will be pulled from github and built. As part of this a small patch will be applied to the MFEM `CMakeLists.txt` file. This process should only need to happen once as future builds will use the cached version of MFEM in `subprojects` and the cached build files of `MFEM` in `build`. ```bash ./mk ``` @@ -10,3 +18,6 @@ if you want to build with no test suite run ```bash ./mk --noTest ``` + +## Current Status +Currently we are working on implimenting modules such as opacity, equation of state, polytrope, and meshing. Builds may not work on any branches at any time. diff --git a/build-config/meson.build b/build-config/meson.build new file mode 100644 index 0000000..cefc100 --- /dev/null +++ b/build-config/meson.build @@ -0,0 +1 @@ +subdir('mfem') \ No newline at end of file diff --git a/build-config/mfem/disable_mfem_selfcheck.patch b/build-config/mfem/disable_mfem_selfcheck.patch new file mode 100644 index 0000000..8e438c0 --- /dev/null +++ b/build-config/mfem/disable_mfem_selfcheck.patch @@ -0,0 +1,41 @@ +--- subprojects/mfem/CMakeLists.txt 2025-02-12 15:54:52.454728232 -0500 ++++ CMakeLists.txt.bak 2025-02-12 16:08:06.654542689 -0500 +@@ -765,7 +765,7 @@ + if (MFEM_ENABLE_EXAMPLES) + add_subdirectory(examples) #install examples if enabled + else() +- add_subdirectory(examples EXCLUDE_FROM_ALL) ++ # add_subdirectory(examples EXCLUDE_FROM_ALL) + endif() + + # Create a target for all miniapps and, optionally, enable it. +@@ -774,7 +774,7 @@ + if (MFEM_ENABLE_MINIAPPS) + add_subdirectory(miniapps) #install miniapps if enabled + else() +- add_subdirectory(miniapps EXCLUDE_FROM_ALL) ++ # add_subdirectory(miniapps EXCLUDE_FROM_ALL) + endif() + + # Target to build all executables, i.e. everything. +@@ -801,19 +801,7 @@ + add_dependencies(${MFEM_EXEC_PREREQUISITES_TARGET_NAME} copy_data) + endif() + +-# Add 'check' target - quick test +-set(MFEM_CHECK_TARGET_NAME ${MFEM_CUSTOM_TARGET_PREFIX}check) +-if (NOT MFEM_USE_MPI) +- add_custom_target(${MFEM_CHECK_TARGET_NAME} +- ${CMAKE_CTEST_COMMAND} -R \"^ex1_ser\" -C ${CMAKE_CFG_INTDIR} +- USES_TERMINAL) +- add_dependencies(${MFEM_CHECK_TARGET_NAME} ex1) +-else() +- add_custom_target(${MFEM_CHECK_TARGET_NAME} +- ${CMAKE_CTEST_COMMAND} -R \"^ex1p\" -C ${CMAKE_CFG_INTDIR} +- USES_TERMINAL) +- add_dependencies(${MFEM_CHECK_TARGET_NAME} ex1p) +-endif() ++message(STATUS "MFEM Miniapps and Examples disabled by patch!") + + #------------------------------------------------------------------------------- + # Documentation diff --git a/build-config/mfem/meson.build b/build-config/mfem/meson.build new file mode 100644 index 0000000..84f69c8 --- /dev/null +++ b/build-config/mfem/meson.build @@ -0,0 +1,24 @@ +cmake = import('cmake') +patchFile = files('disable_mfem_selfcheck.patch') + +patch_check = run_command('grep', '-q', 'MFEM_CHECK_TARGET_NAME', 'subprojects/mfem/CMakeLists.txt', check: false) +if patch_check.returncode() == 0 + message('Patching MFEM CMakeLists.txt to remove self checks') + run_command('patch', '-p4', '-d', '../../subprojects/mfem', '-i', patchFile[0].full_path(), check: true) +else + message('MFEM CMakeLists.txt already patched') +endif + +mfem_cmake_options = cmake.subproject_options() +mfem_cmake_options.add_cmake_defines({ + 'MFEM_ENABLE_EXAMPLES': 'OFF', + 'MFEM_ENABLE_TESTING': 'OFF', + 'MFEM_ENABLE_MINIAPPS': 'OFF', + 'MFEM_USE_BENCMARK': 'OFF', +}) + +mfem_sp = cmake.subproject( + 'mfem', + options: mfem_cmake_options) +mfem_dep = mfem_sp.dependency('mfem') +add_project_arguments('-I' + meson.current_build_dir() + '/subprojects/mfem/__CMake_build/config', language: 'cpp') \ No newline at end of file diff --git a/meson.build b/meson.build index d6ad9fe..3724fac 100644 --- a/meson.build +++ b/meson.build @@ -3,7 +3,10 @@ project('4DSSE', 'cpp', version: '0.0.1a', default_options: ['cpp_std=c++23'], m # Add default visibility for all C++ targets add_project_arguments('-fvisibility=default', language: 'cpp') +# Build external dependencies first so that all the embedded resources are available to the other targets +subdir('build-config') +# Build the main project subdir('src') if get_option('build_tests') subdir('tests') diff --git a/mk b/mk index 6e3a90f..9bab354 100755 --- a/mk +++ b/mk @@ -1,10 +1,5 @@ #!/bin/bash -# Check if the build directory is present, and remove it -if [ -d "build" ]; then - rm -rf build -fi - # Check for the --noTest flag if [[ "$1" == "--noTest" ]]; then meson setup build -Dbuild_tests=false diff --git a/src/meshIO/meson.build b/src/meshIO/meson.build new file mode 100644 index 0000000..1ea870a --- /dev/null +++ b/src/meshIO/meson.build @@ -0,0 +1,19 @@ +# Define the library +meshIO_sources = files( + 'private/meshIO.cpp', +) + +meshIO_headers = files( + 'public/meshIO.h' +) + +# Define the libmeshIO library so it can be linked against by other parts of the build system +libmeshIO = static_library('meshIO', + meshIO_sources, + include_directories: include_directories('public'), + cpp_args: ['-fvisibility=default'], + dependencies: [mfem_dep], + install : true) + +# Make headers accessible +install_headers(meshIO_headers, subdir : '4DSSE/meshIO') \ No newline at end of file diff --git a/src/meshIO/private/meshIO.cpp b/src/meshIO/private/meshIO.cpp new file mode 100644 index 0000000..8674af3 --- /dev/null +++ b/src/meshIO/private/meshIO.cpp @@ -0,0 +1,37 @@ +#include "mfem.hpp" +#include +#include +#include + +#include "meshIO.h" + + +MeshIO::MeshIO(const std::string &mesh_file) +{ + mesh_file_ = mesh_file; + std::ifstream mesh_stream(mesh_file); + if (!mesh_stream) + { + throw std::runtime_error("Mesh file not found: " + mesh_file); + loaded_ = false; + } + else + { + mesh_ = mfem::Mesh(mesh_stream, 1, 2); + loaded_ = true; + } +} + +MeshIO::~MeshIO() +{ +} + +bool MeshIO::IsLoaded() const +{ + return loaded_; +} + +mfem::Mesh& MeshIO::GetMesh() +{ + return mesh_; +} \ No newline at end of file diff --git a/src/meshIO/public/meshIO.h b/src/meshIO/public/meshIO.h new file mode 100644 index 0000000..5975801 --- /dev/null +++ b/src/meshIO/public/meshIO.h @@ -0,0 +1,42 @@ +#ifndef MESHIO_H +#define MESHIO_H + +#include "mfem.hpp" +#include + +/** + * @brief Class for handling mesh input/output operations. + */ +class MeshIO +{ +private: + bool loaded_; ///< Flag to indicate if the mesh is loaded + std::string mesh_file_; ///< Filename of the mesh file + mfem::Mesh mesh_; ///< The mesh object + +public: + /** + * @brief Constructor that initializes the MeshIO object with a mesh file. + * @param mesh_file The name of the mesh file. + */ + MeshIO(const std::string &mesh_file); + + /** + * @brief Destructor for the MeshIO class. + */ + ~MeshIO(); + + /** + * @brief Get the mesh object. + * @return Reference to the mesh object. + */ + mfem::Mesh& GetMesh(); + + /** + * @brief Check if the mesh is loaded. + * @return True if the mesh is loaded, false otherwise. + */ + bool IsLoaded() const; +}; + +#endif // MESHIO_H \ No newline at end of file diff --git a/src/meson.build b/src/meson.build index db87771..9ae232d 100644 --- a/src/meson.build +++ b/src/meson.build @@ -4,4 +4,5 @@ subdir('resources') # Build the main source code subdir('dobj') subdir('const') -subdir('opatIO') \ No newline at end of file +subdir('opatIO') +subdir('meshIO') \ No newline at end of file diff --git a/src/poly/coeff/meson.build b/src/poly/coeff/meson.build new file mode 100644 index 0000000..112a1ff --- /dev/null +++ b/src/poly/coeff/meson.build @@ -0,0 +1,23 @@ +polyCoeff_sources = files( + 'private/coeff.cpp' +) + +polyCoeff_headers = files( + 'public/coeff.h' +) + +libPolyCoeff = static_library('polyCoeff', + polyCoeff_sources, + include_directories : include_directories('.'), + cpp_args: ['-fvisibility=default'], + dependencies: [mfem_dep], + install: true +) + + +polyCoeff_dep = declare_dependency( + include_directories : include_directories('.'), + link_with : libPolyCoeff, + sources : polyCoeff_sources, + dependencies : [mfem_dep] +) \ No newline at end of file diff --git a/src/poly/coeff/private/polyCoeff.cpp b/src/poly/coeff/private/polyCoeff.cpp new file mode 100644 index 0000000..c4076a0 --- /dev/null +++ b/src/poly/coeff/private/polyCoeff.cpp @@ -0,0 +1,40 @@ +#include "mfem.hpp" +#include + +#include "coeff.h" + +/** + * @brief Computes the xi coefficient function. + * + * @param x Input vector. + * @return double The computed xi coefficient. + */ +double xi_coeff_func(const mfem::Vector &x) +{ + return std::pow(x(0), 2); +} + +/** + * @brief Computes the vector xi coefficient function. + * + * @param x Input vector. + * @param v Output vector to store the computed xi coefficient. + */ +void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v) +{ + v.SetSize(1); + v[0] = -std::pow(x(0), 2); +} + +/** + * @brief Computes the initial guess for theta. + * + * @param x Input vector. + * @param root Root value used in the computation. + * @return double The initial guess for theta. + */ +double theta_initial_guess(const mfem::Vector &x, double root) +{ + double xi = x[0]; + return 1 - std::pow(xi / root, 2); +} \ No newline at end of file diff --git a/src/poly/coeff/public/polyCoeff.h b/src/poly/coeff/public/polyCoeff.h new file mode 100644 index 0000000..46204f1 --- /dev/null +++ b/src/poly/coeff/public/polyCoeff.h @@ -0,0 +1,8 @@ +#include "mfem.hpp" +#include + +double xi_coeff_func(const mfem::Vector &x); + +void vec_xi_coeff_func(const mfem::Vector &x, mfem::Vector &v); + +double theta_initial_guess(const mfem::Vector &x, double root); \ No newline at end of file diff --git a/src/poly/meson.build b/src/poly/meson.build new file mode 100644 index 0000000..d9c272f --- /dev/null +++ b/src/poly/meson.build @@ -0,0 +1,3 @@ +subdir('coeff') +subdir('utils') +subdir('solver') \ No newline at end of file diff --git a/src/poly/solver/meson.build b/src/poly/solver/meson.build new file mode 100644 index 0000000..e69de29 diff --git a/src/poly/utils/meson.build b/src/poly/utils/meson.build new file mode 100644 index 0000000..fc072eb --- /dev/null +++ b/src/poly/utils/meson.build @@ -0,0 +1,24 @@ +polyutils_sources = files( + 'private/polyIO.cpp', + 'private/polyMFEMUtils.cpp' +) + +polyutils_headers = files( + 'public/polyIO.h', + 'public/polyMFEMUtils.h' +) + +libpolyutils = static_library('polyutils', + polyutils_sources, + include_directories : include_directories('.'), + cpp_args: ['-fvisibility=default'], + dependencies: [mfem_dep], + install: true +) + +libpolyutils_dep = declare_dependency( + include_directories : include_directories('.'), + link_with : libpolyutils, + sources : polyutils_sources, + dependencies : [mfem_dep] +) diff --git a/src/poly/utils/private/polyIO.cpp b/src/poly/utils/private/polyIO.cpp new file mode 100644 index 0000000..d540cb0 --- /dev/null +++ b/src/poly/utils/private/polyIO.cpp @@ -0,0 +1,23 @@ +#include "mfem.hpp" +#include +#include + +#include "polyIO.h" + +void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename) { + std::ofstream file(filename); + if (!file.is_open()) { + std::cerr << "Error: Could not open " << filename << " for writing." << std::endl; + return; + } + + file << "xi,u\n"; // CSV header + + for (int i = 0; i < u.Size(); i++) { + double xi = mesh.GetVertex(i)[0]; // Get spatial coordinate + file << xi << "," << u[i] << "\n"; // Write to CSV + } + + file.close(); + std::cout << "Solution written to " << filename << std::endl; +} \ No newline at end of file diff --git a/src/poly/utils/private/polyMFEMUtils.cpp b/src/poly/utils/private/polyMFEMUtils.cpp new file mode 100644 index 0000000..781efab --- /dev/null +++ b/src/poly/utils/private/polyMFEMUtils.cpp @@ -0,0 +1,175 @@ +#include "mfem.hpp" +#include +#include +#include + +#include "polyMFEMUtils.h" + +NonlinearPowerIntegrator::NonlinearPowerIntegrator( + mfem::FunctionCoefficient &coeff, + double n) : coeff_(coeff), polytropicIndex(n) { + +} + +void NonlinearPowerIntegrator::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + int dof = el.GetDof(); + elvect.SetSize(dof); + elvect = 0.0; + + mfem::Vector shape(dof); + + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + mfem::IntegrationPoint ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + double weight = ip.weight * Trans.Weight(); + + el.CalcShape(ip, shape); + + double u_val = 0.0; + for (int j = 0; j < dof; j++) { + u_val += elfun(j) * shape(j); + } + double u_safe = std::max(u_val, 0.0); + double u_nl = std::pow(u_safe, polytropicIndex); + + double coeff_val = coeff_.Eval(Trans, ip); + double x2_u_nl = coeff_val * u_nl; + + for (int i = 0; i < dof; i++){ + elvect(i) += shape(i) * x2_u_nl * weight; + } + } +} + +void NonlinearPowerIntegrator::AssembleElementGrad ( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + + const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3); + int dof = el.GetDof(); + elmat.SetSize(dof); + elmat = 0.0; + mfem::Vector shape(dof); + + for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) { + mfem::IntegrationPoint ip = ir->IntPoint(iqp); + Trans.SetIntPoint(&ip); + double weight = ip.weight * Trans.Weight(); + + el.CalcShape(ip, shape); + + double u_val = 0.0; + + for (int j = 0; j < dof; j++) { + u_val += elfun(j) * shape(j); + } + double coeff_val = coeff_.Eval(Trans, ip); + + + // Calculate the Jacobian + double u_safe = std::max(u_val, 0.0); + double d_u_nl = coeff_val * polytropicIndex * std::pow(u_safe, polytropicIndex - 1); + double x2_d_u_nl = d_u_nl; + + for (int i = 0; i < dof; i++) { + for (int j = 0; j < dof; j++) { + elmat(i, j) += shape(i) * x2_d_u_nl * shape(j) * weight; + } + } + + } +} + +BilinearIntegratorWrapper::BilinearIntegratorWrapper( + mfem::BilinearFormIntegrator *integratorInput + ) : integrator(integratorInput) { } + +BilinearIntegratorWrapper::~BilinearIntegratorWrapper() { + delete integrator; +} + +void BilinearIntegratorWrapper::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + int dof = el.GetDof(); + mfem::DenseMatrix elMat(dof); + integrator->AssembleElementMatrix(el, Trans, elMat); + elvect.SetSize(dof); + elvect = 0.0; + for (int i = 0; i < dof; i++) + { + double sum = 0.0; + for (int j = 0; j < dof; j++) + { + sum += elMat(i, j) * elfun(j); + } + elvect(i) = sum; + } +} + +void BilinearIntegratorWrapper::AssembleElementGrad(const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + int dof = el.GetDof(); + elmat.SetSize(dof, dof); + elmat = 0.0; + integrator->AssembleElementMatrix(el, Trans, elmat); +} + +CompositeNonlinearIntegrator::CompositeNonlinearIntegrator() { } + + +CompositeNonlinearIntegrator::~CompositeNonlinearIntegrator() { + for (size_t i = 0; i < integrators.size(); i++) { + delete integrators[i]; + } +} + +void CompositeNonlinearIntegrator::add_integrator(mfem::NonlinearFormIntegrator *integrator) { + integrators.push_back(integrator); +} + +void CompositeNonlinearIntegrator::AssembleElementVector( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::Vector &elvect) { + int dof = el.GetDof(); + elvect.SetSize(dof); + elvect = 0.0; + mfem::Vector temp(dof); + + for (size_t i = 0; i < integrators.size(); i++) { + temp= 0.0; + integrators[i]->AssembleElementVector(el, Trans, elfun, temp); + elvect.Add(1.0, temp); + } +} + +void CompositeNonlinearIntegrator::AssembleElementGrad( + const mfem::FiniteElement &el, + mfem::ElementTransformation &Trans, + const mfem::Vector &elfun, + mfem::DenseMatrix &elmat) { + int dof = el.GetDof(); + elmat.SetSize(dof, dof); + elmat = 0.0; + mfem::DenseMatrix temp(dof); + temp.SetSize(dof, dof); + for (size_t i = 0; i < integrators.size(); i++) { + temp = 0.0; + integrators[i] -> AssembleElementGrad(el, Trans, elfun, temp); + elmat.Add(1.0, temp); + } +} \ No newline at end of file diff --git a/src/poly/utils/public/polyIO.h b/src/poly/utils/public/polyIO.h new file mode 100644 index 0000000..9ca7ceb --- /dev/null +++ b/src/poly/utils/public/polyIO.h @@ -0,0 +1,16 @@ +#ifndef POLY_IO_H +#define POLY_IO_H + +#include "mfem.hpp" +#include + +/** + * @brief Writes the solution to a CSV file. + * + * @param u The GridFunction containing the solution. + * @param mesh The mesh associated with the solution. + * @param filename The name of the CSV file to write to. + */ +void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); + +#endif // POLY_IO_H \ No newline at end of file diff --git a/src/poly/utils/public/polyMFEMUtils.h b/src/poly/utils/public/polyMFEMUtils.h new file mode 100644 index 0000000..ac07264 --- /dev/null +++ b/src/poly/utils/public/polyMFEMUtils.h @@ -0,0 +1,128 @@ +#include "mfem.hpp" +#include + + +void write_solution_to_csv(const mfem::GridFunction &u, const mfem::Mesh &mesh, const std::string &filename); + +/** + * @brief A class for nonlinear power integrator. + */ +class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator { +private: + mfem::FunctionCoefficient coeff_; + double polytropicIndex; +public: + /** + * @brief Constructor for NonlinearPowerIntegrator. + * + * @param coeff The function coefficient. + * @param n The polytropic index. + */ + NonlinearPowerIntegrator(mfem::FunctionCoefficient &coeff, double n); + + /** + * @brief Assembles the element vector. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elvect The element vector to be assembled. + */ + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + + /** + * @brief Assembles the element gradient. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elmat The element matrix to be assembled. + */ + virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; + +/** + * @brief A wrapper class for bilinear integrator. + */ +class BilinearIntegratorWrapper : public mfem::NonlinearFormIntegrator +{ +private: + mfem::BilinearFormIntegrator *integrator; +public: + /** + * @brief Constructor for BilinearIntegratorWrapper. + * + * @param integratorInput The bilinear form integrator input. + */ + BilinearIntegratorWrapper(mfem::BilinearFormIntegrator *integratorInput); + + /** + * @brief Destructor for BilinearIntegratorWrapper. + */ + virtual ~BilinearIntegratorWrapper(); + + /** + * @brief Assembles the element vector. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elvect The element vector to be assembled. + */ + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + + /** + * @brief Assembles the element gradient. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elmat The element matrix to be assembled. + */ + virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; + +/** + * @brief A class for composite nonlinear integrator. + */ +class CompositeNonlinearIntegrator: public mfem::NonlinearFormIntegrator { + private: + std::vector integrators; + public: + /** + * @brief Constructor for CompositeNonlinearIntegrator. + */ + CompositeNonlinearIntegrator(); + + /** + * @brief Destructor for CompositeNonlinearIntegrator. + */ + virtual ~CompositeNonlinearIntegrator(); + + /** + * @brief Adds an integrator to the composite integrator. + * + * @param integrator The nonlinear form integrator to add. + */ + void add_integrator(mfem::NonlinearFormIntegrator *integrator); + + /** + * @brief Assembles the element vector. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elvect The element vector to be assembled. + */ + virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override; + + /** + * @brief Assembles the element gradient. + * + * @param el The finite element. + * @param Trans The element transformation. + * @param elfun The element function. + * @param elmat The element matrix to be assembled. + */ + virtual void AssembleElementGrad(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override; +}; \ No newline at end of file diff --git a/src/resources/mesh/sphere.msh b/src/resources/mesh/sphere.msh new file mode 100644 index 0000000..416c655 Binary files /dev/null and b/src/resources/mesh/sphere.msh differ diff --git a/subprojects/mfem.wrap b/subprojects/mfem.wrap new file mode 100644 index 0000000..43c167b --- /dev/null +++ b/subprojects/mfem.wrap @@ -0,0 +1,5 @@ +[wrap-git] +url = https://github.com/mfem/mfem.git +revision = master + +[cmake] \ No newline at end of file diff --git a/tests/meshIO/meshIOTest.cpp b/tests/meshIO/meshIOTest.cpp new file mode 100644 index 0000000..07eec4f --- /dev/null +++ b/tests/meshIO/meshIOTest.cpp @@ -0,0 +1,38 @@ +#include +#include "meshIO.h" +#include +#include +#include "mfem.hpp" + +std::string EXAMPLE_FILENAME = std::string(getenv("MESON_SOURCE_ROOT")) + "/src/resources/mesh/sphere.msh"; + +double ComputeMeshVolume(mfem::Mesh &mesh) +{ + double totalVolume = 0.0; + for (int i = 0; i < mesh.GetNE(); i++) // Loop over all elements + { + totalVolume += mesh.GetElementVolume(i); + } + return totalVolume; +} + + +class meshIOTest : public ::testing::Test {}; + +TEST_F(meshIOTest, DefaultConstructor) { + EXPECT_NO_THROW(MeshIO meshIO(EXAMPLE_FILENAME)); +} + +TEST_F(meshIOTest, IsLoaded) { + MeshIO meshIO(EXAMPLE_FILENAME); + EXPECT_EQ(meshIO.IsLoaded(), true); +} + +TEST_F(meshIOTest, GetMesh) { + MeshIO meshIO(EXAMPLE_FILENAME); + mfem::Mesh& mesh = meshIO.GetMesh(); + EXPECT_EQ(mesh.GetNE(), 670); + EXPECT_EQ(mesh.GetNV(), 201); + double volume = ComputeMeshVolume(mesh); + EXPECT_DOUBLE_EQ(volume, 3.9357596288315868); +} \ No newline at end of file diff --git a/tests/meshIO/meson.build b/tests/meshIO/meson.build new file mode 100644 index 0000000..6db5297 --- /dev/null +++ b/tests/meshIO/meson.build @@ -0,0 +1,22 @@ +# Test files for mesh +test_sources = [ + 'meshIOTest.cpp', +] + +foreach test_file : test_sources + exe_name = test_file.split('.')[0] + message('Building test: ' + exe_name) + + # Create an executable target for each test + test_exe = executable( + exe_name, + test_file, + dependencies: [gtest_dep, mfem_dep], + include_directories: include_directories('../../src/meshIO/public'), + link_with: libmeshIO, # Link the dobj library + install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly + ) + + # Add the executable as a test + test(exe_name, test_exe) +endforeach diff --git a/tests/meson.build b/tests/meson.build index 45693ff..162daab 100644 --- a/tests/meson.build +++ b/tests/meson.build @@ -6,6 +6,7 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true) subdir('dobj') subdir('const') subdir('opatIO') +subdir('meshIO') # Subdirectories for sandbox tests subdir('dobj_sandbox') diff --git a/utils/meshGeneration/generateMesh.py b/utils/meshGeneration/generateMesh.py new file mode 100644 index 0000000..eb76bb2 --- /dev/null +++ b/utils/meshGeneration/generateMesh.py @@ -0,0 +1,24 @@ +import pygmsh +import meshio + +import argparse + +def generate_spherical_mesh(radius=1, meshSize=0.1): + with pygmsh.geo.Geometry() as geo: + # Create a spherical (ball) geometry centered at (0,0,0) + sphere = geo.add_ball([0, 0, 0], radius) + # Generate a 2D mesh (i.e. surface mesh) of the sphere + myMesh = geo.generate_mesh(dim=3) + return myMesh + +def write_mfem_mesh(meshData, filename): + meshio.write(filename, meshData, file_format='gmsh22') + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Generate a spherical mesh') + parser.add_argument('--radius', type=float, default=1.0, help='Radius of the sphere') + parser.add_argument('--meshSize', type=float, default=0.1, help='Mesh size') + args = parser.parse_args() + + meshData = generate_spherical_mesh(args.radius, args.meshSize) + write_mfem_mesh(meshData, 'sphere.msh') \ No newline at end of file diff --git a/utils/meshGeneration/readme.md b/utils/meshGeneration/readme.md new file mode 100644 index 0000000..27c4d15 --- /dev/null +++ b/utils/meshGeneration/readme.md @@ -0,0 +1,8 @@ +# spherical mesh generation +A simple script to generate the base spherical mesh which 4DSSE uses to solve equations. +This will produce a unit sphere mesh centered on (0, 0, 0) which can be rescaled within the +4DSSE program. + +## Dependecies +- pygmsh +- meshio \ No newline at end of file diff --git a/utils/meshGeneration/sphere.msh b/utils/meshGeneration/sphere.msh new file mode 100644 index 0000000..416c655 Binary files /dev/null and b/utils/meshGeneration/sphere.msh differ diff --git a/utils/opatio/src/opatio/opat/opat.py b/utils/opatio/src/opatio/opat/opat.py index 741a761..eea33da 100644 --- a/utils/opatio/src/opatio/opat/opat.py +++ b/utils/opatio/src/opatio/opat/opat.py @@ -61,6 +61,46 @@ defaultHeader = Header( ) class OpatIO: + """ + @brief Class for handling OPAT file input/output operations. + This class provides methods to validate, manipulate, and save OPAT files. It includes functionalities to validate character arrays, 1D arrays, and 2D arrays, compute checksums, set header information, add tables, and save the OPAT file in both ASCII and binary formats. + Attributes: + header (Header): The header of the OPAT file. + tables (List[Tuple[Tuple[float, float], OPATTable]]): A list of tables in the OPAT file. + Methods: + validate_char_array_size(s: str, nmax: int) -> bool: + Validate the size of a character array. + validate_logKappa(logKappa): + Validate the logKappa array. + validate_1D(arr, name: str): + Validate a 1D array. + compute_checksum(data: bytes) -> bytes: + Compute the SHA-256 checksum of the given data. + set_version(version: int) -> int: + Set the version of the OPAT file. + set_source(source: str) -> str: + Set the source information of the OPAT file. + set_comment(comment: str) -> str: + Set the comment of the OPAT file. + add_table(X: float, Z: float, logR: Iterable[float], logT: Iterable[float], logKappa: Iterable[Iterable[float]]): + Add a table to the OPAT file. + _header_bytes() -> bytes: + Convert the header to bytes. + _table_bytes(table: OPATTable) -> Tuple[bytes, bytes]: + Convert a table to bytes. + _tableIndex_bytes(tableIndex: TableIndex) -> bytes: + Convert a table index to bytes. + __repr__() -> str: + Get the string representation of the OpatIO object. + _format_table_as_string(table: OPATTable, X: float, Z: float) -> str: + Format a table as a string. + print_table_indexes(table_indexes: List[TableIndex]) -> str: + Print the table indexes. + save_as_ascii(filename: str) -> str: + Save the OPAT file as an ASCII file. + save(filename: str) -> str: + Save the OPAT file as a binary file. + """ def __init__(self): self.header: Header = defaultHeader self.tables: List[Tuple[Tuple[float, float], OPATTable]] = []