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/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..74c7770 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 +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/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/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