diff --git a/.gitignore b/.gitignore index 4440c8b..8baa052 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ /build/CMake* *~ /build/* /build-o3/ -/build* \ No newline at end of file +/build* + +.vscode/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 598fdb5..5c0db67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,263 +1,263 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief Main configuration file # # @section LICENSE # # Copyright © 2018 Till Junge # # µSpectre is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License as # published by the Free Software Foundation, either version 3, or (at # your option) any later version. # # µSpectre is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with GNU Emacs; see the file COPYING. If not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # ============================================================================= cmake_minimum_required(VERSION 3.0.0) project(µSpectre) set(CMAKE_CXX_STANDARD 14) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) set(BUILD_SHARED_LIBS ON) set(MUSPECTRE_PYTHON_MAJOR_VERSION 3) add_compile_options(-Wall -Wextra -Weffc++) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) set(MAKE_DOC_TARGET "OFF" CACHE BOOL "If on, a target dev_doc (which builds the documentation) is added") set(MAKE_TESTS "ON" CACHE BOOL "If on, several ctest targets will be built automatically") set(MAKE_EXAMPLES "ON" CACHE BOOL "If on, the executables in the bin folder will be compiled") set(MAKE_BENCHMARKS "ON" CACHE BOOL "If on, the benchmarks will be compiled") set(MPI_PARALLEL "OFF" CACHE BOOL "If on, MPI-parallel solvers become available") set(NETCDF "OFF" CACHE BOOL "If on, netCDF support becomes available") -set(MULIB "OFF" CACHE BOOL "If on, mulib input reader is build. turns on NETCDF") +set(MULIB "ON" CACHE BOOL "If on, mulib input reader is build. turns on NETCDF") set(RUNNING_IN_CI "OFF" CACHE INTERNAL "changes output format for tests") if(${MULIB}) set(NETCDF "ON" CACHE BOOL "must be on now" FORCE) add_definitions(-DWITH_MULIB) endif(${MULIB}) if(${MAKE_TESTS}) enable_testing() find_package(Boost COMPONENTS unit_test_framework REQUIRED) endif(${MAKE_TESTS}) if(${MPI_PARALLEL}) add_definitions(-DWITH_MPI) find_package(MPI) if (NOT ${MPI_FOUND}) message(SEND_ERROR "You chose MPI but CMake cannot find the MPI package") endif(NOT ${MPI_FOUND}) endif(${MPI_PARALLEL}) include(CheckCXXSourceCompiles) if(${NETCDF}) set(NETCDF_CXX "YES") find_package(NetCDF REQUIRED) if (${NETCDF}) # the following checks whether std::filesystem exists and replaces # it by boost::filesystem if necessary set (CMAKE_REQUIRED_LIBRARIES stdc++fs) check_cxx_source_compiles( "#include int main() { std::experimental::filesystem::path(\"/bin/ban/bun\"); }" HAS_STD_FILESYSTEM ) endif (${NETCDF}) if (NOT HAS_STD_FILESYSTEM) add_definitions (-DNO_FILESYSTEM) endif (NOT HAS_STD_FILESYSTEM) endif(${NETCDF}) include(muspectreTools) string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") # using Clang add_compile_options(-Wno-missing-braces) if ("debug" STREQUAL "${build_type}") add_compile_options(-O0) endif() elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # using GCC add_compile_options(-Wno-non-virtual-dtor) add_compile_options(-march=native) if (("relwithdebinfo" STREQUAL "${build_type}") OR ("release" STREQUAL "${build_type}" )) add_compile_options(-march=native) endif() if ("debug" STREQUAL "${build_type}" ) add_compile_options(-O0) endif() elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") # using Intel C++ elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC") # using Visual Studio C++ endif() # Do not trust old gcc. the std::optional has memory bugs if(${CMAKE_COMPILER_IS_GNUCC}) if(${CMAKE_CXX_COMPILER_VERSION} VERSION_LESS 6.0.0) add_definitions(-DNO_EXPERIMENTAL) endif() endif() add_external_package(Eigen3 VERSION 3.3.0 CONFIG) add_external_package(pybind11 VERSION 2.2 CONFIG) include_directories( ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR} ) if(APPLE) include_directories(${CMAKE_INSTALL_PREFIX}/include ${Boost_INCLUDE_DIRS}) endif() #build tests (these are before we add -Werror to the compile options) if (${MAKE_TESTS}) ############################################################################## # build library tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc") if (${NETCDF}) file (GLOB TEST_SRCS ${TEST_SRCS} "${CMAKE_SOURCE_DIR}/tests/netcdf_test_*.cc") endif (${NETCDF}) add_executable(main_test_suite tests/main_test_suite.cc ${TEST_SRCS}) target_link_libraries(main_test_suite ${Boost_LIBRARIES} muSpectre) muSpectre_add_test(main_test_suite TYPE BOOST main_test_suite --report_level=detailed) # build header tests file( GLOB HEADER_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/header_test_*.cc") foreach(header_test ${HEADER_TEST_SRCS}) get_filename_component(header_test_name ${header_test} NAME_WE) string(SUBSTRING ${header_test_name} 12 -1 test_name) list(APPEND header_tests ${test_name}) add_executable(${test_name} tests/main_test_suite.cc ${header_test}) target_link_libraries(${test_name} ${Boost_LIBRARIES} Eigen3::Eigen) target_include_directories(${test_name} INTERFACE ${muSpectre_INCLUDES}) muSpectre_add_test(${test_name} TYPE BOOST ${test_name} --report_level=detailed) endforeach(header_test ${HEADER_TEST_SRCS}) add_custom_target(header_tests) add_dependencies(header_tests ${header_tests}) ############################################################################## # copy python test file( GLOB PY_TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/python_*.py") foreach(pytest ${PY_TEST_SRCS}) get_filename_component(pytest_name ${pytest} NAME) configure_file( ${pytest} "${CMAKE_BINARY_DIR}/${pytest_name}" COPYONLY) endforeach(pytest ${PY_TEST_SRCS}) find_package(PythonInterp ${MUSPECTRE_PYTHON_MAJOR_VERSION} REQUIRED) muSpectre_add_test(python_binding_test TYPE PYTHON python_binding_tests.py) if(${MPI_PARALLEL}) ############################################################################ # add MPI tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/mpi_test_*.cc") add_executable(mpi_main_test_suite tests/mpi_main_test_suite.cc ${TEST_SRCS}) target_link_libraries(mpi_main_test_suite ${Boost_LIBRARIES} muSpectre) muSpectre_add_test(mpi_main_test_suite1 TYPE BOOST MPI_NB_PROCS 1 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(mpi_main_test_suite2 TYPE BOOST MPI_NB_PROCS 2 mpi_main_test_suite --report_level=detailed) muSpectre_add_test(python_mpi_binding_test1 TYPE PYTHON MPI_NB_PROCS 1 python_mpi_binding_tests.py) muSpectre_add_test(python_mpi_binding_test2 TYPE PYTHON MPI_NB_PROCS 2 python_mpi_binding_tests.py) endif(${MPI_PARALLEL}) endif(${MAKE_TESTS}) ################################################################################ # compile the library -add_compile_options( -Werror) +# add_compile_options( -Werror) add_subdirectory( ${CMAKE_SOURCE_DIR}/src/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/language_bindings/ ) if (${MAKE_DOC_TARGET}) add_subdirectory( ${CMAKE_SOURCE_DIR}/doc/ ) endif() ################################################################################ if (${MAKE_EXAMPLES}) #compile executables set(binaries ${CMAKE_SOURCE_DIR}/bin/demonstrator1.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator_dynamic_solve.cc ${CMAKE_SOURCE_DIR}/bin/demonstrator2.cc ${CMAKE_SOURCE_DIR}/bin/hyper-elasticity.cc ${CMAKE_SOURCE_DIR}/bin/small_case.cc) if (${MPI_PARALLEL}) set (binaries ${binaries} ${CMAKE_SOURCE_DIR}/bin/demonstrator_mpi.cc ) endif (${MPI_PARALLEL}) foreach(binaryfile ${binaries}) get_filename_component(binaryname ${binaryfile} NAME_WE) add_executable(${binaryname} ${binaryfile}) target_link_libraries(${binaryname} ${Boost_LIBRARIES} muSpectre) endforeach(binaryfile ${binaries}) #or copy them file (GLOB pybins "${CMAKE_SOURCE_DIR}/bin/*.py") foreach(pybin ${pybins}) get_filename_component(binaryname ${pybin} NAME_WE) configure_file( ${pybin} "${CMAKE_BINARY_DIR}/${binaryname}.py" COPYONLY) endforeach(pybin ${pybins}) endif (${MAKE_EXAMPLES}) ################################################################################ # compile benchmarks if(${MAKE_BENCHMARKS}) file(GLOB benchmarks "${CMAKE_SOURCE_DIR}/benchmarks/benchmark*cc") foreach(benchmark ${benchmarks}) get_filename_component(benchmark_name ${benchmark} NAME_WE) add_executable(${benchmark_name} ${benchmark}) target_link_libraries(${benchmark_name} ${BOOST_LIBRARIES} muSpectre) endforeach(benchmark ${benchmark}) endif(${MAKE_BENCHMARKS}) diff --git a/src/io/CMakeLists.txt b/src/io/CMakeLists.txt index 0f311e9..8c42df1 100644 --- a/src/io/CMakeLists.txt +++ b/src/io/CMakeLists.txt @@ -1,39 +1,40 @@ #============================================================================== # file CMakeLists.txt # # @author Till Junge # # @date 20 Sep 2018 # # @brief configuration for compiling µSpectre's input/output facilities # # @section LICENSE # # Copyright © 2018 Till Junge # # µSpectre is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License as # published by the Free Software Foundation, either version 3, or (at # your option) any later version. # # µSpectre is distributed in the hope that it will be useful, but # WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # You should have received a copy of the GNU General Public License # along with GNU Emacs; see the file COPYING. If not, write to the # Free Software Foundation, Inc., 59 Temple Place - Suite 330, # Boston, MA 02111-1307, USA. # set (io_SRC "") if (${NETCDF}) set(io_SRC ${io_SRC} ${CMAKE_CURRENT_SOURCE_DIR}/mulib_input.cc + ${CMAKE_CURRENT_SOURCE_DIR}/mulib_output.cc ) endif (${NETCDF}) target_sources (muSpectre PRIVATE ${io_SRC}) diff --git a/src/io/mulib_input.hh b/src/io/mulib_input.hh index b10a0cc..0f8b00f 100644 --- a/src/io/mulib_input.hh +++ b/src/io/mulib_input.hh @@ -1,113 +1,115 @@ /** * file mulib_input.hh * * @author Till Junge * * @date 20 Sep 2018 * * @brief Class for handling the netcdf input file format for the minimal * muLib configuration. In this configuration, a linear elastic problem * is assumed, and the input file contains a per-pixel material label * map of the two- or three-dimensional microstructure, a list of * material properties in form of the elastic stiffness tensor as well * as a map between labels and material properties. The input file * format is netCDF3. * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MULIB_H #define MULIB_H #include "common/common.hh" #include "common/utilities.hh" #include #include #include #include namespace muSpectre { /** * Forward declaration */ template class CellBase; /** * thin wrapper around a netcdf file for input from µLib */ class MuLibInput { public: //! Default constructor MuLibInput() = delete; /** * Constructor from file path */ MuLibInput(filesystem::path path); //! Copy constructor MuLibInput(const MuLibInput &other) = delete; //! Move constructor MuLibInput(MuLibInput &&other) = default; //! Destructor virtual ~MuLibInput() = default; //! Copy assignment operator MuLibInput& operator=(const MuLibInput &other) = delete; //! Move assignment operator MuLibInput& operator=(MuLibInput &&other) = default; //! returns the spatial dimension of the problem defined in the file Dim_t get_dim() const; //! returns the spatial dimension of the problem defined in the file const std::vector & get_resolutions() const; //! returns the number of materials in the microstructure Dim_t get_nb_materials() const {return this->nb_materials;} /** * fills the materials for the provided Cell and initialises it */ template void setup_cell(CellBase & cell); + //! return input folder path + auto input_folder() const {return path.parent_path();}; protected: const filesystem::path path; netCDF::NcFile file; std::vector resolutions{}; Dim_t nb_materials{}; netCDF::NcVar mat_params{}; netCDF::NcVar micro_structure{}; std::vector mat_labels{}; private: }; } // muSpectre #endif /* MULIB_H */ diff --git a/src/io/mulib_output.cc b/src/io/mulib_output.cc index d6e120f..48f32b4 100644 --- a/src/io/mulib_output.cc +++ b/src/io/mulib_output.cc @@ -1,79 +1,149 @@ /** * file mulib_output.cc * * @author Martin Doskar * * @date 2 Nov 2018 * * @brief Implementation of the muLib output fixture * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "io/mulib_input.hh" -#include "solver/solvers.hh" +#include "io/mulib_output.hh" +#include "solver/solver_common.hh" #include "../external/base64.h" -#include "../external/json.h" - +#include "../external/json.hpp" + +#include +#include +#include +#include +#include +#include +#include using json = nlohmann::json; constexpr auto MESH_INDEX = 1; namespace muSpectre { - MuLibOutput::MuLibOutput(const MuLibInput & input) + MuLibOutput::MuLibOutput(const MuLibInput & input) : + outFolder{input.input_folder()}, resolutions{input.get_resolutions()} { - // Extract out folder from the inputs - // Extract resolution + // Verify resolution + if(!std::all_of(std::begin(this->resolutions), std::end(this->resolutions), [](const auto& r){ return r > 0;})) + throw std::invalid_argument{ "MuLibOutput::MuLibOutput: Non-positive resolution" }; + + if (!filesystem::exists(this->outFolder)) + filesystem::create_directory(this->outFolder); } - void MuLibOutput::store_resuls(const std::vector & results) const + void MuLibOutput::store_results(const std::vector & results) const { + // Define auxiliary constants + const auto nDim = static_cast(this->resolutions.size()); + const auto nTensorEntries = nDim*nDim; + const auto auxRes = std::array{ + static_cast(nDim > 0 ? this->resolutions.at(0) : 1), + static_cast(nDim > 1 ? this->resolutions.at(1) : 1), + static_cast(nDim > 2 ? this->resolutions.at(2) : 1) }; + const auto nLoci = std::accumulate(std::begin(auxRes), std::end(auxRes), 1, std::multiplies()); + const auto nLC = (nDim == 1) ? 1 : ((nDim == 2) ? 3 : ((nDim == 3) ? 6 : -1)); + + // Common part of JSON results files (-1 entries must be provided later) const auto outputTemplate = json{ {"MeshIndex", MESH_INDEX}, - {"Index", -1}, - {"TimeSteps", 1}, {"Location", "Cells"}, {"Compression", nullptr}, {"Encoding", { {"DataType", "Float64" }, - {"OriginalLength", -1 }, + {"OriginalLength", nLoci }, {"Offset", 0 }, - {"Length", -1 } + {"Length", nLoci } }} }; - auto data = std::vector{ 1.0, 2.0, 3.0 }; - - auto a = myJSONobj; - a["Index"] = 1; - a["FieldName"] = "Beta"; - a["ComponentName"] = "Beta"; - a["Encoding"]["Length"] = data.length(); - a["Encoding"]["OriginalLength"] = data.length(); - a["Data"] = base64::encode_vector_into_base64_string(data); - - + auto resultStrain = outputTemplate; + resultStrain["FieldName"] = "Strain"; + auto resultStress = outputTemplate; + resultStress["FieldName"] = "Stress"; + + // Save outputs load-case-wise + for(auto iLC = 0; iLC < nLC; ++iLC){ + + // Parse data into the correct format + auto strainData = std::vector>(nTensorEntries,std::vector{}); + auto stressData = std::vector>(nTensorEntries,std::vector{}); + + std::for_each(std::begin(strainData),std::end(strainData),[&nLoci](auto& v){ v.reserve(nLoci);}); + std::for_each(std::begin(stressData),std::end(stressData),[&nLoci](auto& v){ v.reserve(nLoci);}); + + for(auto iZ = 0; iZ < auxRes.at(2); ++iZ){ + for(auto iY = 0; iY < auxRes.at(1); ++iY){ + for(auto iX = 0; iX < auxRes.at(0); ++iX){ + const auto idx = iX*auxRes.at(2)*auxRes.at(1) + iY*auxRes.at(1) + iZ; + for (auto iE = 0; iE < nTensorEntries; ++iE){ + strainData.at(iE).push_back(static_cast(results.at(iLC).grad(idx*nTensorEntries+iE))); + stressData.at(iE).push_back(static_cast(results.at(iLC).stress(idx*nTensorEntries+iE))); + } + } + } + } + + // Write results JSON file + for (auto iE = 0; iE < nTensorEntries; ++iE){ + + const auto BCindex = 1; + const auto strainIndex = 1000 + (BCindex + 1)*100 + iLC*10 + iE; + const auto stressIndex = 2000 + (BCindex + 1)*100 + iLC*10 + iE; + + resultStrain["Index"] = strainIndex; + resultStrain["TimeSteps"] = iLC+1; + resultStrain["ComponentName"] = "e(" + std::to_string((iE % nDim) + 1) + ',' + std::to_string((iE / nDim) + 1) + ')'; + resultStrain["Data"] = base64::encode_vector_into_base64_string(strainData.at(iE)); + + resultStress["Index"] = stressIndex; + resultStress["TimeSteps"] = iLC+1; + resultStress["ComponentName"] = "s(" + std::to_string((iE % nDim) + 1) + ',' + std::to_string((iE / nDim) + 1) + ')'; + resultStress["Data"] = base64::encode_vector_into_base64_string(stressData.at(iE)); + + auto strainFile = std::ofstream{ this->outFolder / filesystem::path{ std::to_string(strainIndex) + ".result.json" } }; + if (strainFile.fail()) + throw std::runtime_error{ "muLibOutput::store_results: Unable to open strain result file." }; + strainFile << std::setw(4) << resultStrain; + strainFile.close(); + + auto stressFile = std::ofstream{ this->outFolder / filesystem::path{ std::to_string(stressIndex) + ".result.json" } }; + if (stressFile.fail()) + throw std::runtime_error{ "muLibOutput::store_results: Unable to open stress result file." }; + stressFile << std::setw(4) << resultStress; + stressFile.close(); + + } + + } } } \ No newline at end of file diff --git a/src/io/mulib_output.hh b/src/io/mulib_output.hh index c7af6b9..ea308b8 100644 --- a/src/io/mulib_output.hh +++ b/src/io/mulib_output.hh @@ -1,61 +1,61 @@ /** * file mulib_output.hh * * @author Martin Doskar * * @date 2 Nov 2018 * * @brief Class for extracting and saving local fields in the JSON format * suitable for muLib.Web visualisation. File names are fixed and * follow the convention of muLib.Web format. * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MULIB_OUT_H #define MULIB_OUT_H #include "common/common.hh" #include "common/utilities.hh" namespace muSpectre{ /** * Forward declaration */ class MuLibInput; struct OptimizeResult; /** * Thin class for storing homogenization results in JSON file for muLib */ class MuLibOutput { public: - MuLibOutput(const MuLibInput & input) explicit; + MuLibOutput(const MuLibInput & input); - void store_results(const OptimizeResult & results) const {}; + void store_results(const std::vector & results) const; private: filesystem::path outFolder{}; std::vector resolutions{}; }; } #endif /* MULIB_OUT_H */ \ No newline at end of file diff --git a/src/solver/mulib_solver.cc b/src/solver/mulib_solver.cc index 11ca97e..c817a1b 100644 --- a/src/solver/mulib_solver.cc +++ b/src/solver/mulib_solver.cc @@ -1,180 +1,185 @@ /** * @file muLib_solver.cc * * @author Till Junge * * @date 21 Sep 2018 * * @brief implementation of the mulib solver * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_linear_elastic_generic.hh" #include "solver/mulib_solver.hh" #include "cell/cell_factory.hh" #include "solver/solver_cg.hh" #include "solver/solvers.hh" +#include "io/mulib_output.hh" namespace muSpectre { namespace internal { template struct VoigtVec { }; template <> struct VoigtVec { static decltype(auto) get() { Eigen::Matrix mat{}; mat << 0, 0, 1, 1, 0, 1, 1, 0; return mat; } }; template <> struct VoigtVec { static decltype(auto) get() { Eigen::Matrix mat{}; mat << 0, 0, 1, 1, 2, 2, 1, 2, 0, 2, 0, 1, 2, 1, 2, 0, 1, 0; return mat; } }; } // internal template void mulib_worker(MuLibInput & file, Real newton_tol, Real equil_tol, Real cg_tol, Uint maxiter, Dim_t verbose) { Ccoord_t resolutions{}; for (auto && tup: akantu::zip(resolutions, file.get_resolutions())) { std::get<0>(tup) = std::get<1>(tup); } constexpr Rcoord_t lengths{CcoordOps::get_cube(1.)}; constexpr auto form{Formulation::small_strain}; auto cell {make_cell(resolutions, lengths, form)}; file.template setup_cell(cell); std::cout << newton_tol << equil_tol << verbose; using Delta_t = Eigen::MatrixXd; using Delta_vec = std::vector; Delta_vec Deltas{}; for (int i{0}; i < Dim; ++i) { for (int j{i}; j < Dim; ++j) { Eigen::MatrixXd Delta(Dim, Dim); Delta.setZero(); Delta(i, j) +=.5; Delta(j, i) +=.5; Deltas.push_back(Delta); } } SolverCG cg{cell, cg_tol, maxiter, bool(verbose)}; auto results = de_geus(cell, Deltas, cg, newton_tol, verbose); Delta_vec mean_delta_stresses{}; T4Mat C{}; C.setZero(); { size_t counter{}; for (int i{0}; i < Dim; ++i) { for (int j{i}; j < Dim; ++j, ++counter) { auto get_mean = [](const auto & vec) { Eigen::Map matrix(vec.data(), Dim*Dim, vec.rows()/(Dim*Dim)); auto& mat{ matrix.rowwise().mean()}; Eigen::Matrix mean; Eigen::Map>(mean.data()) = mat; return mean; }; const auto & opt_result = results[counter]; auto stress{get_mean(opt_result.stress)}; std::cout << std::endl<< "mean stress:" << std::endl << stress << std::endl; std::cout << "mean strain:" << std::endl << get_mean(opt_result.grad) << std::endl; for (int k{0}; k < Dim; ++k) { for (int l{0}; l < Dim; ++l) { get(C, i,j,k,l) = get(C, j,i,k,l) = stress(k,l); } } } } } std::cout << "stiffness tensor :" << std::endl << C << std::endl; Eigen::Matrix C_voigt{}; C_voigt.setZero(); const auto redirect{internal::VoigtVec::get()}; for (int I{0}; I < vsize(Dim); ++I) { const Dim_t & i{redirect(I,0)}; const Dim_t & j{redirect(I,1)}; for (int J{0}; J < vsize(Dim); ++J) { const Dim_t & k{redirect(J,0)}; const Dim_t & l{redirect(J,1)}; C_voigt(I,J) = get(C, i,j,k,l); } } std::cout << "stiffness tensor voigt :" << std::endl << C_voigt << std::endl; + // Store outputs + auto output = MuLibOutput{file}; + output.store_results(results); + } void mulib(const filesystem::path & path, Real newton_tol, Real equil_tol, Real cg_tol, Uint maxiter, Dim_t verbose) { MuLibInput file{path}; const auto dim{file.get_dim()}; switch (dim) { case twoD: { return mulib_worker(file, newton_tol, equil_tol, cg_tol, maxiter, verbose); break; } case threeD: { return mulib_worker(file, newton_tol, equil_tol, cg_tol, maxiter, verbose); break; } default: throw std::runtime_error("only two- and tree-dimensional cases considered"); break; } } } // muSpectre diff --git a/src/solver/mulib_solver.hh b/src/solver/mulib_solver.hh index 4274db0..8ef31df 100644 --- a/src/solver/mulib_solver.hh +++ b/src/solver/mulib_solver.hh @@ -1,46 +1,47 @@ /** * @file muLib_solver.hh * * @author Till Junge * * @date 21 Sep 2018 * * @brief Solver taking a mulib input file and returning homogenised properties * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MULIB_SOLVER_H #define MULIB_SOLVER_H #include "io/mulib_input.hh" +#include "io/mulib_output.hh" namespace muSpectre { void mulib(const filesystem::path & path, Real newton_tol, Real equil_tol, Real cg_tol, Uint maxiter, Dim_t verbose = 0); } // muSpectre #endif /* MULIB_SOLVER_H */