diff --git a/CMakeLists.txt b/CMakeLists.txt index 3a8e022..8208d07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,431 +1,432 @@ #################################################### # # # SpecMiCP : CMakeLists.txt # # # #################################################### project(specmicp) cmake_minimum_required(VERSION 2.8.8) # For an explanation of the options see the INSTALL file # External Package ######################################################################### set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/") # Boost # ----- find_package(Boost REQUIRED) include_directories(${Boost_INCLUDE_DIR}) # Eigen # ----- find_package(Eigen3 REQUIRED) # This module comes from the Eigen3 Package include_directories(${EIGEN3_INCLUDE_DIR}) # Eigen unsuported # GMRES.h is really the file we are using/looking for # If it doesn't exist then the solver will not be included in the list of the parse solvers if(EXISTS "${EIGEN3_INCLUDE_DIR}/unsupported/Eigen/src/IterativeSolvers/GMRES.h") add_definitions(-DEIGEN_UNSUPPORTED_FOUND) INCLUDE_DIRECTORIES("${EIGEN3_INCLUDE_DIR}/unsupported/") endif() # Compilation flags ######################################################################### if(UNIX) # Generic options SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native -std=c++11 -fuse-ld=gold") SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -pedantic") SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O2 -DNDEBUG") # OpenMP if(SPECMICP_USE_OPENMP) SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp") add_definitions(-DSPECMICP_USE_OPENMP) endif() message(STATUS "c++ flags ${CMAKE_CXX_FLAGS}") else() message(WARNING "not tested !") endif() # Activate this option to use finite difference jacobian # add_definitions(-DSPECMICP_DEBUG_EQUATION_FD_JACOBIAN) # Directories ######################################################################### # 3rd party # ========= # JsonCPP # ------- set(JSONCPP_DIR ${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) add_custom_target(3rdparty_header SOURCES ${JSONCPP_DIR}/json/json.h ${JSONCPP_DIR}/json/json-forwards.h ) include_directories(${PROJECT_SOURCE_DIR}/3rdparty/jsoncpp/) # Source # ======= set(PROJECT_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src) # Includes # ========= include_directories(${PROJECT_SOURCE_DIR}) add_custom_target(common SOURCES ${PROJECT_SOURCE_DIR}/common.hpp ${PROJECT_SOURCE_DIR}/database.hpp ) # Modules ######################################################################### # MiCPSolver # ============== # Mixed complementarity Solver # include only module set(MICPSOLVER_DIR ${PROJECT_SOURCE_DIR}/micpsolver) add_custom_target(micpsolver_inc SOURCES ${MICPSOLVER_DIR}/micpsolver_base.hpp ${MICPSOLVER_DIR}/micpsolver_base.inl ${MICPSOLVER_DIR}/micpsolver.hpp ${MICPSOLVER_DIR}/micpsolver.inl ${MICPSOLVER_DIR}/micpsolver_min.hpp ${MICPSOLVER_DIR}/micpsolver_min.inl ${MICPSOLVER_DIR}/micpsolverold.hpp ${MICPSOLVER_DIR}/micpsolverold.inl ${MICPSOLVER_DIR}/micpsolver_structs.hpp ${MICPSOLVER_DIR}/micpprog.hpp ${MICPSOLVER_DIR}/ncp_function.hpp ${MICPSOLVER_DIR}/estimate_cond_number.hpp ) # SpecMiCP # ============ # The speciation solver set(SPECMICP_DIR ${PROJECT_SOURCE_DIR}/specmicp) set(KINETICS_DIR ${SPECMICP_DIR}/kinetics) set(ADIMENSIONAL_DIR ${SPECMICP_DIR}/adimensional) set(PROBLEM_SOLVER_DIR ${SPECMICP_DIR}/problem_solver) # headers only module add_custom_target(specmic_inc SOURCES ${SPECMICP_DIR}/simulation_box.hpp ${SPECMICP_DIR}/reduced_system_structs.hpp ${SPECMICP_DIR}/reduced_system_solver_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_structs.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solution.hpp ${ADIMENSIONAL_DIR}/adimensional_system_solver_structs.hpp ${PROBLEM_SOLVER_DIR}/formulation.hpp ${KINETICS_DIR}/kinetic_model.hpp ${KINETICS_DIR}/kinetic_variables.hpp ${KINETICS_DIR}/kinetic_system_solver_structs.hpp ) # Source for the library set(SPECMICP_LIBFILE ${SPECMICP_DIR}/equilibrium_data.cpp ${SPECMICP_DIR}/reduced_system.cpp ${SPECMICP_DIR}/reduced_system_solver.cpp ${SPECMICP_DIR}/reaction_path.cpp ${ADIMENSIONAL_DIR}/adimensional_system.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solver.cpp ${ADIMENSIONAL_DIR}/adimensional_system_solution_extractor.cpp ${ADIMENSIONAL_DIR}/equilibrium_curve.cpp ${PROBLEM_SOLVER_DIR}/dissolver.cpp # kinetic ${KINETICS_DIR}/kinetic_system.cpp ${KINETICS_DIR}/kinetic_system_solver.cpp ) add_library(specmicp ${SPECMICP_LIBFILE}) # Database # ============ set(DATABASE_DIR ${PROJECT_SOURCE_DIR}/database) set(DATABASE_LIBFILE ${DATABASE_DIR}/data_container.cpp ${DATABASE_DIR}/database.cpp ${DATABASE_DIR}/reader.cpp ${DATABASE_DIR}/selector.cpp ${DATABASE_DIR}/mineral_selector.cpp ${DATABASE_DIR}/canonicalizer.cpp ${DATABASE_DIR}/switch_basis.cpp ${JSONCPP_DIR}/jsoncpp.cpp ) add_custom_target(data_incl SOURCES ${DATABASE_DIR}/common_def.hpp ${DATABASE_DIR}/module.hpp ) add_library(specmicp_database ${DATABASE_LIBFILE}) # ODEint # ========== # ODE integration, headers-only set(ODEINT_DIR ${PROJECT_SOURCE_DIR}/odeint) add_custom_target(odeint_incl SOURCES ${ODEINT_DIR}/odeint_types.hpp ${ODEINT_DIR}/butcher_tableau.hpp ${ODEINT_DIR}/runge_kutta_step.hpp ${ODEINT_DIR}/runge_kutta_step.inl ${ODEINT_DIR}/runge_kutta_step_structs.hpp ) # DFPMSolver # ============== ## Parabolic solver for finite element problem, header only set(DFPMSOLVER_DIR ${PROJECT_SOURCE_DIR}/dfpmsolver) add_custom_target(dfpmsolver_incl SOURCES ${DFPMSOLVER_DIR}/driver.hpp ${DFPMSOLVER_DIR}/driver.inl ${DFPMSOLVER_DIR}/driver_structs.hpp ${DFPMSOLVER_DIR}/dfpm_program.hpp # # Parabolic solver ${DFPMSOLVER_DIR}/parabolic_driver.hpp ${DFPMSOLVER_DIR}/parabolic_driver.inl ${DFPMSOLVER_DIR}/parabolic_program.hpp ${DFPMSOLVER_DIR}/parabolic_structs.hpp ) # DFPM # ======== # Dynaflow+- , Finite element solver, mostly headers set (DFPM_DIR ${PROJECT_SOURCE_DIR}/dfpm) add_custom_target(dfpm_incl SOURCES ${DFPM_DIR}/types.hpp ${DFPM_DIR}/1dtransport/diffusion_parameters.hpp ${DFPM_DIR}/mesh.hpp ${DFPM_DIR}/meshes/mesh1dfwd.hpp ${DFPM_DIR}/meshes/mesh1d.hpp ${DFPM_DIR}/meshes/uniform_mesh1d.hpp ${DFPM_DIR}/meshes/axisymmetric_uniform_mesh1d.hpp + ${DFPM_DIR}/meshes/axisymmetric_mesh1d.hpp ) set(DFPMLIB ${DFPM_DIR}/1dtransport/diffusion.cpp ) add_library(dfpm ${DFPMLIB}) # ReactMiCP # ============= # Reactive transport solver set(REACTMICP_DIR ${PROJECT_SOURCE_DIR}/reactmicp) set(REACTMICP_SOLVER_DIR ${REACTMICP_DIR}/solver/) set(REACTMICP_SYSTEMS_DIR ${REACTMICP_DIR}/systems/) add_custom_target(reactmicp_incl SOURCES ${REACTMICP_DIR}/common.hpp # ${REACTMICP_DIR}/micpsolvers/parabolicprogram.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver_structs.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver.hpp ## remove or fix ${REACTMICP_DIR}/micpsolvers/micpsolver.inl ## remove or fix # ${REACTMICP_DIR}/systems/secondary_variables/secondary.hpp ${REACTMICP_DIR}/systems/secondary_variables/secondary.inl ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver_structs.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_parameters_base.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_parameters.hpp ${REACTMICP_DIR}/systems/saturated_diffusion/boundary_conditions.hpp # Flexible reactive tranport solver # --------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/variables_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/transport_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/chemistry_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/upscaling_stagger_base.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/stagger_structs.hpp ${REACTMICP_SOLVER_DIR}/staggers_base/decl.inl ${REACTMICP_SOLVER_DIR}/staggers_base/staggers_base.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/variablesfwd.hpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/input_data.hpp ) set(REACTMICP_LIBFILES ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/reactive_transport_neutrality_solver.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_program.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/transport_neutrality_variables.cpp ${REACTMICP_DIR}/systems/saturated_diffusion/variables.cpp ${REACTMICP_DIR}/equilibrium_curve/chemistry.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_extractor.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_coupler.cpp ${REACTMICP_DIR}/equilibrium_curve/eqcurve_solid_transport.cpp # Flexible reactive transport solver # ---------------------------------- ${REACTMICP_SOLVER_DIR}/reactive_transport_solver.cpp # ${REACTMICP_SYSTEMS_DIR}/saturated_react/variables.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_program.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/transport_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/equilibrium_stagger.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/react_solver.cpp ${REACTMICP_SYSTEMS_DIR}/saturated_react/init_variables.cpp ) add_library(reactmicp ${REACTMICP_LIBFILES}) # Utils # ========= # Stuff that may be used by the other modules set(UTILS_DIR ${PROJECT_SOURCE_DIR}/utils) add_custom_target(utils_inc SOURCES ${UTILS_DIR}/log.hpp # sparse solvers # -------------- ${UTILS_DIR}/sparse_solvers/sparse_solver.hpp ${UTILS_DIR}/sparse_solvers/sparse_solver_base.hpp ${UTILS_DIR}/sparse_solvers/sparse_solver_structs.hpp ${UTILS_DIR}/sparse_solvers/sparse_qr.hpp ${UTILS_DIR}/sparse_solvers/sparse_lu.hpp ${UTILS_DIR}/sparse_solvers/sparse_bicgstab.hpp ${UTILS_DIR}/sparse_solvers/sparse_gmres.hpp ${UTILS_DIR}/options_handler.hpp ${UTILS_DIR}/perfs_handler.hpp ${UTILS_DIR}/moving_average.hpp ) # Physics # =========== # Stuff that may be used by the other modules, but is related to physics set(PHYSICS_DIR ${PROJECT_SOURCE_DIR}/physics) add_custom_target(physics_inc SOURCES ${PHYSICS_DIR}/constants.hpp ${PHYSICS_DIR}/laws.hpp ${PHYSICS_DIR}/units.hpp ) # Database ######################################################################### # Not code, just necessary data set(DATA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/data) set(DATABASE_LIST ${DATA_DIR}/specmicp_database.js ${DATA_DIR}/cemdata_specmicp.js ) add_custom_target(data SOURCES ${DATABASE_LIST}) file(INSTALL ${DATABASE_LIST} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/data/) # Documentation ######################################################################### # "common" documentation # ----------------------- set( DOCS_LIST ${CMAKE_CURRENT_SOURCE_DIR}/README.md ${CMAKE_CURRENT_SOURCE_DIR}/INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/TODO ${CMAKE_CURRENT_SOURCE_DIR}/COPYING ) add_custom_target(docs SOURCES ${DOCS_LIST}) # Doxygen documentation # --------------------- # add a target to generate API documentation with Doxygen find_package(Doxygen) if(DOXYGEN_FOUND) # Configure doxygen configure_file(${CMAKE_CURRENT_SOURCE_DIR}/doc/Doxyfile.in ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile @ONLY) # Citations for the documentations add_custom_target(citations SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib) file(INSTALL ${CMAKE_CURRENT_SOURCE_DIR}/doc/citations.bib DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/doc/) # Target to build the documentations add_custom_target(doc ${DOXYGEN_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/doc/Doxyfile WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMENT "Generating API documentation with Doxygen" VERBATIM ) endif(DOXYGEN_FOUND) # The following modules have their own CMakeLists.txt # Python interface ######################################################################### add_subdirectory( cython ) # Tests ######################################################################### add_subdirectory( tests ) # Examples ######################################################################### add_subdirectory( examples ) diff --git a/src/dfpm/meshes/axisymmetric_mesh1d.hpp b/src/dfpm/meshes/axisymmetric_mesh1d.hpp new file mode 100644 index 0000000..baebe08 --- /dev/null +++ b/src/dfpm/meshes/axisymmetric_mesh1d.hpp @@ -0,0 +1,125 @@ +#ifndef SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP +#define SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP + +#include "mesh1d.hpp" + +namespace specmicp { +namespace mesh { + + +//! \brief A generic Axisymmetric 1D mesh +class AxisymmetricMesh1D: public Mesh1D +{ +public: + const index_t radius = 0; + const index_t cell_vol_0 = 1; + const index_t cell_vol_1 = 2; + + //! \brief Build an axisymmetric 1D mesh + //! + //! \param radius is a vector of the position of the nodes + //! \param height is the height of the sample (or depth) + AxisymmetricMesh1D(Vector radius, scalar_t height): + Mesh1D(radius.rows()), + m_height(height), + m_data(radius.rows(), 3), + m_face_radius(radius.rows()-1) + { + // radius + m_data.col(0) = radius; + // Radius of the faces + for (index_t element=0; element( + std::make_shared(radius, height)); + +} + +//! \brief Factory method to build a pointer to a uniform axisymmetric 1D mesh +//! +//! \param nb_nodes number of nodes in the mesh +//! \param radius of the first node (biggest radius in the mesh) +//! \param dx length of an element +//! \param height is the height of the sample (or depth) +//! +//! Note: this method buil an AxisymmetricMesh1D instance while +//! the specmicp::mesh::axisymmetric_uniform_mesh1d method build a +//! specmicp::mesh::AxisymmetricUniformMesh1D instance. +inline Mesh1DPtr uniform_axisymmetric_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height) +{ + Vector radius_s(nb_nodes); + for (index_t node=0; node( + std::make_shared(radius_s, height)); + +} + +} // end namespace mesh +} // end namespace specmicp + +#endif // SPECMICP_DFPM_MESHES_AXISYMMETRICMESH1D_HPP diff --git a/src/dfpm/meshes/axisymmetric_uniform_mesh1d.hpp b/src/dfpm/meshes/axisymmetric_uniform_mesh1d.hpp index 086928d..46bb873 100644 --- a/src/dfpm/meshes/axisymmetric_uniform_mesh1d.hpp +++ b/src/dfpm/meshes/axisymmetric_uniform_mesh1d.hpp @@ -1,104 +1,104 @@ /*------------------------------------------------------- - Module : reactmicp/meshes - File : axisymmetric_uniform_mesh1d.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP #define SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP #include "mesh1d.hpp" // \file axisymmetric_uniform_mesh1d.hpp Uniform axisymmetric 1D mesh namespace specmicp { namespace mesh { //! \brief A uniform 1D mesh class AxisymmetricUniformMesh1D: public Mesh1D { public: AxisymmetricUniformMesh1D(index_t nb_nodes, scalar_t radius, scalar_t height): Mesh1D(nb_nodes), m_radius(radius), m_radius_int(0.0), m_dx(radius/(nb_nodes-1)), m_height(height) {} AxisymmetricUniformMesh1D(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height): Mesh1D(nb_nodes), m_radius(radius), m_radius_int(radius-(nb_nodes-1)*dx), m_dx(dx), m_height(height) {} scalar_t get_radius_node(index_t node) { return m_radius-m_dx*node; } //! Return the radius of a face of an element scalar_t get_radius_face(index_t element) { return get_radius_node(get_node(element, 1))+m_dx/2; } scalar_t get_dx(index_t _) {return m_dx;} scalar_t get_face_area(index_t element) { return 2*M_PI*get_radius_face(element)*m_height;} scalar_t get_volume_element(index_t element) {return M_PI*m_height*( std::pow(get_radius_node(get_node(element, 0)),2) - std::pow(get_radius_node(get_node(element, 1)),2)) ;} scalar_t get_volume_cell(index_t node) { if (node ==0) return M_PI*m_height*(std::pow(m_radius,2)-std::pow(get_radius_face(0),2)); else if (node == nb_nodes()-1) if (m_radius_int == 0.0) return M_PI*m_height*std::pow(m_dx/2.0,2); else return M_PI*m_height*(std::pow(get_radius_face(node-1),2) - std::pow(m_radius_int,2)); else return M_PI*m_height*(std::pow(get_radius_face(node-1),2) - std::pow(get_radius_face(node),2)); } scalar_t get_volume_cell_element(index_t element, index_t enode) { if (enode == 0) { return M_PI*m_height*(std::pow(get_radius_node(get_node(element,enode)),2) - std::pow(get_radius_face(element),2)); } else { return M_PI*m_height*(std::pow(get_radius_face(element),2) - std::pow(get_radius_node(get_node(element,enode)),2)); } } private: scalar_t m_radius; scalar_t m_radius_int; scalar_t m_dx; scalar_t m_height; }; //! \brief Factory method to build a pointer to a uniform 1D mesh -inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t cross_section) +inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t height) { return std::static_pointer_cast( - std::make_shared(nb_nodes, radius, cross_section)); + std::make_shared(nb_nodes, radius, height)); } //! \brief Factory method to build a pointer to a uniform 1D mesh -inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t cross_section) +inline Mesh1DPtr axisymmetric_uniform_mesh1d(index_t nb_nodes, scalar_t radius, scalar_t dx, scalar_t height) { return std::static_pointer_cast( - std::make_shared(nb_nodes, radius, dx, cross_section)); + std::make_shared(nb_nodes, radius, dx, height)); } } // end namespace mesh } // end namespace specmicp #endif // SPECMICP_REACTMICP_MESHES_AXISYMMETRICMESH1D_HPP diff --git a/tests/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp b/tests/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp index 0eb250a..91f858b 100644 --- a/tests/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp +++ b/tests/reactmicp/meshes/test_axisymmetric_uniform_mesh_1d.cpp @@ -1,42 +1,86 @@ #include "catch.hpp" #include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp" +#include "dfpm/meshes/axisymmetric_mesh1d.hpp" + using namespace specmicp::mesh; TEST_CASE("Axisymmetric uniform mesh 1d", "[Reactmicp, mesh, axisymmetric, uniform, 1D]") { SECTION("test uniform mesh 1d") { double nb_node = 6; double radius = 5; double dx = 1; double height = 2; AxisymmetricUniformMesh1D mesh(nb_node, radius, height); REQUIRE(mesh.nb_nodes() == nb_node); REQUIRE(mesh.nb_elements() == nb_node - 1); REQUIRE(mesh.get_dx(0) == dx); REQUIRE(mesh.get_radius_face(0) == 4.5); REQUIRE(mesh.get_radius_face(1) == 3.5); REQUIRE(std::abs(mesh.get_face_area(0) - 2*M_PI*4.5*height) < 1e-10); REQUIRE(std::abs(mesh.get_volume_element(0) - M_PI*height*(5*5-4*4)) < 1e-10); REQUIRE(std::abs(mesh.get_volume_cell(0) - M_PI*height*(5*5-4.5*4.5)) < 1e-10); REQUIRE(std::abs(mesh.get_volume_element(1) - M_PI*height*(4*4-3*3)) < 1e-10); REQUIRE(std::abs(mesh.get_volume_cell(1) - M_PI*height*(4.5*4.5-3.5*3.5)) < 1e-10); int check_nb_node = 0; for (int node: mesh.range_nodes()) {++check_nb_node;} REQUIRE(check_nb_node == nb_node); int check_nb_element = 0; for (int element: mesh.range_elements()) {++check_nb_element;} REQUIRE(check_nb_element == mesh.nb_elements()); int check_nen = 0; for (int node: mesh.range_nen()) {++check_nen;} REQUIRE(check_nen == 2); } } + + +TEST_CASE("Axisymmetric mesh 1d", "[Reactmicp, mesh, axisymmetric, 1D]") { + + SECTION("test uniform mesh 1d") { + + double nb_node = 6; + double radius = 5; + double dx = 1; + double height = 2; + + Mesh1DPtr uniform_mesh = axisymmetric_uniform_mesh1d(nb_node, radius, dx, height); + Mesh1DPtr generic_mesh = uniform_axisymmetric_mesh1d(nb_node, radius, dx, height); + + CHECK(generic_mesh->nb_nodes() == nb_node); + CHECK(uniform_mesh->nb_nodes() == generic_mesh->nb_nodes()); + CHECK(uniform_mesh->nb_elements() == generic_mesh->nb_elements()); + + CHECK(generic_mesh->get_dx(0) == dx); + + CHECK(generic_mesh->get_face_area(0) == Approx(2*M_PI*4.5*height)); + CHECK(generic_mesh->get_face_area(0) == uniform_mesh->get_face_area(0)); + + CHECK(generic_mesh->get_volume_element(0) == Approx(M_PI*height*(5*5-4*4))); + CHECK(generic_mesh->get_volume_cell(0) == Approx(M_PI*height*(5*5-4.5*4.5))); + CHECK(generic_mesh->get_volume_element(1) == Approx(M_PI*height*(4*4-3*3))); + CHECK(generic_mesh->get_volume_cell(1) == Approx(M_PI*height*(4.5*4.5-3.5*3.5))); + + int check_nb_node = 0; + for (int node: generic_mesh->range_nodes()) {++check_nb_node;} + CHECK(check_nb_node == nb_node); + + int check_nb_element = 0; + for (int element: generic_mesh->range_elements()) {++check_nb_element;} + CHECK(check_nb_element == generic_mesh->nb_elements()); + + int check_nen = 0; + for (int node: generic_mesh->range_nen()) {++check_nen;} + CHECK(check_nen == 2); + } +} +