diff --git a/CMakeLists.txt b/CMakeLists.txt index fa08220..cde684a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,212 +1,213 @@ # ============================================================================= # 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(RUNNING_IN_CI "OFF" CACHE INTERNAL "changes output format for tests") 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(muspectreTools) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" STREQUAL "AppleClang") # using Clang add_compile_options(-Wno-missing-braces) elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # using GCC add_compile_options(-Wno-non-virtual-dtor) string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if ("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) +add_external_package(CGAL VERSION 4.12 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 tests file( GLOB TEST_SRCS "${CMAKE_SOURCE_DIR}/tests/test_*.cc") 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) ############################################################################## # 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_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/cmake/FindCGAL.cmake b/cmake/FindCGAL.cmake new file mode 100644 index 0000000..1091e8d --- /dev/null +++ b/cmake/FindCGAL.cmake @@ -0,0 +1,77 @@ +# - Find CGAL +# Find the CGAL includes and client library +# This module defines +# CGAL_INCLUDE_DIR, where to find CGAL.h +# CGAL_LIBRARIES, the libraries needed to use CGAL. +# CGAL_FOUND, If false, do not try to use CGAL. + +if(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) + set(CGAL_FOUND TRUE) + +else(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) + + FIND_PATH(CGAL_INCLUDE_DIR CGAL/basic.h + ${CGAL_ROOT}/include + /usr/include + /usr/local/include + $ENV{ProgramFiles}/CGAL/*/include + $ENV{SystemDrive}/CGAL/*/include + ) + + find_library(CGAL_LIBRARIES NAMES CGAL libCGAL + PATHS + ${CGAL_ROOT}/lib + /usr/lib + /usr/local/lib + /usr/lib/CGAL + /usr/lib64 + /usr/local/lib64 + /usr/lib64/CGAL + $ENV{ProgramFiles}/CGAL/*/lib + $ENV{SystemDrive}/CGAL/*/lib + ) + +# set(Boost_DEBUG ON) + find_package(Boost COMPONENTS thread REQUIRED) + if(Boost_FOUND) + set(BOOST_THREAD_LIBRARIES ${Boost_LIBRARIES}) + endif(Boost_FOUND) + + # check boost version we may need other components + if("${Boost_VERSION}" VERSION_GREATER "104900") + find_package(Boost COMPONENTS thread system REQUIRED) + if(Boost_FOUND) + set(BOOST_THREAD_LIBRARIES ${Boost_LIBRARIES}) + endif(Boost_FOUND) + endif("${Boost_VERSION}" VERSION_GREATER "104900") + + find_library(GMP_LIBRARIES NAMES gmp libgmp + PATHS + ${GMP_ROOT}/lib + /usr/lib + /usr/local/lib + /usr/lib/gmp + /usr/lib64 + /usr/local/lib64 + /usr/lib64/gmp + $ENV{ProgramFiles}/gmp/*/lib + $ENV{SystemDrive}/gmp/*/lib + ) + +message(STATUS "CGAL_INCLUDE_DIR=${CGAL_INCLUDE_DIR}") +message(STATUS "CGAL_LIBRARIES=${CGAL_LIBRARIES}") +message(STATUS "BOOST_THREAD_LIBRARIES=${BOOST_THREAD_LIBRARIES}") +message(STATUS "GMP_LIBRARIES=${GMP_LIBRARIES}") + + if(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) + set(CGAL_FOUND TRUE) + message(STATUS "Found CGAL: ${CGAL_INCLUDE_DIR}, ${CGAL_LIBRARIES}, ${BOOST_THREAD_LIBRARIES}, ${GMP_LIBRARIES}") + INCLUDE_DIRECTORIES(${CGAL_INCLUDE_DIR} $ENV{CGAL_CFG}) + else(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) + set(CGAL_FOUND FALSE) + message(STATUS "CGAL not found.") + endif(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) + + mark_as_advanced(CGAL_INCLUDE_DIR CGAL_LIBRARIES BOOST_THREAD_LIBRARIES GMP_LIBRARIES) + +endif(CGAL_INCLUDE_DIR AND CGAL_LIBRARIES AND BOOST_THREAD_LIBRARIES AND GMP_LIBRARIES) diff --git a/src/cell/CMakeFiles/CMakeDirectoryInformation.cmake b/src/cell/CMakeFiles/CMakeDirectoryInformation.cmake new file mode 100644 index 0000000..ad50ee4 --- /dev/null +++ b/src/cell/CMakeFiles/CMakeDirectoryInformation.cmake @@ -0,0 +1,16 @@ +# CMAKE generated file: DO NOT EDIT! +# Generated by "Unix Makefiles" Generator, CMake Version 3.11 + +# Relative path conversion top directories. +set(CMAKE_RELATIVE_PATH_TOP_SOURCE "/home/afalsafi/program-files/munew2/muSpectre") +set(CMAKE_RELATIVE_PATH_TOP_BINARY "/home/afalsafi/program-files/munew2/muSpectre") + +# Force unix paths in dependencies. +set(CMAKE_FORCE_UNIX_PATHS 1) + + +# The C and CXX include file regular expressions for this directory. +set(CMAKE_C_INCLUDE_REGEX_SCAN "^.*$") +set(CMAKE_C_INCLUDE_REGEX_COMPLAIN "^$") +set(CMAKE_CXX_INCLUDE_REGEX_SCAN ${CMAKE_C_INCLUDE_REGEX_SCAN}) +set(CMAKE_CXX_INCLUDE_REGEX_COMPLAIN ${CMAKE_C_INCLUDE_REGEX_COMPLAIN}) diff --git a/src/cell/CMakeFiles/progress.marks b/src/cell/CMakeFiles/progress.marks new file mode 100644 index 0000000..573541a --- /dev/null +++ b/src/cell/CMakeFiles/progress.marks @@ -0,0 +1 @@ +0 diff --git a/src/cell/CMakeLists.txt b/src/cell/CMakeLists.txt index 284ce60..3dba371 100644 --- a/src/cell/CMakeLists.txt +++ b/src/cell/CMakeLists.txt @@ -1,34 +1,35 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief configuration for cell implementations # # @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 (cells_SRC ${CMAKE_CURRENT_SOURCE_DIR}/cell_base.cc + ${CMAKE_CURRENT_SOURCE_DIR}/cell_split.cc ) target_sources(muSpectre PRIVATE ${cells_SRC}) diff --git a/src/cell/cell_base.cc b/src/cell/cell_base.cc index 3e69ebd..12ef290 100644 --- a/src/cell/cell_base.cc +++ b/src/cell/cell_base.cc @@ -1,426 +1,485 @@ -/** +/* * @file cell_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief Implementation for cell base class * * Copyright © 2017 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 "cell/cell_base.hh" #include "common/ccoord_operations.hh" #include "common/iterators.hh" #include "common/tensor_algebra.hh" +#include "common/common.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template - CellBase::CellBase(Projection_ptr projection_) + CellBase::CellBase(Projection_ptr projection_, SplittedCell is_cell_splitted) :subdomain_resolutions{projection_->get_subdomain_resolutions()}, subdomain_locations{projection_->get_subdomain_locations()}, domain_resolutions{projection_->get_domain_resolutions()}, pixels(subdomain_resolutions, subdomain_locations), domain_lengths{projection_->get_domain_lengths()}, fields{std::make_unique()}, F{make_field("Gradient", *this->fields)}, P{make_field("Piola-Kirchhoff-1", *this->fields)}, - projection{std::move(projection_)} + projection{std::move(projection_)}, + is_cell_splitted{is_cell_splitted} { } /** * turns out that the default move container in combination with * clang segfaults under certain (unclear) cicumstances, because the * move constructor of the optional appears to be busted in gcc * 7.2. Copying it (K) instead of moving it fixes the issue, and * since it is a reference, the cost is practically nil */ template CellBase::CellBase(CellBase && other): subdomain_resolutions{std::move(other.subdomain_resolutions)}, subdomain_locations{std::move(other.subdomain_locations)}, domain_resolutions{std::move(other.domain_resolutions)}, pixels{std::move(other.pixels)}, domain_lengths{std::move(other.domain_lengths)}, fields{std::move(other.fields)}, F{other.F}, P{other.P}, K{other.K}, // this seems to segfault under clang if it's not a move materials{std::move(other.materials)}, - projection{std::move(other.projection)} + projection{std::move(other.projection)}, + is_cell_splitted{is_cell_splitted} { } /* ---------------------------------------------------------------------- */ template typename CellBase::Material_t & CellBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_strain_vector() -> Vector_ref { return this->get_strain().eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_stress_vector() const -> ConstVector_ref { return this->get_stress().eigenvec(); } /* ---------------------------------------------------------------------- */ template void CellBase:: set_uniform_strain(const Eigen::Ref & strain) { this->F.get_map() = strain; } /* ---------------------------------------------------------------------- */ template auto CellBase::evaluate_stress() -> ConstVector_ref { if (not this->initialised) { this->initialise(); } for (auto & mat: this->materials) { mat->compute_stresses(this->F, this->P, this->get_formulation()); } return this->P.const_eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_stress_tangent() -> std::array { if (not this->initialised) { this->initialise(); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(this->F, this->P, this->K.value(), this->get_formulation()); } const TangentField_t & k = this->K.value(); return std::array{ this->P.const_eigenvec(), k.const_eigenvec()}; } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_projected_directional_stiffness (Eigen::Ref delF) -> Vector_ref { // the following const_cast should be safe, as long as the // constructed delF_field is const itself const TypedField delF_field ("Proxied raw memory for strain increment", *this->fields, Eigen::Map(const_cast(delF.data()), delF.size()), this->F.get_nb_components()); if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"δP; temp output for directional stiffness"}; auto & delP = this->get_managed_field(out_name); auto Kmap{this->K.value().get().get_map()}; auto delPmap{delP.get_map()}; MatrixFieldMap delFmap(delF_field); for (auto && tup: akantu::zip(Kmap, delFmap, delPmap)) { auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return Vector_ref(this->project(delP).data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template std::array CellBase::get_strain_shape() const { return this->projection->get_strain_shape(); } /* ---------------------------------------------------------------------- */ template void CellBase::apply_projection(Eigen::Ref vec) { TypedField field("Proxy for projection", *this->fields, vec, this->F.get_nb_components()); this->projection->apply_projection(field); } /* ---------------------------------------------------------------------- */ template typename CellBase::FullResponse_t - CellBase::evaluate_stress_tangent(StrainField_t & grad) { + CellBase:: + evaluate_stress_tangent(StrainField_t & grad) { if (this->initialised == false) { this->initialise(); } //! High level compatibility checks if (grad.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(grad, this->P, this->K.value(), - this->get_formulation()); + this->get_formulation(), this->is_cell_splitted); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::directional_stiffness(const TangentField_t &K, const StrainField_t &delF, StressField_t &delP) { for (auto && tup: akantu::zip(K.get_map(), delF.get_map(), delP.get_map())){ auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return this->project(delP); } /* ---------------------------------------------------------------------- */ template typename CellBase::Vector_ref CellBase::directional_stiffness_vec(const Eigen::Ref &delF) { if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); Vector_ref(in_tempref.data(), this->get_nb_dof()) = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return Vector_ref(out_tempref.data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template Eigen::ArrayXXd CellBase:: directional_stiffness_with_copy (Eigen::Ref delF) { if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); in_tempref.eigen() = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return out_tempref.eigen(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_strain() { if (this->initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename CellBase::StressField_t & CellBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ template const typename CellBase::TangentField_t & CellBase::get_tangent(bool create) { if (!this->K) { if (create) { this->K = make_field("Tangent Stiffness", *this->fields); } else { throw std::runtime_error ("K does not exist"); } } return this->K.value(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_managed_field(std::string unique_name) { if (!this->fields->check_field_exists(unique_name)) { return make_field(unique_name, *this->fields); } else { return static_cast(this->fields->at(unique_name)); } } /* ---------------------------------------------------------------------- */ template void CellBase::initialise(FFT_PlanFlags flags) { // check that all pixels have been assigned exactly one material this->check_material_coverage(); for (auto && mat: this->materials) { mat->initialise(); } // resize all global fields (strain, stress, etc) this->fields->initialise(this->subdomain_resolutions, this->subdomain_locations); // initialise the projection and compute the fft plan this->projection->initialise(flags); this->initialised = true; } /* ---------------------------------------------------------------------- */ template void CellBase::save_history_variables() { for (auto && mat: this->materials) { mat->save_history_variables(); } } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::begin() { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::end() { return this->pixels.end(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_adaptor() -> Adaptor { return Adaptor(*this); } /* ---------------------------------------------------------------------- */ template void CellBase::check_material_coverage() { auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector*> assignments(nb_pixels, nullptr); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); auto& assignment{assignments.at(index)}; if (assignment != nullptr) { std::stringstream err{}; err << "Pixel " << pixel << "is already assigned to material '" << assignment->get_name() << "' and cannot be reassigned to material '" << mat->get_name(); throw std::runtime_error(err.str()); } else { assignments[index] = mat.get(); } } } - // find and identify unassigned pixels std::vector unassigned_pixels; for (size_t i = 0; i < assignments.size(); ++i) { if (assignments[i] == nullptr) { - unassigned_pixels.push_back( - CcoordOps::get_ccoord(this->subdomain_resolutions, - this->subdomain_locations, i)); + unassigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, + this->subdomain_locations, i)); } } if (unassigned_pixels.size() != 0) { std::stringstream err {}; err << "The following pixels have were not assigned a material: "; for (auto & pixel: unassigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } } + template + Rcoord_t CellBase::get_pixel_lengths() const{ + auto nb_pixels = this->get_domain_resolutions(); + auto length_pixels = this->get_domain_lengths(); + Rcoord_t ret_val; + for (int i = 0 ; i < DimS ; i++ ){ + ret_val[i] = length_pixels[i] / nb_pixels[i]; + } + return ret_val; + } + + template + Rcoord_t CellBase::get_pixel_coordinate(Ccoord_t & pixel) const{ + auto pixel_length = this->get_pixel_lengths(); + Rcoord_t pixel_coordinate; + for (int i = 0 ; i < DimS; i++){ + pixel_coordinate[i] = pixel_length[i] * (pixel[i] + 0.5); + //0.5 is added to shift to the center of the pixel; + } + // return pixel_length * pixel; + return pixel_coordinate; + } + + template + bool CellBase::is_inside(Rcoord_t point){ + auto length_pixels = this->get_domain_lengths(); + Dim_t counter = 1; + for(int i = 0; i < DimS ; i++){ + if(point[i] < length_pixels[i]) { + counter++; + } + } + return counter == DimS; + } + + + template + bool CellBase::is_inside(Ccoord_t pixel){ + auto nb_pixels = this->get_domain_resolutions(); + Dim_t counter = 1; + for(int i = 0; i < DimS ; i++){ + if(pixel[i] < nb_pixels[i]) { + counter++; + } + } + return counter == DimS; + } + template + Real CellBase::get_pixel_volume(){ + auto pixel_length = get_pixel_lengths(); + Real Retval = 1.0; + for (int i = 0 ; i < DimS; i++){ + Retval*=pixel_length[i]; + } + return Retval; + } + template class CellBase; template class CellBase; } // muSpectre diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh index ac991ab..0a2be3c 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,511 +1,524 @@ /** * @file cell_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell cell with single * projection operator * * Copyright © 2017 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 CELL_BASE_H #define CELL_BASE_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include #include #include #include #include namespace muSpectre { /** * Cell adaptors implement the matrix-vector multiplication and * allow the system to be used like a sparse matrix in * conjugate-gradient-type solvers */ template class CellAdaptor; /** * Base class for cells that is not templated and therefore can be * in solvers that see cells as runtime-polymorphic objects. This * allows the use of standard * (i.e. spectral-method-implementation-agnostic) solvers, as for * instance the scipy solvers */ class Cell { public: //! sparse matrix emulation using Adaptor = CellAdaptor; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = Eigen::Matrix; //! dynamic matrix type for setting strains using Matrix_t = Eigen::Matrix; //! ref to constant vector using ConstVector_ref = Eigen::Map; //! output vector reference for solvers using Vector_ref = Eigen::Map; //! Default constructor Cell() = default; //! Copy constructor Cell(const Cell &other) = default; //! Move constructor Cell(Cell &&other) = default; //! Destructor virtual ~Cell() = default; //! Copy assignment operator Cell& operator=(const Cell &other) = default; //! Move assignment operator Cell& operator=(Cell &&other) = default; //! for handling double initialisations right bool is_initialised() const {return this->initialised;} //! returns the number of degrees of freedom in the cell virtual Dim_t get_nb_dof() const = 0; //! return the communicator object virtual const Communicator & get_communicator() const = 0; /** * formulation is hard set by the choice of the projection class */ virtual const Formulation & get_formulation() const = 0; /** * returns the material dimension of the problem */ virtual Dim_t get_material_dim() const = 0; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ virtual std::array get_strain_shape() const = 0; /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() = 0; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const = 0; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() = 0; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() = 0; /** * applies the projection operator in-place on the input vector */ virtual void apply_projection(Eigen::Ref vec) = 0; /** * freezes all the history variables of the materials */ virtual void save_history_variables() = 0; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) = 0; /** * set uniform strain (typically used to initialise problems */ virtual void set_uniform_strain(const Eigen::Ref &) = 0; //! get a sparse matrix view on the cell virtual Adaptor get_adaptor() = 0; protected: bool initialised{false}; //!< to handle double initialisation right private: }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase: public Cell { public: using Parent = Cell; using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< physical coordinates type //! global field collection using FieldCollection_t = GlobalFieldCollection; //! the collection is handled in a `std::unique_ptr` using Collection_ptr = std::unique_ptr; //! polymorphic base material type using Material_t = MaterialBase; //! materials handled through `std::unique_ptr`s using Material_ptr = std::unique_ptr; //! polymorphic base projection type using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; //! expected type for strain fields using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; //! iterator type over all cell pixel's using iterator = typename CcoordOps::Pixels::iterator; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = typename Parent::Vector_t; //! ref to constant vector using ConstVector_ref = typename Parent::ConstVector_ref; //! output vector reference for solvers using Vector_ref = typename Parent::Vector_ref; //! sparse matrix emulation using Adaptor = Parent::Adaptor; //! Default constructor CellBase() = delete; //! constructor using sizes and resolution - CellBase(Projection_ptr projection); + CellBase(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::no); //! Copy constructor CellBase(const CellBase &other) = delete; //! Move constructor CellBase(CellBase &&other); //! Destructor virtual ~CellBase() = default; //! Copy assignment operator CellBase& operator=(const CellBase &other) = delete; //! Move assignment operator CellBase& operator=(CellBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this cell */ Material_t & add_material(Material_ptr mat); /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() override; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const override; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() override; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() override; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc. (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) override; //! return the template param DimM (required for polymorphic use of `Cell` Dim_t get_material_dim() const override final {return DimM;} /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const override final; /** * applies the projection operator in-place on the input vector */ void apply_projection(Eigen::Ref vec) override final; /** * set uniform strain (typically used to initialise problems */ void set_uniform_strain(const Eigen::Ref &) override; /** - * evaluates all materials + * evaluate all materials */ - FullResponse_t evaluate_stress_tangent(StrainField_t & F); + virtual FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF or G:K:δε) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * vectorized version for eigen solvers, no copy, but only works * when fields have ArrayStore=false */ Vector_ref directional_stiffness_vec(const Eigen::Ref & delF); /** * Evaluate directional stiffness into a temporary array and * return a copy. This is a costly and wasteful interface to * directional_stiffness and should only be used for debugging or * in the python interface */ Eigen::ArrayXXd directional_stiffness_with_copy (Eigen::Ref delF); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); //! returns a ref to the cell's strain field StrainField_t & get_strain(); //! returns a ref to the cell's stress field const StressField_t & get_stress() const; //! returns a ref to the cell's tangent stiffness field const TangentField_t & get_tangent(bool create = false); + + //! returns a ref to a temporary field managed by the cell StrainField_t & get_managed_field(std::string unique_name); /** * general initialisation; initialises the projection and * fft_engine (i.e. infrastructure) but not the materials. These * need to be initialised separately */ void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, this function * calls this update for each material in the cell */ void save_history_variables() override final; iterator begin(); //!< iterator to the first pixel iterator end(); //!< iterator past the last pixel //! number of pixels in the cell size_t size() const {return pixels.size();} //! return the subdomain resolutions of the cell const Ccoord & get_subdomain_resolutions() const { return this->subdomain_resolutions;} //! return the subdomain locations of the cell const Ccoord & get_subdomain_locations() const { return this->subdomain_locations;} //! return the domain resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->domain_resolutions;} //! return the sizes of the cell const Rcoord & get_domain_lengths() const {return this->domain_lengths;} + //! return the real space coordinates of the center of pixel + Rcoord get_pixel_coordinate(Ccoord & pixel) const; + //! return the physical dimension of a pixel + Rcoord get_pixel_lengths() const; + //! check if the pixel is inside of the cell + bool is_inside(Rcoord point); + //! check if the point is inside of the cell + bool is_inside(Ccoord pixel); + //calculate the volume of each pixel: + Real get_pixel_volume(); /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const override final { return this->projection->get_formulation();} /** * get a reference to the projection object. should only be * required for debugging */ Eigen::Map get_projection() { return this->projection->get_operator();} //! returns the spatial size constexpr static Dim_t get_sdim() {return DimS;}; //! return a sparse matrix adaptor to the cell Adaptor get_adaptor() override; //! returns the number of degrees of freedom in the cell Dim_t get_nb_dof() const override {return this->size()*ipow(DimS, 2);}; //! return the communicator object virtual const Communicator & get_communicator() const override { return this->projection->get_communicator(); } protected: //! make sure that every pixel is assigned to one and only one material - void check_material_coverage(); + virtual void check_material_coverage(); const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions const Ccoord & domain_resolutions; //!< the cell's domain resolutions CcoordOps::Pixels pixels; //!< helper to iterate over the pixels const Rcoord & domain_lengths; //!< the cell's lengths Collection_ptr fields; //!< handle for the global fields of the cell StrainField_t & F; //!< ref to strain field StressField_t & P; //!< ref to stress field //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; //! container of the materials present in the cell std::vector materials{}; Projection_ptr projection; //!< handle for the projection operator - + bool initialised{false}; //!< to handle double initialisation right + SplittedCell is_cell_splitted ; private: }; /** * lightweight resource handle wrapping a `muSpectre::Cell` or * a subclass thereof into `Eigen::EigenBase`, so it can be * interpreted as a sparse matrix by Eigen solvers */ template class CellAdaptor: public Eigen::EigenBase> { public: using Scalar = double; //!< sparse matrix traits using RealScalar = double; //!< sparse matrix traits using StorageIndex = int; //!< sparse matrix traits enum { ColsAtCompileTime = Eigen::Dynamic, MaxColsAtCompileTime = Eigen::Dynamic, RowsAtCompileTime = Eigen::Dynamic, MaxRowsAtCompileTime = Eigen::Dynamic, IsRowMajor = false }; //! constructor CellAdaptor(Cell & cell):cell{cell}{} //!returns the number of logical rows Eigen::Index rows() const {return this->cell.get_nb_dof();} //!returns the number of logical columns Eigen::Index cols() const {return this->rows();} //! implementation of the evaluation template Eigen::Product operator*(const Eigen::MatrixBase& x) const { return Eigen::Product (*this, x.derived()); } Cell & cell; //!< ref to the cell }; } // muSpectre namespace Eigen { namespace internal { //! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` through a //! specialization of `Eigen::internal::generic_product_impl`: template struct generic_product_impl // GEMV stands for matrix-vector : generic_product_impl_base > { //! undocumented typedef typename Product::Scalar Scalar; //! undocumented template static void scaleAndAddTo(Dest& dst, const CellAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/) { // This method should implement "dst += alpha * lhs * rhs" inplace, // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, dst.noalias() += const_cast(lhs).cell. evaluate_projected_directional_stiffness(rhs); } }; } } #endif /* CELL_BASE_H */ diff --git a/src/cell/cell_factory.hh b/src/cell/cell_factory.hh index c826b78..161f775 100644 --- a/src/cell/cell_factory.hh +++ b/src/cell/cell_factory.hh @@ -1,161 +1,178 @@ /** * @file cell_factory.hh * * @author Till Junge * * @date 15 Dec 2017 * * @brief Cell factories to help create cells with ease * * Copyright © 2017 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 CELL_FACTORY_H #define CELL_FACTORY_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "cell/cell_base.hh" +#include "cell/cell_split.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/projection_small_strain.hh" #include "fft/fftw_engine.hh" #ifdef WITH_MPI #include "common/communicator.hh" #include "fft/fftwmpi_engine.hh" #endif #include namespace muSpectre { /** * Create a unique ptr to a Projection operator (with appropriate * FFT_engine) to be used in a cell constructor */ template > inline std::unique_ptr> cell_input(Ccoord_t resolutions, Rcoord_t lengths, Formulation form) { auto fft_ptr{ std::make_unique(resolutions, dof_for_formulation(form, DimS))}; switch (form) { case Formulation::finite_strain: { using Projection = ProjectionFiniteStrainFast; return std::make_unique(std::move(fft_ptr), lengths); break; } case Formulation::small_strain: { using Projection = ProjectionSmallStrain; return std::make_unique(std::move(fft_ptr), lengths); break; } default: { throw std::runtime_error("unknow formulation"); break; } } } /** * convenience function to create a cell (avoids having to build * and move the chain of unique_ptrs */ template , typename FFTEngine=FFTWEngine> inline Cell make_cell(Ccoord_t resolutions, Rcoord_t lengths, Formulation form) { auto && input = cell_input(resolutions, lengths, form); auto cell{Cell{std::move(input)}}; return cell; } + template , + typename FFTEngine=FFTWEngine> + inline + Cell make_cell_split(Ccoord_t resolutions, + Rcoord_t lengths, + Formulation form) { + + auto && input = cell_input(resolutions, lengths, + form); + Cell cell (std::move(input), SplittedCell::yes); + //CellBase cell_base {cell}; + return std::unique_ptr>(CellBase {cell}); + } + + #ifdef WITH_MPI /** * Create a unique ptr to a parallel Projection operator (with appropriate * FFT_engine) to be used in a cell constructor */ template > inline std::unique_ptr> parallel_cell_input(Ccoord_t resolutions, Rcoord_t lengths, Formulation form, const Communicator & comm) { auto fft_ptr{std::make_unique(resolutions, dof_for_formulation(form, DimM), comm)}; switch (form) { case Formulation::finite_strain: { using Projection = ProjectionFiniteStrainFast; return std::make_unique(std::move(fft_ptr), lengths); break; } case Formulation::small_strain: { using Projection = ProjectionSmallStrain; return std::make_unique(std::move(fft_ptr), lengths); break; } default: { throw std::runtime_error("unknown formulation"); break; } } } /** * convenience function to create a parallel cell (avoids having to build * and move the chain of unique_ptrs */ template , typename FFTEngine=FFTWMPIEngine> inline Cell make_parallel_cell(Ccoord_t resolutions, Rcoord_t lengths, Formulation form, const Communicator & comm) { auto && input = parallel_cell_input(resolutions, lengths, form, comm); auto cell{Cell{std::move(input)}}; return cell; } #endif /* WITH_MPI */ } // muSpectre #endif /* CELL_FACTORY_H */ diff --git a/src/cell/cell_split.cc b/src/cell/cell_split.cc new file mode 100644 index 0000000..2c79756 --- /dev/null +++ b/src/cell/cell_split.cc @@ -0,0 +1,173 @@ +/** + * @file cell_split.cc + * + * @author Ali Falsafi + * + * @date 19 Apr 2018 + * + * @brief Implementation for cell base class + * + * Copyright © 2017 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 "cell/cell_base.hh" +#include "cell/cell_split.hh" +#include "common/ccoord_operations.hh" +#include "common/iterators.hh" +#include "common/tensor_algebra.hh" + +#include +#include +#include + +namespace muSpectre { + + /* ---------------------------------------------------------------------- */ + + template + CellSplit::CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted ) + :Parent(std::move(projection), is_cell_splitted) {} + + /* ---------------------------------------------------------------------- */ + + template + std::vector CellSplit::get_intersected_ratios(){ + auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); + std::vector pixel_assigned_ratio(nb_pixels, 0.0); + for (auto & mat: this->materials) { + for (auto & pixel: *mat) { + auto index = CcoordOps::get_index(this->subdomain_resolutions, + this->subdomain_locations, + pixel); + pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel); + } + } + return pixel_assigned_ratio; + } + + + template + template + void CellSplit::complete_material_assignment(MaterialType& material, InternalArgs ... args){ + auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); + std::vector pixel_assigned_ratio(nb_pixels, 0.0); + for (auto & mat: this->materials) { + for (auto & pixel: *mat) { + auto index = CcoordOps::get_index(this->subdomain_resolutions, + this->subdomain_locations, + pixel); + pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel); + } + } + for (size_t i = 0; i < nb_pixels; ++i) { + if (pixel_assigned_ratio.at(i) < 1.0) { + material.add_pixel_split(CcoordOps::get_ccoord(this->subdomain_resolutions, + this->subdomain_locations, i), 1-pixel_assigned_ratio.at(i)); + } + } + } + + + + template + void CellSplit::check_material_coverage(){ + auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); + std::vector pixel_assigned_ratio(nb_pixels, 0.0); + for (auto & mat: this->materials) { + for (auto & pixel: *mat) { + auto index = CcoordOps::get_index(this->subdomain_resolutions, + this->subdomain_locations, + pixel); + pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel); + } + } + std::vector> over_assigned_pixels; + std::vector> under_assigned_pixels; + for (size_t i = 0; i < nb_pixels; ++i) { + if (pixel_assigned_ratio.at(i) > 1.0) { + over_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, + this->subdomain_locations, i)); + }else if (pixel_assigned_ratio.at(i) < 1.0) { + under_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, + this->subdomain_locations, i)); + } + } + if (over_assigned_pixels.size() != 0) { + std::stringstream err {}; + err << "Execesive material is assigned to the following pixels: "; + for (auto & pixel: over_assigned_pixels) { + err << pixel << ", "; + } + err << "and that cannot be handled"; + throw std::runtime_error(err.str()); + } + if (under_assigned_pixels.size() != 0) { + std::stringstream err {}; + err << "Insufficient material is assigned to the following pixels: "; + for (auto & pixel: under_assigned_pixels) { + err << pixel << ", "; + } + err << "and that cannot be handled"; + throw std::runtime_error(err.str()); + } + } + /* ---------------------------------------------------------------------- */ + //this piece of code handles the evaluation of stress an dtangent matrix + //in case the cells have materials in which pixels are partially composed of + //diffferent materials. + + template + typename CellSplit::FullResponse_t + CellSplit::evaluate_stress_tangent(StrainField_t & F){ + return evaluate_split_stress_tangent(F); + } + + template + typename CellSplit::FullResponse_t + CellSplit::evaluate_split_stress_tangent(StrainField_t & F){ + if (this->initialised == false) { + this->initialise(); + } + //! High level compatibility checks + if (F.size() != this->F.size()) { + throw std::runtime_error("Size mismatch"); + } + constexpr bool create_tangent{true}; + this->get_tangent(create_tangent); + + // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution + //of the partial influence of different materials assigend to each pixel. Therefore, this values should be + // initiialised as zero filled tensors + this->set_p_k_zero(); + + for (auto & mat: this->materials) { + mat->compute_stresses_tangent(F, this->P, this->K.value(), + this->get_formulation(), this->is_cell_splitted); + } + return std::tie(this->P, this->K.value()); + } + + template + void CellSplit::set_p_k_zero(){ + this->P.set_zero(); + this->K.value().get().set_zero(); + } + + template class CellSplit; + template class CellSplit; +} //muspectre diff --git a/src/cell/cell_split.hh b/src/cell/cell_split.hh new file mode 100644 index 0000000..25f516e --- /dev/null +++ b/src/cell/cell_split.hh @@ -0,0 +1,114 @@ +/** + * @file cell_split.hh + * + * @author Ali Falsafi + * + * @date 19 Apr 2018 + * + * @brief Base class representing a unit cell able to handle + * split material assignments + * + * Copyright © 2017 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 CELL_SPLIT_H +#define CELL_SPLIT_H + +#include "common/common.hh" +#include "common/ccoord_operations.hh" +#include "common/field.hh" +#include "common/utilities.hh" +#include "materials/material_base.hh" +#include "fft/projection_base.hh" +#include "cell/cell_traits.hh" +#include "cell/cell_base.hh" + +#include +#include +#include +#include +namespace muSpectre { + + //! DimS spatial dimension (dimension of problem + //! DimM material_dimension (dimension of constitutive law) + /* This class handles the cells that has splitly assigned material to their pixels */ + template + class CellSplit: public CellBase { + public: + using Parent = CellBase; //!< base class + //! global field collection + using FieldCollection_t = GlobalFieldCollection; + using Projection_t = ProjectionBase; + //! projections handled through `std::unique_ptr`s + using Projection_ptr = std::unique_ptr; + using StrainField_t = + TensorField; + //! expected type for stress fields + using StressField_t = + TensorField; + //! expected type for tangent stiffness fields + using TangentField_t = + TensorField; + //! combined stress and tangent field + using FullResponse_t = + std::tuple; + + + //! Default constructor + CellSplit() = delete; + + //! constructor using sizes and resolution + CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::yes); + + //! Copy constructor + CellSplit(const CellSplit &other) = delete; + + //! Move constructor + CellSplit(CellSplit &&other) = default; + + //! Destructor + virtual ~CellSplit() = default; + + //! Copy assignment operator + CellSplit& operator=(const CellSplit &other) = delete; + + //! Move assignment operator + CellSplit& operator=(CellSplit &&other) = default; + + + //completes the assignmnet of material with a specific material so all the under-assigned pixels would be assigned to a material + template + void complete_material_assignment(MaterialType& material, InternalArgs ... args); + + //returns the assigend ratios to each pixel + std::vector get_intersected_ratios(); + + protected: + void check_material_coverage() override final; + void set_p_k_zero(); + //full resppnse is consisted of the stresses and tangent matrix + FullResponse_t evaluate_stress_tangent(StrainField_t & F) override final; + FullResponse_t evaluate_split_stress_tangent(StrainField_t & F); + private: + }; + + + +} + +#endif /* CELL_SPLIT_H */ diff --git a/src/common/CMakeLists.txt b/src/common/CMakeLists.txt index 0545f5b..2d2c383 100644 --- a/src/common/CMakeLists.txt +++ b/src/common/CMakeLists.txt @@ -1,36 +1,37 @@ # ============================================================================= # file CMakeLists.txt # # @author Till Junge # # @date 08 Jan 2018 # # @brief configuration for files in common # # @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 (common_SRC #${CMAKE_CURRENT_SOURCE_DIR}/units.cc ${CMAKE_CURRENT_SOURCE_DIR}/common.cc + ${CMAKE_CURRENT_SOURCE_DIR}/intersection_octree.cc ${CMAKE_CURRENT_SOURCE_DIR}/voigt_conversion.cc ) target_sources(muSpectre PRIVATE ${common_SRC}) diff --git a/src/common/common.hh b/src/common/common.hh index 9db6790..d30229f 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,311 +1,317 @@ /** * @file common.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Small definitions of commonly used types throughout µSpectre * * @section LICENSE * * Copyright © 2017 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 #include #include #include #include #include #ifndef COMMON_H #define COMMON_H namespace muSpectre { /** * Eigen uses signed integers for dimensions. For consistency, µSpectre uses them througout the code. needs to represent -1 for eigen */ using Dim_t = int; constexpr Dim_t oneD{1}; //!< constant for a one-dimensional problem constexpr Dim_t twoD{2}; //!< constant for a two-dimensional problem constexpr Dim_t threeD{3}; //!< constant for a three-dimensional problem constexpr Dim_t firstOrder{1}; //!< constant for vectors constexpr Dim_t secondOrder{2}; //!< constant second-order tensors constexpr Dim_t fourthOrder{4}; //!< constant fourth-order tensors //@{ //! @anchor scalars //! Scalar types used for mathematical calculations using Uint = unsigned int; using Int = int; using Real = double; using Complex = std::complex; //@} //! Ccoord_t are cell coordinates, i.e. integer coordinates template using Ccoord_t = std::array; //! Real space coordinates template using Rcoord_t = std::array; /** * Allows inserting `muSpectre::Ccoord_t` and `muSpectre::Rcoord_t` * into `std::ostream`s */ template std::ostream & operator << (std::ostream & os, const std::array & index) { os << "("; for (size_t i = 0; i < dim-1; ++i) { os << index[i] << ", "; } os << index.back() << ")"; return os; } //! element-wise division template Rcoord_t operator/(const Rcoord_t & a, const Rcoord_t & b) { Rcoord_t retval{a}; for (size_t i = 0; i < dim; ++i) { retval[i]/=b[i]; } return retval; } //! element-wise division template Rcoord_t operator/(const Rcoord_t & a, const Ccoord_t & b) { Rcoord_t retval{a}; for (size_t i = 0; i < dim; ++i) { retval[i]/=b[i]; } return retval; } //! convenience definitions constexpr Real pi{3.1415926535897932384626433}; //! compile-time potentiation required for field-size computations template constexpr R ipow(R base, I exponent) { static_assert(std::is_integral::value, "Type must be integer"); R retval{1}; for (I i = 0; i < exponent; ++i) { retval *= base; } return retval; } /** * Copyright banner to be printed to the terminal by executables * Arguments are the executable's name, year of writing and the name * + address of the copyright holder */ void banner(std::string name, Uint year, std::string cpy_holder); /** * Planner flags for FFT (follows FFTW, hopefully this choice will * be compatible with alternative FFT implementations) * @enum muSpectre::FFT_PlanFlags */ enum class FFT_PlanFlags { estimate, //!< cheapest plan for slowest execution measure, //!< more expensive plan for fast execution patient //!< very expensive plan for fastest execution }; //! continuum mechanics flags enum class Formulation{ finite_strain, //!< causes evaluation in PK1(F) small_strain, //!< causes evaluation in σ(ε) small_strain_sym //!< symmetric storage as vector ε }; + //! split cell flags + enum class SplittedCell{ + yes, + no + }; + /** * compile time computation of voigt vector */ template constexpr Dim_t vsize(Dim_t dim) { if (sym) { return (dim * (dim - 1) / 2 + dim); } else { return dim*dim; } } //! compute the number of degrees of freedom to store for the strain //! tenor given dimension dim constexpr Dim_t dof_for_formulation(const Formulation form, const Dim_t dim) { switch (form) { case Formulation::small_strain_sym: { return vsize(dim); break; } default: return ipow(dim, 2); break; } } //! inserts `muSpectre::Formulation`s into `std::ostream`s std::ostream & operator<<(std::ostream & os, Formulation f); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of stress measure they provide, //! and µSpectre will handle conversions enum class StressMeasure { Cauchy, //!< Cauchy stress σ PK1, //!< First Piola-Kirchhoff stress PK2, //!< Second Piola-Kirchhoff stress Kirchhoff, //!< Kirchhoff stress τ Biot, //!< Biot stress Mandel, //!< Mandel stress no_stress_ //!< only for triggering static_asserts }; //! inserts `muSpectre::StressMeasure`s into `std::ostream`s std::ostream & operator<<(std::ostream & os, StressMeasure s); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of strain measure they require and //! µSpectre will provide it enum class StrainMeasure { Gradient, //!< placement gradient (δy/δx) Infinitesimal, //!< small strain tensor .5(∇u + ∇uᵀ) GreenLagrange, //!< Green-Lagrange strain .5(Fᵀ·F - I) Biot, //!< Biot strain Log, //!< logarithmic strain Almansi, //!< Almansi strain RCauchyGreen, //!< Right Cauchy-Green tensor LCauchyGreen, //!< Left Cauchy-Green tensor no_strain_ //!< only for triggering static_assert }; //! inserts `muSpectre::StrainMeasure`s into `std::ostream`s std::ostream & operator<<(std::ostream & os, StrainMeasure s); /* ---------------------------------------------------------------------- */ /** * all isotropic elastic moduli to identify conversions, such as E * = µ(3λ + 2µ)/(λ+µ). For the full description, see * https://en.wikipedia.org/wiki/Lam%C3%A9_parameters * Not all the conversions are implemented, so please add as needed */ enum class ElasticModulus { Bulk, //!< Bulk modulus K K = Bulk, //!< alias for ``ElasticModulus::Bulk`` Young, //!< Young's modulus E E = Young, //!< alias for ``ElasticModulus::Young`` lambda, //!< Lamé's first parameter λ Shear, //!< Shear modulus G or µ G = Shear, //!< alias for ``ElasticModulus::Shear`` mu = Shear, //!< alias for ``ElasticModulus::Shear`` Poisson, //!< Poisson's ratio ν nu = Poisson, //!< alias for ``ElasticModulus::Poisson`` Pwave, //!< P-wave modulus M M=Pwave, //!< alias for ``ElasticModulus::Pwave`` no_modulus_}; //!< only for triggering static_asserts /** * define comparison in order to exploit that moduli can be * expressed in terms of any two other moduli in any order (e.g. K * = K(E, ν) = K(ν, E) */ constexpr inline bool operator<(ElasticModulus A, ElasticModulus B) { return static_cast(A) < static_cast(B); } /* ---------------------------------------------------------------------- */ /** Compile-time function to g strain measure stored by muSpectre depending on the formulation **/ constexpr StrainMeasure get_stored_strain_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StrainMeasure::Gradient; break; } case Formulation::small_strain: { return StrainMeasure::Infinitesimal; break; } default: return StrainMeasure::no_strain_; break; } } /** Compile-time function to g stress measure stored by muSpectre depending on the formulation **/ constexpr StressMeasure get_stored_stress_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StressMeasure::PK1; break; } case Formulation::small_strain: { return StressMeasure::Cauchy; break; } default: return StressMeasure::no_stress_; break; } } /* ---------------------------------------------------------------------- */ /** Compile-time functions to get the stress and strain measures after they may have been modified by choosing a formulation. For instance, a law that expecs a Green-Lagrange strain as input will get the infinitesimal strain tensor instead in a small strain computation **/ constexpr StrainMeasure get_formulation_strain_type(Formulation form, StrainMeasure expected) { switch (form) { case Formulation::finite_strain: { return expected; break; } case Formulation::small_strain: { return get_stored_strain_type(form); break; } default: return StrainMeasure::no_strain_; break; } } } // muSpectre #ifndef EXPLICITLY_TURNED_ON_CXX17 #include "common/utilities.hh" #endif #endif /* COMMON_H */ diff --git a/src/common/intersection_octree.cc b/src/common/intersection_octree.cc new file mode 100644 index 0000000..a6cbd29 --- /dev/null +++ b/src/common/intersection_octree.cc @@ -0,0 +1,191 @@ +/** + * @file intersection_octree.cc + * + * @author Ali Falsafi + * + * @date May 2018 + * + * @brief common operations on pixel addressing + * + * Copyright © 2018 Ali Falsafi + * + * µ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 +#include + +#include "common/ccoord_operations.hh" +#include "common/common.hh" +#include "cell/cell_base.hh" +#include "materials/material_base.hh" +#include "common/intersection_octree.hh" +#include "common/intersection_volume_calculator.hh" + +namespace muSpectre { + + /* ---------------------------------------------------------------------- */ + + template + Node::Node(const Rcoord &new_origin, const Ccoord &new_lengths, + int depth, RootNode_t& root, bool is_root ) + :root_node(root), + origin(new_origin), + Clengths(new_lengths), + depth(depth), + is_pixel((depth == root.max_depth)), + children_no (((is_pixel)? 0 : pow(2, DimS))) + { + for (int i = 0 ; i < DimS ; i++){ + this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; + } + if (not is_root){ + this->check_node(); + } + } + + template + void Node::check_node() { + Real total_volume = 1.0, intersected_volume = 0.0, intersection_ratio = 0.0; + constexpr Real machine_prescision = 1e-13; + for (int i = 0; i < DimS; i++){ + total_volume *= this->Rlengths[i]; + } + + // this volume should be calculated by CGAL as the intersection volume of the precipitate and the Node + intersected_volume = total_volume ; + intersected_volume = intersection_volume_calculator(this->root_node.precipitate_vertices, + this->origin, this->Rlengths); + + intersection_ratio = intersected_volume/ total_volume; + + if (intersection_ratio > (1.0 - machine_prescision) && intersection_ratio < (1.0 + machine_prescision)) { + // all the pixels in this node should be assigned to the precipitate material + Real pix_num = pow(2, (this->root_node.max_depth - this->depth)); + Ccoord origin_point, pixels_number; + for (int i = 0; i < DimS; i++ ) { + origin_point[i] = this->origin.at(i) / this->root_node.pixel_lengths.at(i); + pixels_number[i] = pix_num; + } + + CcoordOps::Pixels pixels (pixels_number, origin_point); + for (auto && pix:pixels){ + this->root_node.intersected_pixels.push_back(pix); + this->root_node.intersection_ratios.push_back(intersection_ratio); + } + }else if (intersection_ratio > machine_prescision) { + this->split_node(intersection_ratio); + } + } + + + + template + void Node::split_node(Real intersection_ratio){ + if (this->depth == this->root_node.max_depth) { + Ccoord pixel; + for (int i = 0 ; i < DimS ; i++){ + pixel[i] = this->origin[i] / this->root_node.pixel_lengths[i]; + } + this->root_node.intersected_pixels.push_back(pixel); + this->root_node.intersection_ratios.push_back(intersection_ratio); + } else { + this->divide_node(); + } + } + + + template + void Node::divide_node(){ + Rcoord new_origin; + Ccoord new_length; + this->children.reserve(children_no); + CcoordOps::Pixels sub_nodes (CcoordOps::get_cube(2)); + for (auto && sub_node: sub_nodes ) { + for (int i = 0 ; i< DimS ; i++){ + new_length[i] = (int) this->Clengths[i] *0.5; + new_origin[i] = this->origin[i] + sub_node[i] * Rlengths[i] * 0.5; + } + this->children.emplace_back(new_origin, new_length, this->depth+1, this->root_node, false); + } + } + + + template + RootNode::RootNode(CellBase& cell, + std::vector vert_precipitate) + /*Calling parent constructing method simply by argumnets of (coordinates of origin, 2^depth_max in each direction, 0 ,*this)*/ + + :Node(Rcoord{}, + CcoordOps::get_cube(pow(2, + std::ceil(std::log2(cell.get_domain_resolutions().at(std::distance(std::max_element + (cell.get_domain_resolutions().begin(), + cell.get_domain_resolutions().end()) + ,cell.get_domain_resolutions().begin()) + ) + ) + ) + ) + ), + 0, *this, true), + + cell(cell), + cell_length (cell.get_domain_lengths()), + pixel_lengths (cell.get_pixel_lengths()), + cell_resolution (cell.get_domain_resolutions()), + max_resolution (this->cell_resolution.at(std::distance(std::max_element(this->cell_resolution.begin(), + this->cell_resolution.end()) + ,this->cell_resolution.begin()))), + max_depth (std::ceil(std::log2(this->max_resolution))), + precipitate_vertices (vert_precipitate) + { + for (int i = 0 ; i < DimS ; i++){ + this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; + } + this->check_node(); + } + + + template + std::vector> RootNode::get_intersected_pixels(){ + return this->intersected_pixels; + } + + template + std::vector RootNode::get_intersection_ratios(){ + return this->intersection_ratios; + } + + + template + void RootNode:: check_root_node() + { + for (int i = 0 ; i < DimS ; i++){ + this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; + } + this->check_node(); + } + + template class RootNode; + template class RootNode; + + template class Node; + template class Node; + + +} //muSpctre + diff --git a/src/common/intersection_octree.hh b/src/common/intersection_octree.hh new file mode 100644 index 0000000..d9a39fe --- /dev/null +++ b/src/common/intersection_octree.hh @@ -0,0 +1,145 @@ +/** + * @file intersection_octree.hh + * + * @author Ali Falsafi + * + * @date May 2018 + * + * @brief common operations on pixel addressing + * + * Copyright © 2018 Ali Falsafi + * + * µ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 INTERSECTION_OCTREE_H +#define INTERSECTION_OCTREE_H + +#include + +#include "common/ccoord_operations.hh" +#include "common/common.hh" +#include "cell/cell_base.hh" +#include "materials/material_base.hh" +#include "common/intersection_volume_calculator.hh" + + + +namespace muSpectre { + + + template + class RootNode; + + template + class Node + { + public: + using Rcoord = Rcoord_t; //!< physical coordinates type + using Ccoord = Ccoord_t; //!< cell coordinates type + using RootNode_t = RootNode; + //! Default constructor + Node() = delete; + + // Construct by origin, lenghts, and depth + Node(const Rcoord &new_origin, const Ccoord &new_lenghts, + int depth, RootNode_t& root, bool is_root ); + + // Construct by cell + Node(const CellBase & cell, RootNode_t& root); + + //! Copy constructor + Node(const Node &other) = delete; + + //! Move constructor + Node(Node &&other) = default; + + //! Destructor + ~Node() = default; + + //This function checks the status of the node and orders its devision into + //smaller nodes or asssign material to it + virtual void check_node(); + + //This function gives the ratio of the node which happens to be inside the precipitate + //and assign materila to it if node is a pixel or divide it furhter if it is not + void split_node(Real ratio); + + //this function constructs children of a node + void divide_node(); + + protected: + //////// + RootNode_t& root_node; + Rcoord origin, Rlengths{}; + Ccoord Clengths{}; + int depth; + bool is_pixel; + int children_no; + std::vector children{}; + }; + + + + template + class RootNode: public Node + { + friend class Node; + public: + using Rcoord = Rcoord_t; //!< physical coordinates type + using Ccoord = Ccoord_t; //!< cell coordinates type + using Parent = Node; //!< base class + + //!Default Constructor + RootNode() = delete ; + + //! Constructing a root node for a cell and a preticipate inside that cell + RootNode(CellBase& cell, + std::vector vert_precipitate); + + //! Copy constructor + RootNode(const RootNode &other) = delete; + + //! Move constructor + RootNode(RootNode &&other) = default; + + //! Destructor + ~RootNode() = default; + + //returns the pixels which have intersection raio with the preipitate + std::vector get_intersected_pixels (); + + //return the intersection ratios of corresponding to the pixels returned by get_intersected_pixels() + std::vector get_intersection_ratios (); + + //checking rootnode: + void check_root_node(); + + protected: + CellBase& cell; + Rcoord cell_length, pixel_lengths; + Ccoord cell_resolution; + int max_resolution ; + int max_depth; + std::vector precipitate_vertices{}; + std::vector intersected_pixels{}; + std::vector intersection_ratios{}; + }; + +} //muSpectre + +#endif /* INTERSECTION_OCTREE_H */ diff --git a/src/common/intersection_volume_calculator.hh b/src/common/intersection_volume_calculator.hh new file mode 100644 index 0000000..107d024 --- /dev/null +++ b/src/common/intersection_volume_calculator.hh @@ -0,0 +1,206 @@ +/** + * @file intersection_volume_calculator.hh + * + * @author Ali Falsafi + * + * @date 04 June 2018 + * + * @brief common operations on pixel addressing + * + * Copyright © 2018 Ali Falsafi + * + * µ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 INTERSECTION_VOLUME_CALCULATOR_H +#define INTERSECTION_VOLUME_CALCULATOR_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "common/ccoord_operations.hh" +#include "common/common.hh" + +#include +#include +#include + +typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Nef_polyhedron_3 Nef_polyhedron; +typedef Polyhedron::HalfedgeDS HalfedgeDS; +typedef Kernel::Point_3 Point_3; +typedef Kernel::Vector_3 Vector_3; +typedef Kernel::Tetrahedron_3 Tetrahedron; +typedef Polyhedron::Vertex_iterator Vertex_iterator; +typedef Polyhedron::Facet_iterator Facet_iterator; +typedef Polyhedron::Halfedge_around_facet_circulator Hafc; +typedef CGAL::Aff_transformation_3 Transformation_mat; +using Dim_t = int; +using Real = double; + + + +Polyhedron node_polyhedron_generator(std::array origin, std::array lengths) { + std::vector node_vertices ; + node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] )); + node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] )); + node_vertices.push_back(Point_3(origin[0] , origin[1] + lengths[1], origin[2] )); + node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] + lengths[2])); + node_vertices.push_back(Point_3(origin[0] , origin[1] + lengths[1], origin[2] + lengths[2])); + node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] + lengths[2])); + node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] )); + node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] + lengths[2])); + + Polyhedron node_poly; + CGAL::convex_hull_3(node_vertices.begin(), node_vertices.end(), node_poly); + return node_poly; +} + + +Polyhedron convex_polyhedron_generator(std::vector> convex_poly_vertices) { + std::vector poly_vertices ; + for (auto && point:convex_poly_vertices ){ + poly_vertices.push_back(Point_3(point[0], point[1], point[2])); + } + Polyhedron convex_poly; + CGAL::convex_hull_3(poly_vertices.begin(), poly_vertices.end(), convex_poly); + return convex_poly; +} + + + +template +class Correction{ +public: + std::array correct_origin (std::array array); + std::array correct_length (std::array array ); + std::vector> correct_vector (std::vector> vector); +}; + + + +template<> +class Correction <3> { +public: + std::array correct_origin (std::array array){ + return array; + } + + std::array correct_length (std::array array){ + return array; + } + std::vector> correct_vector(std::vector> vertices){ + std::vector> corrected_convex_poly_vertices; + return vertices; + } +}; + + +template<> +class Correction <2> { +public: + std::vector> correct_vector(std::vector> vertices){ + std::vector> corrected_convex_poly_vertices; + for (auto && vertice : vertices){ + corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 0.0}); + corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 1.0}); + } + return corrected_convex_poly_vertices; + } + std::array correct_origin (std::array array){ + return std::array{array.at(0),array.at(1),0.0}; + } + + std::array correct_length (std::array array){ + return std::array{array.at(0),array.at(1),1.0}; + } +}; + + + +template +Real intersection_volume_calculator(std::vector> convex_poly_vertices, + std::array origin, std::array lengths){ + + Correction correction; + + std::vector> corrected_convex_poly_vertices (correction.correct_vector(convex_poly_vertices)); + std::array corrected_origin(correction.correct_origin(origin)); + std::array corrected_lengths(correction.correct_length(lengths)); + + + std::vector result_vertex ; + std::vector a_result_facet; + Vector_3 result_vector ; + Point_3 center_point = Point_3( 0.0, 0.0, 0.0) ; + + Point_3 origin_coordinate_point (0,0,0), tet_vertices[4]; + Real return_volume = 0.0; + Polyhedron node_poly, convex_poly, result_poly; + Tetrahedron tetra; + node_poly = node_polyhedron_generator(corrected_origin, corrected_lengths); + convex_poly = convex_polyhedron_generator(corrected_convex_poly_vertices); + Nef_polyhedron convex_poly_nef, node_poly_nef, result_nef ; + convex_poly_nef = Nef_polyhedron (convex_poly); + node_poly_nef = Nef_polyhedron (node_poly); + + //computing the intersection: + result_nef = convex_poly_nef * node_poly_nef ; + + if(result_nef.is_simple() and not result_nef.is_empty()) { + result_nef.convert_to_Polyhedron(result_poly); + + //extracting vertex + int v_id = 0 ; + for ( Vertex_iterator v = result_poly.vertices_begin(); v != result_poly.vertices_end(); ++v){ + result_vertex.push_back(v->point()); + result_vector = v->point() - origin_coordinate_point; + center_point = center_point + result_vector; + v->id() = v_id++; + } + + //calculating the center point of the resultant polyhedron: + Transformation_mat scale (1,0,0,0,0,1,0,0,0,0,1,0,v_id); + center_point = scale(center_point); + tet_vertices[0] = center_point; + + int f_id = 0; + for ( Facet_iterator f = result_poly.facets_begin(); f != result_poly.facets_end(); ++f){ + a_result_facet.clear(); + f->id() = f_id++; + Hafc circ = f->facet_begin(); + //assignig corresponding vertex to facets: + do { + a_result_facet.push_back(circ-> vertex()-> id()); + } while ( ++circ != f->facet_begin()); + //making tetreahedrals to calculate volumes: + for(int ii = 0 ; ii<3; ii++) tet_vertices[ii+1] = result_vertex[a_result_facet[ii]]; + //calculating volumes: + tetra = Tetrahedron (tet_vertices[0], tet_vertices[1], tet_vertices[2], tet_vertices[3]); + return_volume += CGAL::to_double(tetra.volume()); + } + } + return return_volume; +} + +#endif //INTERSECTION_VOLUME_CALCULATOR diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc index 6287595..7a956c1 100644 --- a/src/materials/material_base.cc +++ b/src/materials/material_base.cc @@ -1,79 +1,86 @@ /** * @file material_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief implementation of materi * * Copyright © 2017 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_base.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialBase::MaterialBase(std::string name) - :name(name) { + :name(name),assigned_ratio{make_field("assigned ratio", + this->internal_fields)}, + assigned_ratio_mapped{assigned_ratio}{ static_assert((DimM == oneD)|| (DimM == twoD)|| (DimM == threeD), "only 1, 2, or threeD supported"); } /* ---------------------------------------------------------------------- */ template const std::string & MaterialBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel(const Ccoord &ccoord) { this->internal_fields.add_pixel(ccoord); } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses(const Field_t & F, Field_t & P, - Formulation form) { + Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses(StrainField_t::check_ref(F), StressField_t::check_ref(P), - form); + form, is_cell_splitted); + } + + template + Real MaterialBase::get_assigned_ratio(Ccoord pixel){ + return this->assigned_ratio_mapped[pixel]; } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, - Formulation form) { + Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses_tangent(StrainField_t::check_ref(F), StressField_t::check_ref(P), TangentField_t::check_ref(K), - form); + form, is_cell_splitted); } template class MaterialBase<2, 2>; template class MaterialBase<2, 3>; template class MaterialBase<3, 3>; } // muSpectre diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 9d1365c..64be10c 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,161 +1,165 @@ /** * @file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * Copyright © 2017 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 MATERIAL_BASE_H #define MATERIAL_BASE_H #include "common/common.hh" #include "common/field.hh" #include "common/field_collection.hh" #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) and /** * @a DimM is the material dimension (i.e., the dimension of constitutive * law; even for e.g. two-dimensional problems the constitutive law could * live in three-dimensional space for e.g. plane strain or stress problems) */ template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for cell-wide fields, like stress, strain, etc using GFieldCollection_t = GlobalFieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = LocalFieldCollection; using iterator = typename MFieldCollection_t::iterator; //!< pixel iterator //! polymorphic base class for fields only to be used for debugging using Field_t = internal::FieldBase; //! Full type for stress fields using StressField_t = TensorField; //! Full type for strain fields using StrainField_t = StressField_t; //! Full type for tangent stiffness fields fields using TangentField_t = TensorField; using Ccoord = Ccoord_t; //!< cell coordinates type //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) = delete; //! Destructor virtual ~MaterialBase() = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) = delete; /** * take responsibility for a pixel identified by its cell coordinates * WARNING: this won't work for materials with additional info per pixel * (as, e.g. for eigenstrain), we need to pass more parameters. Materials * of this tye need to overload add_pixel */ virtual void add_pixel(const Ccoord & ccooord); //! allocate memory, etc, but also: wipe history variables! virtual void initialise() = 0; /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, the virtual * base implementation does nothing, but materials with history * variables need to implement this */ virtual void save_history_variables() {}; //! return the material's name const std::string & get_name() const; //! spatial dimension for static inheritance constexpr static Dim_t sdim() {return DimS;} //! material dimension for static inheritance constexpr static Dim_t mdim() {return DimM;} //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, - Formulation form) = 0; + Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses, mostly for debugging and * testing. Has runtime-cost associated with compatibility-checking and * conversion of the Field_t arguments that can be avoided by using the * version with strongly typed field references */ void compute_stresses(const Field_t & F, Field_t & P, - Formulation form); + Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, - Formulation form) = 0; + Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses and tangent moduli, mostly for * debugging and testing. Has runtime-cost associated with * compatibility-checking and conversion of the Field_t arguments that can * be avoided by using the version with strongly typed field references */ void compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, - Formulation form); - + Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); + // this function return the ratio of which the input pixel is consisted of this material + Real get_assigned_ratio(Ccoord pixel); //! iterator to first pixel handled by this material inline iterator begin() {return this->internal_fields.begin();} //! iterator past the last pixel handled by this material inline iterator end() {return this->internal_fields.end();} //! number of pixels assigned to this material inline size_t size() const {return this->internal_fields.size();} protected: const std::string name; //!< material's name (for output and debugging) MFieldCollection_t internal_fields{};//!< storage for internal variables - + using ScalarField_t = ScalarField, Real> ; + ScalarField_t & assigned_ratio; //!< field holding the assigning ratio of the material + using ScalarFieldMap_t = ScalarFieldMap, Real, true>; + ScalarFieldMap_t assigned_ratio_mapped; private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 9fb5bea..3665355 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,718 +1,880 @@ /** * @file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * Copyright © 2017 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 MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H #include "common/common.hh" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.hh" #include #include #include #include namespace muSpectre { // Forward declaration for factory function template class CellBase; /** * material traits are used by `muSpectre::MaterialMuSpectre` to * break the circular dependence created by the curiously recurring * template parameter. These traits must define * - these `muSpectre::FieldMap`s: * - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a * constant second-order `muSpectre::TensorField` * - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a * writable secord-order `muSpectre::TensorField` * - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a * writable fourth-order `muSpectre::TensorField` * - `strain_measure`: the expected strain type (will be replaced by the * small-strain tensor ε * `muspectre::StrainMeasure::Infinitesimal` in small * strain computations) * - `stress_measure`: the measure of the returned stress. Is used by * `muspectre::MaterialMuSpectre` to transform it into * Cauchy stress (`muspectre::StressMeasure::Cauchy`) in - * small-strain computations and into first + * small-strain computations and into first * Piola-Kirchhoff stress `muspectre::StressMeasure::PK1` * in finite-strain computations - * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing + * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing * internal variables */ template struct MaterialMuSpectre_traits { }; template class MaterialMuSpectre; /** * Base class for most convenient implementation of materials */ template class MaterialMuSpectre: public MaterialBase { public: /** * type used to determine whether the - * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only + * `muSpectre::MaterialMuSpectre::iterable_proxy evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; //!< base class //! global field collection using GFieldCollection_t = typename Parent::GFieldCollection_t; //! expected type for stress fields using StressField_t = typename Parent::StressField_t; //! expected type for strain fields using StrainField_t = typename Parent::StrainField_t; //! expected type for tangent stiffness fields using TangentField_t = typename Parent::TangentField_t; - + using Ccoord = Ccoord_t; //!< cell coordinates type //! traits for the CRTP subclass using traits = MaterialMuSpectre_traits; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) = delete; //! Destructor virtual ~MaterialMuSpectre() = default; //! Factory template static Material & make(CellBase & cell, ConstructorArgs &&... args); //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete; //! allocate memory, etc virtual void initialise() override; + // this fucntion is speicalized to assign partial material to a pixel + //void add_pixel_split(const Ccoord & local_ccoord, Real ratio = 1.0); + template + void add_pixel_split(const Ccoord_t & pixel, Real ratio, InternalArgs ... args); + using Parent::compute_stresses; using Parent::compute_stresses_tangent; //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, - Formulation form) override final; + Formulation form, + SplittedCell is_cell_splitted) override; //! computes stress and tangent modulus virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, - Formulation form) override final; + Formulation form, + SplittedCell is_cell_splitted)override; protected: //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 - template + template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P) __attribute__ ((visibility ("default"))); //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 - template + template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K) __attribute__ ((visibility ("default"))); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate over //! that does or does not include tangent moduli template class iterable_proxy; /** * inheriting classes with internal variables need to overload this function */ typename traits::InternalVariables& get_internals() { return static_cast(*this).get_internals();} bool is_initialised{false}; //!< to handle double initialisation right private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre:: MaterialMuSpectre(std::string name) - :Parent(name) { + :Parent(name){ using stress_compatible = typename traits::StressMap_t:: template is_compatible; using strain_compatible = typename traits::StrainMap_t:: template is_compatible; using tangent_compatible = typename traits::TangentMap_t:: template is_compatible; static_assert((stress_compatible::value && stress_compatible::explain()), "The material's declared stress map is not compatible " "with the stress field. More info in previously shown " "assert."); static_assert((strain_compatible::value && strain_compatible::explain()), "The material's declared strain map is not compatible " "with the strain field. More info in previously shown " "assert."); static_assert((tangent_compatible::value && tangent_compatible::explain()), "The material's declared tangent map is not compatible " "with the tangent field. More info in previously shown " "assert."); } /* ---------------------------------------------------------------------- */ template template Material & MaterialMuSpectre:: make(CellBase & cell, ConstructorArgs && ... args) { auto mat = std::make_unique(args...); auto & mat_ref = *mat; cell.add_material(std::move(mat)); return mat_ref; } - + template + template + void MaterialMuSpectre:: + add_pixel_split(const Ccoord &local_ccoord, Real ratio, InternalArgs ...args){ + auto & this_mat = static_cast(*this); + this_mat.add_pixel(local_ccoord, args...); + this->assigned_ratio.push_back(ratio); + } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: initialise() { if (!this->is_initialised) { this->internal_fields.initialise(); this->is_initialised = true; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses(const StrainField_t &F, StressField_t &P, - Formulation form) { + Formulation form, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { - this->template compute_stresses_worker(F, P); + switch (is_cell_splitted){ + case (SplittedCell::no): { + this->compute_stresses_worker(F, P); + break; + } + case (SplittedCell::yes): { + this->compute_stresses_worker(F, P); + break; + } + } break; } case Formulation::small_strain: { - this->template compute_stresses_worker(F, P); + switch (is_cell_splitted){ + case (SplittedCell::no): { + this->compute_stresses_worker(F, P); + break; + } + case (SplittedCell::yes): { + this->compute_stresses_worker(F, P); + break; + } + } break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, - Formulation form) { + Formulation form, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { - this->template compute_stresses_worker(F, P, K); + switch (is_cell_splitted){ + case (SplittedCell::no): { + this->compute_stresses_worker(F, P, K); + break; + } + case (SplittedCell::yes): { + this->compute_stresses_worker(F, P, K); + break; + } + } break; } case Formulation::small_strain: { - this->template compute_stresses_worker(F, P, K); + switch (is_cell_splitted){ + case (SplittedCell::no): { + this->compute_stresses_worker(F, P, K); + break; + } + case (SplittedCell::yes): { + this->compute_stresses_worker(F, P, K); + break; + } + } break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template - template + template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] - (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t Stresses, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ - get_formulation_strain_type(Form, traits::strain_measure)}; - + get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple - auto & F = std::get<0>(Strains); - auto && strain = MatTB::convert_strain(F); + auto & grad = std::get<0>(Strains); + auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli - Stresses = + + // for the case that cells do not contain any splitt cells + auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ + Stresses = + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress_tangent(std::move(strain), + internals...);}, + internal_variables); + }; + + // for the case that cells contain splitt cells + auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ + auto stress_tgt_contributions = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); + auto && stress = std::get<0>(Stresses); + auto && tangent = std::get<1>(Stresses); + stress += ratio * std::get<0>(stress_tgt_contributions); + tangent += ratio * std::get<1>(stress_tgt_contributions); + }; + + switch (is_cell_splitted){ + case SplittedCell::no : { + non_split_pixel(); + break; + } + case SplittedCell::yes : { + split_pixel(); + break; + } + } }; auto constitutive_law_finite_strain = [this] - (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t Stresses, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); - // TODO: Figure this out: I can't std::move(internals...), - // because if there are no internals, compilation fails with "no - // matching function for call to ‘move()’'. These are tuples of - // lvalue references, so it shouldn't be too bad, but still - // irksome. - // return value contains a tuple of rvalue_refs to both stress // and tangent moduli - auto stress_tgt = + + // for the case that cells do not contain splitt cells + auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ + auto stress_tgt = + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress_tangent(std::move(strain), + internals...);}, + internal_variables); + Stresses = MatTB::PK1_stress + (std::move(grad), + std::move(std::get<0>(stress_tgt)), + std::move( std::get<1>(stress_tgt))); + }; + + // for the case that cells contain splitt cells + auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ + auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); - auto & stress = std::get<0>(stress_tgt); - auto & tangent = std::get<1>(stress_tgt); - Stresses = MatTB::PK1_stress - (std::move(grad), - std::move(stress), - std::move(tangent)); + auto && stress_tgt_contributions = MatTB::PK1_stress + (std::move(grad), + std::move(std::get<0>(stress_tgt)), + std::move(std::get<1>(stress_tgt))); + auto && stress = std::get<0>(Stresses); + auto && tangent = std::get<1>(Stresses); + stress += ratio * std::get<0>(stress_tgt_contributions); + tangent += ratio * std::get<1>(stress_tgt_contributions); + }; + + switch (is_cell_splitted){ + case SplittedCell::no : { + non_split_pixel(); + break; + } + case SplittedCell::yes : { + split_pixel(); + break; + } + } }; iterable_proxy fields{*this, F, P, K}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for tangent reference, check iterator" "value_type"); + static_assert(std::is_same(arglist))>>::value, + "Type mismatch for ratio reference, expected a real number"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template - template + template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] - (Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t Stresses,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple - auto & F = std::get<0>(Strains); - auto && strain = MatTB::convert_strain(F); + auto & grad = std::get<0>(Strains); + auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli auto & sigma = std::get<0>(Stresses); + + //for the case that cell does not contain any splitted pixel + auto non_split_pixel = [&strain, &this_mat,ratio, &sigma, &internal_variables](){ sigma = - apply([&strain, &this_mat] (auto && ... internals) { - return - this_mat.evaluate_stress(std::move(strain), - internals...);}, - internal_variables); + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); + }; + + // for the case that cells contain splitt cells + auto split_pixel = [&strain, &this_mat, &ratio, &sigma, &internal_variables](){ + sigma += ratio * + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); + }; + + switch (is_cell_splitted){ + case SplittedCell::no : { + non_split_pixel(); + break; + } + case SplittedCell::yes : { + split_pixel(); + break; + } + } }; auto constitutive_law_finite_strain = [this] - (Strains_t Strains, Stresses_t && Stresses, auto && internal_variables) { + (Strains_t Strains, Stresses_t && Stresses,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple - auto & F = std::get<0>(Strains); - auto && strain = MatTB::convert_strain(F); + auto & grad = std::get<0>(Strains); + auto && strain = MatTB::convert_strain(grad); // TODO: Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. - // return value contains a tuple of rvalue_refs to both stress - // and tangent moduli - auto && stress = - apply([&strain, &this_mat] (auto && ... internals) { - return - this_mat.evaluate_stress(std::move(strain), - internals...);}, - internal_variables); + //for the case that cell does not contain any splitted pixel + auto non_split_pixel = [ &strain, &this_mat, &ratio, &Stresses, &internal_variables, &grad] (){ + auto stress = + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); auto & P = get<0>(Stresses); P = MatTB::PK1_stress - (F, stress); + (grad, stress); + }; + + // for the case that cells contain splitted cells + auto split_pixel = [&strain, &this_mat, ratio, &Stresses, internal_variables, &grad](){ + auto stress = + apply([&strain, &this_mat] (auto && ... internals) { + return + this_mat.evaluate_stress(std::move(strain), + internals...);}, + internal_variables); + + auto && P = get<0>(Stresses); + P += ratio * MatTB::PK1_stress + (grad, stress); + }; + + switch (is_cell_splitted){ + case SplittedCell::no : { + non_split_pixel(); + break; + } + case SplittedCell::yes : { + split_pixel(); + break; + } + } }; iterable_proxy fields{*this, F, P}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue - * references (see value_tye in the class definition of + * references (see value_type in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); + static_assert(std::is_same(arglist))>>::value, + "Type mismatch for ratio reference, expected a real number"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K) :material{mat}, strain_field{F}, stress_tup{P,K}, internals{material.get_internals()}{}; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) :material{mat}, strain_field{F}, stress_tup{P}, internals{material.get_internals()}{}; //! Expected type for strain fields using StrainMap_t = typename traits::StrainMap_t; //! Expected type for stress fields using StressMap_t = typename traits::StressMap_t; //! Expected type for tangent stiffness fields using TangentMap_t = typename traits::TangentMap_t; //! expected type for strain values using Strain_t = typename traits::StrainMap_t::reference; //! expected type for stress values using Stress_t = typename traits::StressMap_t::reference; //! expected type for tangent stiffness values using Tangent_t = typename traits::TangentMap_t::reference; //! tuple of intervnal variables, depends on the material using InternalVariables = typename traits::InternalVariables; //! tuple containing a stress and possibly a tangent stiffness field using StressFieldTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness field map using StressMapTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness value ref using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! Copy constructor iterable_proxy(const iterable_proxy &other) = default; //! Move constructor iterable_proxy(iterable_proxy &&other) = default; //! Destructor virtual ~iterable_proxy() = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) = default; /** * dereferences into a tuple containing strains, and internal * variables, as well as maps to the stress and potentially * stiffness maps where to write the response of a pixel */ class iterator { public: //! type to refer to internal variables owned by a CRTP material using InternalReferences = MatTB::ReferenceTuple_t; //! return type to be unpacked per pixel my the constitutive law using value_type = - std::tuple, Stress_tTup, InternalReferences>; + std::tuple, Stress_tTup, Real, InternalReferences>; using iterator_category = std::forward_iterator_tag; //!< stl conformance //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterator(const iterable_proxy & it, bool begin = true) : it{it}, strain_map{it.strain_field}, stress_map {it.stress_tup}, index{begin ? 0:it.material.internal_fields.size()}{} //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const iterable_proxy & it; //!< ref to the proxy StrainMap_t strain_map; //!< map onto the global strain field //! map onto the global stress field and possibly tangent stiffness StressMapTup stress_map; size_t index; //!< index or pixel currently referred to private: }; //! returns iterator to first pixel if this material iterator begin() {return std::move(iterator(*this));} //! returns iterator past the last pixel in this material iterator end() {return std::move(iterator(*this, false));} protected: MaterialMuSpectre & material; //!< reference to the proxied material const StrainField_t & strain_field; //!< cell's global strain field //! references to the global stress field and perhaps tangent stiffness StressFieldTup stress_tup; //! references to the internal variables InternalVariables & internals; private: }; /* ---------------------------------------------------------------------- */ template template bool MaterialMuSpectre::iterable_proxy::iterator:: operator!=(const iterator & other) const { return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy:: iterator & MaterialMuSpectre::iterable_proxy::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy::iterator:: value_type MaterialMuSpectre::iterable_proxy::iterator:: operator*() { const Ccoord_t pixel{ this->it.material.internal_fields.get_ccoord(this->index)}; auto && strain = std::make_tuple(this->strain_map[pixel]); - + auto && ratio = this->it.material.get_assigned_ratio(pixel); auto && stresses = apply([&pixel] (auto && ... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...);}, this->stress_map); auto && internal = this->it.material.get_internals(); const auto id{this->index}; auto && internals = apply([id] (auto && ... internals_) { return InternalReferences{internals_[id]...};}, internal); return std::make_tuple(std::move(strain), std::move(stresses), + std::move(ratio), std::move(internals)); } - } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/tests/test_patch_split_octree_cell_newton_cg_solver.cc b/tests/test_patch_split_octree_cell_newton_cg_solver.cc new file mode 100644 index 0000000..6db1aa6 --- /dev/null +++ b/tests/test_patch_split_octree_cell_newton_cg_solver.cc @@ -0,0 +1,190 @@ +/** + * @file test_patch_split_octree_cell_newton_cg_solver.cc + * + * @author Ali Faslfi + * + * @date 12 Jun 2018 + * + * @brief Tests for split cells and octree material assignment + * + * Copyright © 2017 Till Ali Faslafi + * + * µ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 "tests.hh" +//#include "solver/solvers.hh" +//#include "solver/solver_cg.hh" +//#include "solver/solver_eigen.hh" +#include "solver/deprecated_solvers.hh" +#include "solver/deprecated_solver_cg.hh" +#include "solver/deprecated_solver_cg_eigen.hh" +#include "fft/fftw_engine.hh" +#include "fft/projection_finite_strain_fast.hh" +#include "materials/material_linear_elastic1.hh" +#include "common/iterators.hh" +#include "common/ccoord_operations.hh" +#include "common/common.hh" +#include "cell/cell_factory.hh" +#include "cell/cell_split.hh" +#include "common/intersection_octree.hh" + +#include +#include + +namespace muSpectre { + + + BOOST_AUTO_TEST_SUITE(split_cell_octree_tests); + + BOOST_AUTO_TEST_CASE(octree_intersection) { + + constexpr Dim_t dim{twoD}; + + using Rcoord = Rcoord_t; + using Ccoord = Ccoord_t; + + using Mat_t = MaterialLinearElastic1; + const Real contrast {10}; + const Real Young_soft{1.0030648180242636}, Poisson_soft{0.29930675909878679}; + const Real Young_hard{contrast * Young_soft}, Poisson_hard{0.29930675909878679}; + const Real machine_precision = 1e-14; + + constexpr Ccoord resolutions_split{11,11}; + constexpr Rcoord lengths_split{1.1, 1.1}; + auto fft_ptr_split_1{std::make_unique>(resolutions_split, ipow(dim, 2))}; + auto proj_ptr_split_1{std::make_unique>(std::move(fft_ptr_split_1), lengths_split)}; + + auto fft_ptr_split_2{std::make_unique>(resolutions_split, ipow(dim, 2))}; + auto proj_ptr_split_2{std::make_unique>(std::move(fft_ptr_split_2), lengths_split)}; + + CellSplit sys_split_1(std::move(proj_ptr_split_1), SplittedCell::yes); + CellSplit sys_split_2(std::move(proj_ptr_split_2), SplittedCell::yes); + + auto& Material_hard_split_1 = Mat_t::make(sys_split_1, "hard", 10*Young_hard, Poisson_hard); + auto& Material_soft_split_1 = Mat_t::make(sys_split_1, "soft", Young_soft, Poisson_soft); + + auto& Material_hard_split_2 = Mat_t::make(sys_split_2, "hard", 10*Young_hard, Poisson_hard); + auto& Material_soft_split_2 = Mat_t::make(sys_split_2, "soft", Young_soft, Poisson_soft); + + //square precipitate vertices: + constexpr Rcoord center{0.55, 0.55}; + constexpr Rcoord halfside{0.11, 0.11}; + constexpr Real halfdiagon{sqrt( pow(halfside[0], 2) + + pow(halfside[1], 2) )}; + + //making the first precipitate + std::vector precipitate_1_vertices ; + precipitate_1_vertices.push_back ({center[0] - halfside[0], center[1] - halfside[1]}); + precipitate_1_vertices.push_back ({center[0] + halfside[0], center[1] - halfside[1]}); + precipitate_1_vertices.push_back ({center[0] - halfside[0], center[1] + halfside[1]}); + precipitate_1_vertices.push_back ({center[0] + halfside[0], center[1] + halfside[1]}); + + RootNode precipitate_1 (sys_split_1, precipitate_1_vertices); + + //Extracting the intersected pixels and their correspondent intersection ratios: + auto precipitate_1_intersects = precipitate_1.get_intersected_pixels(); + auto precipitate_1_intersection_ratios = precipitate_1.get_intersection_ratios(); + + //making the second precipitate + std::vector precipitate_2_vertices ; + precipitate_2_vertices.push_back ({center[0] - halfdiagon, center[1] }); + precipitate_2_vertices.push_back ({center[0] + halfdiagon, center[1] }); + precipitate_2_vertices.push_back ({center[0] , center[1] - halfdiagon}); + precipitate_2_vertices.push_back ({center[0] , center[1] + halfdiagon}); + + RootNode precipitate_2 (sys_split_2, precipitate_2_vertices); + auto precipitate_2_intersects = precipitate_2.get_intersected_pixels(); + auto precipitate_2_intersection_ratios = precipitate_2.get_intersection_ratios(); + + //assign precipitate materials: + typedef std::vector::iterator Ccoord_iter_1, Ccoord_iter_2; + typedef std::vector::iterator Real_iter_1, Real_iter_2; + + for (std::pair it(precipitate_1_intersects.begin(), precipitate_1_intersection_ratios.begin()) ; + it.first != precipitate_1_intersects.end(); ++it.first , ++it.second ){ + Material_hard_split_1.add_pixel_split(*it.first, *it.second); + Material_soft_split_1.add_pixel_split(*it.first, 1.0-*it.second); + } + + for (std::pair it(precipitate_2_intersects.begin(), precipitate_2_intersection_ratios.begin()) ; + it.first != precipitate_2_intersects.end(); ++it.first , ++it.second ){ + Material_hard_split_2.add_pixel_split(*it.first, *it.second); + Material_soft_split_2.add_pixel_split(*it.first, 1.0-*it.second); + } + + //complete material assignment + //sys_split_1.complete_material_assignment(Material_soft_split); + std::vector assigned_ratio_1 = sys_split_1.get_intersected_ratios(); + int iterator = 0 ; + for (auto && tup: akantu::enumerate(sys_split_1)) { + auto && pixel = std::get<1>(tup); + if (assigned_ratio_1[iterator] != 1.0){ + Material_soft_split_1.add_pixel_split(pixel, 1.0-assigned_ratio_1[iterator]); + } + iterator++; + } + + std::vector assigned_ratio_2 = sys_split_2.get_intersected_ratios(); + iterator = 0 ; + for (auto && tup: akantu::enumerate(sys_split_2)) { + auto && pixel = std::get<1>(tup); + if (assigned_ratio_2[iterator] != 1.0){ + Material_soft_split_2.add_pixel_split(pixel, 1.0-assigned_ratio_2[iterator]); + } + iterator++; + } + //checking whether the octree and intersected volume worked correctly: + Real area_preticipitate_1 = 0.0 ,area_preticipitate_2 = 0.0 ; + for (auto&& precepetiate_area:precipitate_1_intersection_ratios){ + area_preticipitate_1 += precepetiate_area; + } + for (auto&& precepetiate_area:precipitate_2_intersection_ratios){ + area_preticipitate_2 += precepetiate_area; + } + BOOST_CHECK_LE(abs(area_preticipitate_2 - area_preticipitate_1), area_preticipitate_1 * machine_precision); + ////////////// + + //Finding the equilbrium with Sppectral method: + sys_split_1.initialise(); + sys_split_2.initialise(); + + Eigen::Rotation2D rotation_1_2( pi/4); + Eigen::Rotation2D rotation_2_1(-pi/4); + Grad_t delF0_1, delF0_2; + delF0_1 << 1.0, 0.0, 0.0, 0.0; + delF0_2 = rotation_1_2 * delF0_1; + + + constexpr Real cg_tol_1{1e-8}, newton_tol_1{1e-5}; + constexpr Uint maxiter_1{CcoordOps::get_size(resolutions_split)*ipow(dim, secondOrder)*10}; + constexpr bool verbose_1{false}; + + GradIncrements grads_1; grads_1.push_back(delF0_1); + + DeprecatedSolverCG cg_1{sys_split_1, cg_tol_1, maxiter_1, bool(verbose_1)}; + Eigen::ArrayXXd res1{deprecated_de_geus(sys_split_1, grads_1, cg_1, newton_tol_1, verbose_1)[0].grad}; + + DeprecatedSolverCGEigen cg_2{sys_split_2, cg_tol_1, maxiter_1, bool(verbose_1)}; + Eigen::ArrayXXd res2{deprecated_de_geus(sys_split_2, grads_1, cg_2, newton_tol_1, verbose_1)[0].grad}; + + + BOOST_CHECK_LE(abs(0), cg_tol_1); + } + + BOOST_AUTO_TEST_SUITE_END(); + +} // muSpectre