diff --git a/CMakeLists.txt b/CMakeLists.txt index a8b28c0..f399f9b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,58 +1,66 @@ cmake_minimum_required(VERSION 3.7.0) project(µSpectre_prototype) set(CMAKE_CXX_STANDARD 14) set(CMAKE_EXPORT_COMPILE_COMMANDS ON) add_compile_options(-Wall -Wextra -Werror) enable_testing() set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") # using Clang elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") # using GCC string( TOLOWER "${CMAKE_BUILD_TYPE}" build_type ) if ("release" STREQUAL "${build_type}" ) add_compile_options(-march=native) add_compile_options(-mavx) add_compile_options(-mfma) add_compile_options(-msse) 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() find_package(Boost COMPONENTS unit_test_framework REQUIRED) find_package(Eigen3 3.3 REQUIRED NO_MODULE) include_directories( tests src ${EIGEN3_INCLUDE_DIRS} ) #find_package(MPI REQUIRED) find_package(FFTW REQUIRED) #find_package(FFTWMPI REQUIRED) add_subdirectory( "${CMAKE_SOURCE_DIR}/src/" ) #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} Eigen3::Eigen materials common systems) add_test(main_test_suite main_test_suite) +file( GLOB binaries "bin/*.cc") +foreach(binaryfile ${binaries}) + get_filename_component(binaryname ${binaryfile} NAME_WE) + add_executable(${binaryname} ${binaryfile}) + target_link_libraries(${binaryname} ${Boost_LIBRARIES} Eigen3::Eigen materials common systems) +endforeach(binaryfile ${binarys}) + + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 851311e..f8ad1e0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,10 +1,11 @@ include_directories( common materials system + external ) add_subdirectory(common) add_subdirectory(materials) add_subdirectory(system) diff --git a/src/materials/material_hyper_elastic.cc b/src/materials/material_hyper_elastic.cc index 3d35fd9..2f03131 100644 --- a/src/materials/material_hyper_elastic.cc +++ b/src/materials/material_hyper_elastic.cc @@ -1,126 +1,126 @@ /** * file material_hyper_elastic.cc * * @author Till Junge * * @date 01 May 2017 * * @brief implementation for hyperelastic material law * * @section LICENCE * * Copyright (C) 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 "material_hyper_elastic.hh" #include "materials/materials_toolbox.hh" #include namespace muSpectre { //----------------------------------------------------------------------------// template MaterialHyperElastic::MaterialHyperElastic(Real Young, Real Poisson) : Young(Young), Poisson(Poisson), lambda(Young*Poisson/(1+Poisson)/(1-2*Poisson)), mu(Young/2/(1+Poisson)), C(lambda*Tensors::outer(Tensors::I2(),Tensors::I2()) + 2*mu*Tensors::I4S()) { } //----------------------------------------------------------------------------// template void MaterialHyperElastic::initialize() { auto&& nb_pixel = this->indices.size(); this->gradF = Eigen::Tensor(m_dim, m_dim, nb_pixel); this->stressS = Eigen::Tensor(m_dim, m_dim, nb_pixel); this->stressP = Eigen::Tensor(m_dim, m_dim, nb_pixel); this->stiffnessK = Eigen::Tensor(m_dim, m_dim, m_dim, m_dim, nb_pixel); } //----------------------------------------------------------------------------// template void MaterialHyperElastic:: compute_First_Piola_Kirchhoff_stress(const SecondArrayMap & grad_array, SecondArrayMap & stress_array, FourthArrayMap & stiffness_array) { auto&& nb_pixel = this->indices.size(); // make local contiguous copy of gradient and compute Green-Lagrange strain for (size_t i = 0; i < nb_pixel; ++i) { auto F_input {this->get_const_view(i, grad_array)}; SecondMap F_copy(&this->gradF(0, 0, i)); F_copy = F_input; } // do the math in a single loop for (size_t i = 0; i < nb_pixel; ++i) { SecondMap F(&this->gradF(0, 0, i)); SecondMap S(&this->stressS(0, 0, i)); SecondMap P(&this->stressP(0, 0, i)); SecondTens E = .5*(F.transpose()*F-SecondMap::Identity()); using fix_siz = Eigen::Sizes; using fix_ten = Eigen::TensorFixedSize; using fix_map = Eigen::TensorMap; fix_map E_t(&E(0, 0), m_dim, m_dim); fix_map S_t(&S(0, 0), m_dim, m_dim); std::array, 2> prod_dims {Eigen::IndexPair{2, 1}, Eigen::IndexPair{3, 0}}; S_t = this->C.contract(E_t, prod_dims); // First Piola-Kirchhoff stress P = F*S; // product dims for dot42 with transposed second argument std::array, 1> dims{Eigen::IndexPair{3, 1}}; fix_map F_t(&F(0, 0), m_dim, m_dim); FourthMap K(&this->stiffnessK(0, 0, 0, 0, i), m_dim, m_dim, m_dim, m_dim); auto C_FT = this->C.contract(F_t, dims); auto F_C_FT = Tensors::dot24(F_t, C_FT); auto IRT_F_C_FT = Tensors::dotdot(Tensors::I4RT(), F_C_FT); auto IRT_F_C_FT_IRT = Tensors::dotdot< m_dim>(IRT_F_C_FT, Tensors::I4RT()); auto S_I = Tensors::dot24(S_t, Tensors::I4()); - K = S_I+ IRT_F_C_FT_IRT; + K = (S_I+ IRT_F_C_FT_IRT).shuffle(std::array{1,0,2,3}); } // make local contiguous copy of gradient and compute Green-Lagrange strain for (size_t i = 0; i < nb_pixel; ++i) { auto P_output {this->get_view(i, stress_array)}; SecondMap P(&this->stressP(0, 0, i)); P_output = P; auto K_output {this->get_view(i, stiffness_array)}; FourthMap K(&this->stiffnessK(0, 0, 0, 0, i), m_dim, m_dim, m_dim, m_dim); K_output = K; } } template class MaterialHyperElastic<2, 2>; template class MaterialHyperElastic<2, 3>; template class MaterialHyperElastic<3, 3>; } // muSpectre diff --git a/src/system/fft_engine.hh b/src/system/fft_engine.hh index 434ab83..23ca526 100644 --- a/src/system/fft_engine.hh +++ b/src/system/fft_engine.hh @@ -1,65 +1,65 @@ /** * file fft_engine.hh * * @author Till Junge * * @date 11 May 2017 * * @brief abstract base class defining interface for fft_systems * * @section LICENCE * * Copyright (C) 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 "common/common.hh" #include #include #ifndef FFT_ENGINE_H #define FFT_ENGINE_H namespace muSpectre { enum class FFT_PlanFlags {estimate, measure, patient}; template class FFT_Engine { public: using Index = Index_t; FFT_Engine(Index nb_pixels) :nb_pixels(nb_pixels), tot_nb_pixels(std::accumulate(nb_pixels.begin(), nb_pixels.end(), 1, std::multiplies())){}; virtual Real * convolve(const Real * const Ghat, const Real * const arr) = 0; //! out-of place convolution virtual void convolve (const Real * const Ghat, const Real * const in, Real * out) = 0; - + size_t get_tot_nb_pixels() const {return this->tot_nb_pixels;} protected: virtual void init(FFT_PlanFlags flags) = 0; const Index nb_pixels; size_t tot_nb_pixels; }; } // muSpectre #endif /* FFT_ENGINE_H */ diff --git a/src/system/fftw_engine_c2c.cc b/src/system/fftw_engine_c2c.cc index 41a670c..925001f 100644 --- a/src/system/fftw_engine_c2c.cc +++ b/src/system/fftw_engine_c2c.cc @@ -1,162 +1,212 @@ /** * file fftw_engine_c2c.cc * * @author Till Junge * * @date 11 May 2017 * * @brief fft_engine usinge FFTW * * @section LICENCE * * Copyright (C) 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 "system/fftw_engine_c2c.hh" #include #include namespace muSpectre { template FFTW_EngineC2C::FFTW_EngineC2C(Index nb_pixels, FFT_PlanFlags flags) :FFT_Engine(nb_pixels) { this->i_work_space = fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM); this->o_work_space = fftw_alloc_complex(this->tot_nb_pixels* DimM*DimM); this->init(flags); } //----------------------------------------------------------------------------// template FFTW_EngineC2C::~FFTW_EngineC2C() { fftw_destroy_plan(this->plan_fft); fftw_destroy_plan(this->plan_ifft); fftw_free(this->i_work_space); fftw_free(this->o_work_space); fftw_cleanup(); } //----------------------------------------------------------------------------// template void FFTW_EngineC2C::init(FFT_PlanFlags inflags) { const int & rank = DimS; std::array narr; const int * const n = &narr[0]; std::copy(this->nb_pixels.begin(), this->nb_pixels.end(), narr.begin()); int howmany = DimM*DimM; - fftw_complex * in = reinterpret_cast(this->i_work_space); - const int * const inembed = n; + fftw_complex * in = this->i_work_space; + const int * const inembed = nullptr; int istride = howmany; int idist = 1; fftw_complex * out = in; - const int * const onembed = n; + const int * const onembed = nullptr; int ostride = howmany; int odist = idist; int sign = FFTW_FORWARD; unsigned int flags; if (inflags == FFT_PlanFlags::estimate) { flags = FFTW_ESTIMATE; } if (inflags == FFT_PlanFlags::measure) { flags = FFTW_MEASURE; } if (inflags == FFT_PlanFlags::patient) { flags = FFTW_PATIENT; } this->plan_fft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, sign, flags); + if (this->plan_fft == nullptr) { + throw std::runtime_error("Plan failed"); + } sign = FFTW_BACKWARD; - in = reinterpret_cast(this->o_work_space); - out = reinterpret_cast(this->o_work_space); + in = this->o_work_space; + out = in; this->plan_ifft = fftw_plan_many_dft(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, sign, flags); + if (this->plan_ifft == nullptr) { + throw std::runtime_error("Plan failed"); + } } //----------------------------------------------------------------------------// template Real * FFTW_EngineC2C::convolve(const Real * const Ghat, const Real * const arr) { std::copy(arr, arr+this->tot_nb_pixels*DimM*DimM, reinterpret_cast(this->i_work_space)); fftw_execute(this->plan_fft); //----------------------------------------------------------------------------// const std::array,2> prod_dims {Eigen::IndexPair{2, 1}, Eigen::IndexPair{3, 0}}; //set output to zero Eigen::Map> o_space_handle(reinterpret_cast(this->o_work_space), this->tot_nb_pixels* DimM*DimM) ; o_space_handle.setZero(); for (size_t pix_id = 0; pix_id < this->tot_nb_pixels; ++pix_id) { using T4shape = Eigen::Sizes; using T2shape = Eigen::Sizes; using T4type = Eigen::TensorFixedSize; using T2type = Eigen::TensorFixedSize; using T4map = Eigen::TensorMap; using T2map = Eigen::TensorMap; const Real* const G_ptr = &Ghat[pix_id*DimM*DimM*DimM*DimM]; const T4map G(const_cast(G_ptr),DimM, DimM, DimM,DimM); Complex * P_ptr = reinterpret_cast(&this->i_work_space[pix_id*DimM*DimM]); const T2map P(P_ptr, DimM, DimM); Complex* Pout_ptr = reinterpret_cast(&this->o_work_space[pix_id*DimM*DimM]); T2map Pout(Pout_ptr, DimM, DimM); - ////TODO manual contraction to work around complex- - //for (Dim_t i = 0; i < DimM; ++i) { - // for (Dim_t j = 0; j < DimM; ++j) { - // for (Dim_t k = 0; k < DimM; ++k) { - // for (Dim_t l = 0; l < DimM; ++l) { - // Pout(i,j) = G(i,j,k,l)*P(l,k); - // } - // } - // } - //} Eigen::TensorFixedSize Gc = G.template cast(); Pout = Gc.contract(P, prod_dims); } //----------------------------------------------------------------------------// fftw_execute(this->plan_ifft); + Real * ret_arr = reinterpret_cast(this->o_work_space); + for (size_t pix_id = 0; pix_id < this->tot_nb_pixels; ++pix_id) { + ret_arr[pix_id] = this->o_work_space[pix_id][0]/this->tot_nb_pixels; + } + return ret_arr; + } + //----------------------------------------------------------------------------// + template + void FFTW_EngineC2C::convolve(const Real * const Ghat, + const Real * const in, + Real * out) { + + std::copy(in, in+this->tot_nb_pixels*DimM*DimM, + reinterpret_cast(this->i_work_space)); + + fftw_execute(this->plan_fft); + + //----------------------------------------------------------------------------// + const std::array,2> + prod_dims {Eigen::IndexPair{2, 1}, + Eigen::IndexPair{3, 0}}; + + //set output to zero + Eigen::Map> + o_space_handle(reinterpret_cast(this->o_work_space), + this->tot_nb_pixels* DimM*DimM) ; + o_space_handle.setZero(); + + for (size_t pix_id = 0; pix_id < this->tot_nb_pixels; ++pix_id) { - return reinterpret_cast(this->o_work_space); + using T4shape = Eigen::Sizes; + using T2shape = Eigen::Sizes; + using T4type = Eigen::TensorFixedSize; + using T2type = Eigen::TensorFixedSize; + using T4map = Eigen::TensorMap; + using T2map = Eigen::TensorMap; + const Real* const G_ptr = &Ghat[pix_id*DimM*DimM*DimM*DimM]; + const T4map G(const_cast(G_ptr),DimM, DimM, DimM,DimM); + Complex * P_ptr = + reinterpret_cast(&this->i_work_space[pix_id*DimM*DimM]); + const T2map P(P_ptr, DimM, DimM); + + Complex* Pout_ptr = + reinterpret_cast(&this->o_work_space[pix_id*DimM*DimM]); + T2map Pout(Pout_ptr, DimM, DimM); + + Eigen::TensorFixedSize Gc = G.template cast(); + Pout = Gc.contract(P, prod_dims)/static_cast(this->tot_nb_pixels); + } + //----------------------------------------------------------------------------// + fftw_execute(this->plan_ifft); + for (size_t pix_id = 0; pix_id < this->tot_nb_pixels*DimM*DimM; ++pix_id) { + out[pix_id] = this->o_work_space[pix_id][0]; + } } template class FFTW_EngineC2C<1, 1>; + template class FFTW_EngineC2C<2, 1>; template class FFTW_EngineC2C<2, 2>; template class FFTW_EngineC2C<3, 3>; } // muSpectre diff --git a/src/system/fftw_engine_c2c.hh b/src/system/fftw_engine_c2c.hh index f98c3a1..0adaa96 100644 --- a/src/system/fftw_engine_c2c.hh +++ b/src/system/fftw_engine_c2c.hh @@ -1,59 +1,59 @@ /** * file fftw_engine_c2c.hh * * @author Till Junge * * @date 11 May 2017 * * @brief implementaion of fft_engine using redundant c2c fftw transform * * @section LICENCE * * Copyright (C) 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 "system/fft_engine.hh" #include #ifndef FFTW_ENGINE_C2C_H #define FFTW_ENGINE_C2C_H namespace muSpectre { template class FFTW_EngineC2C: public FFT_Engine { public: using Index = Index_t; FFTW_EngineC2C(Index nb_pixels, FFT_PlanFlags flags); virtual ~FFTW_EngineC2C(); virtual Real * convolve(const Real * const Ghat, const Real * const arr) override; - virtual void convolve (const Real * const Ghat, const Real * const in, Real * out) override {if(Ghat == in || out) {}} + virtual void convolve (const Real * const Ghat, const Real * const in, Real * out) override; protected: virtual void init(FFT_PlanFlags flags) override; fftw_complex * i_work_space, * o_work_space; fftw_plan plan_fft; fftw_plan plan_ifft; }; } // muSpectre #endif /* FFTW_ENGINE_C2C_H */ diff --git a/src/system/fftw_engine_r2c.cc b/src/system/fftw_engine_r2c.cc index d605de1..45cee5c 100644 --- a/src/system/fftw_engine_r2c.cc +++ b/src/system/fftw_engine_r2c.cc @@ -1,192 +1,193 @@ /** * file fftw_engine_r2c.cc * * @author Till Junge * * @date 13 May 2017 * * @brief fft_engine using real-data fftw transform * * @section LICENCE * * Copyright (C) 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 "system/fftw_engine_r2c.hh" #include namespace muSpectre { /** here, I'm doing everything out-of plas, so that the real input arrays * can be contiguous */ template FFTW_EngineR2C::FFTW_EngineR2C(Index nb_pixels, FFT_PlanFlags flags) :FFT_Engine(nb_pixels), nb_hermitian(this->tot_nb_pixels/nb_pixels.back() * (nb_pixels.back()/2+1)) { this->work_space = fftw_alloc_complex(this->nb_hermitian* DimM*DimM); this->init(flags); } //----------------------------------------------------------------------------// template FFTW_EngineR2C::~FFTW_EngineR2C() { fftw_destroy_plan(this->plan_fft); fftw_destroy_plan(this->plan_ifft); fftw_free(this->work_space); fftw_cleanup(); } //----------------------------------------------------------------------------// template void FFTW_EngineR2C::init(FFT_PlanFlags inflags) { const int & rank = DimS; std::array narr; const int * const n = &narr[0]; std::copy(this->nb_pixels.begin(), this->nb_pixels.end(), narr.begin()); int howmany = DimM*DimM; //temporary buffer for plan Real * r_work_space = fftw_alloc_real(this->tot_nb_pixels*howmany); Real * in = r_work_space; const int * const inembed = nullptr;//nembed are tricky: they refer to physical layout int istride = howmany; int idist = 1; fftw_complex * out = this->work_space; const int * const onembed = nullptr; int ostride = howmany; int odist = idist; unsigned int flags; if (inflags == FFT_PlanFlags::estimate) { flags = FFTW_ESTIMATE; } if (inflags == FFT_PlanFlags::measure) { flags = FFTW_MEASURE; } if (inflags == FFT_PlanFlags::patient) { flags = FFTW_PATIENT; } this->plan_fft = fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, FFTW_PRESERVE_INPUT | flags); if (this->plan_fft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_complex * i_in = this->work_space; Real * i_out = r_work_space; this->plan_ifft = fftw_plan_many_dft_c2r(rank, n, howmany, i_in, inembed, istride, idist, i_out, onembed, ostride, odist, flags); if (this->plan_ifft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_free(r_work_space); } //----------------------------------------------------------------------------// template FFTW_EngineR2C::index_iterator::index_iterator(const FFTW_EngineR2C & parent, bool begin) :parent(parent), index(begin ? 0: parent.nb_hermitian) { this->half_strides.back() = this->full_strides.back() = 1; this->full_strides[DimS-2] = this->parent.nb_pixels[DimS-1]; this->half_strides[DimS-2]= this->full_strides[DimS-2]/2+1; for (Dim_t i = DimS-3; i >= 0; --i) { //the full strides are cumulative this->full_strides[i] = this->parent.nb_pixels[i+1] * this->full_strides[i+1]; // the half strides are not (operator* divided tmp_index successively) this->half_strides[i] = this->parent.nb_pixels[i+1]; } } //----------------------------------------------------------------------------// template void FFTW_EngineR2C::convolve(const Real *const Ghat, const Real * const in, Real * out) { fftw_execute_dft_r2c(this->plan_fft, const_cast(in), this->work_space); //----------------------------------------------------------------------------// const std::array,2> prod_dims {Eigen::IndexPair{2, 1}, Eigen::IndexPair{3, 0}}; for (auto pix_id:*this) { const size_t & real_pix_id = std::get<0>(pix_id); const size_t & complex_pix_id = std::get<1>(pix_id); using T4shape = Eigen::Sizes; using T2shape = Eigen::Sizes; using T4type = Eigen::TensorFixedSize; using T2type = Eigen::TensorFixedSize; using T4map = Eigen::TensorMap; using T2map = Eigen::TensorMap; const Real* const G_ptr = &Ghat[real_pix_id*DimM*DimM*DimM*DimM]; const T4map G(const_cast(G_ptr),DimM, DimM, DimM,DimM); Complex * P_ptr = reinterpret_cast(&this->work_space[complex_pix_id*DimM*DimM]); T2map P(P_ptr, DimM, DimM); P = 1./static_cast(this->tot_nb_pixels)*G.template cast().contract(P, prod_dims).eval(); - + } fftw_execute_dft_c2r(this->plan_ifft, this->work_space, out); } //----------------------------------------------------------------------------// template typename FFTW_EngineR2C::index_iterator FFTW_EngineR2C::begin() { return index_iterator(*this); } //----------------------------------------------------------------------------// template typename FFTW_EngineR2C::index_iterator FFTW_EngineR2C::end() { return index_iterator(*this, false); } template class FFTW_EngineR2C<1, 1>; + template class FFTW_EngineR2C<2, 1>; template class FFTW_EngineR2C<2, 2>; template class FFTW_EngineR2C<3, 3>; } // muSpectre diff --git a/src/system/system_base.cc b/src/system/system_base.cc index 0474b90..91d8307 100644 --- a/src/system/system_base.cc +++ b/src/system/system_base.cc @@ -1,379 +1,389 @@ /** * file system_base.cc * * @author Till Junge * * @date 09 May 2017 * * @brief Implementation for system base class * * @section LICENCE * * Copyright (C) 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 "system/system_base.hh" #include "iostream" #include "materials/materials_toolbox.hh" #include "common/common.hh" #include namespace muSpectre { namespace Projection { //----------------------------------------------------------------------------// //! computes fft frequecies. see function fftfreq in numpy's helper.py vector fft_freqs(Uint nb_samples) { vector retval(1, nb_samples); int N = (nb_samples-1)/2 + 1; // needs to be signed int for neg freqs for (int i = 0; i < N; ++i) { retval(0, i) = i; } for (int i = N; i < static_cast(nb_samples); ++i) { retval(0, i) = -(static_cast(nb_samples)/2)+i-N; } return retval; } //----------------------------------------------------------------------------// vector fft_freqs(Uint nb_samples, Real length) { return fft_freqs(nb_samples)/length; } //----------------------------------------------------------------------------// template FreqStruc::FreqStruc(const std::array& sizes, const Index_t& nb_pixels) :nb_pixels(nb_pixels), sizes(sizes) { for (size_t i = 0; i < Dim; ++i) { this->xi[i] = Projection::fft_freqs(this->nb_pixels[i]); } } //----------------------------------------------------------------------------// template inline typename FreqStruc::iterator FreqStruc::begin() { return iterator(*this, {0}); } //----------------------------------------------------------------------------// template inline typename FreqStruc::iterator FreqStruc::end() { return iterator(*this, this->nb_pixels); } //----------------------------------------------------------------------------// template FreqStruc::iterator::iterator(const FreqStruc & parent, Index_t index) :parent(parent), index(std::accumulate(index.begin(), index.end(), 1, std::multiplies())){ } //----------------------------------------------------------------------------// template inline Index_t FreqStruc::iterator::get_index() { // implements row-major choice. might have to be templated at some point static_assert(((Dim == twoD) || (Dim == threeD)), "only 2d or 3d"); auto & stride = this->parent.nb_pixels[1]; return Index_t{static_cast(this->index/stride), static_cast(this->index%stride)}; } //----------------------------------------------------------------------------// template<> inline Index_t FreqStruc::iterator::get_index() { // implements row-major choice. might have to be templated at some point auto && kstride = this->parent.nb_pixels[2]; Dim_t && k = this->index%kstride; auto && jstride = this->parent.nb_pixels[1]; Dim_t && j = (this->index/kstride)%jstride; Dim_t && i = (this->index/kstride)/jstride; return Index_t{i, j, k}; } //----------------------------------------------------------------------------// template inline Eigen::Matrix FreqStruc::iterator::operator*() { Eigen::Matrix retval; auto && indices = this->get_index(); if (Dim == twoD) { retval << parent.xi[0](indices[0]), parent.xi[1](indices[1]); } else if (Dim == threeD) { retval << parent.xi[0](indices[0]), parent.xi[1](indices[1]), parent.xi[2](indices[2]); } return retval; } //----------------------------------------------------------------------------// template inline typename FreqStruc::iterator& FreqStruc::iterator::operator++() { this->index++; return *this; } } // Projection //----------------------------------------------------------------------------// template SystemBase::SystemBase(std::array sizes, Index nb_pixels, Formulation form) : sizes{sizes}, nb_pixels{nb_pixels}, nb_pixel(std::accumulate(nb_pixels.begin(), nb_pixels.end(), 1, std::multiplies())), form{form} { for (const auto & size: this->nb_pixels) { if (size%2 != 1) { throw std::runtime_error("Sorry, for now can only handle odd grid sizes"); } } if (form != Formulation::finite_strain) { throw std::runtime_error("Sorry, only finite strain for the time being"); } //resize and allocate this->Ghats = T4Vector(DimM, DimM, DimM, DimM, this->size()); this->Ks = T4Vector(DimM, DimM, DimM, DimM, this->size()); this->grad = T2Vector(DimM, DimM, this->size()); this->stress = T2Vector(DimM, DimM, this->size()); this->Ghats.setZero(); this->build_ghats(); } //----------------------------------------------------------------------------// template void SystemBase::set_fft_engine(FFT_Ptr fft_eng) { this->fft_engine = fft_eng; } //----------------------------------------------------------------------------// template void SystemBase::solve(const T2 & DeltaF, const Real newton_tol, const Real cg_tol, - Uint maxiter) { + Uint maxiter, bool verbose, bool cg_verbose) { if(maxiter == 0) { maxiter = this->size()*10; } - //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) - std::cout << "Algo 5.2 with newton_tol = " << newton_tol << ", cg_tol = " - << cg_tol << " maxiter = " << maxiter << " and ΔF =" < void SystemBase::cg(Vector& delF0, const Vector& rhs, - const Real cg_tol, const Uint maxiter) { + const Real cg_tol, const Uint maxiter, + bool verbose) { //implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p 112) // residual vector: set r_0 <-Ax_0 - b Vector r_vec = this->apply_K(VecMap(delF0.data(), delF0.size())); this->fft_engine->convolve(this->Ghats.data(), r_vec.data(), r_vec.data()); r_vec -= rhs; // search direction: set p_0 = -r_0 Vector p_vec = -r_vec; const Real cg_norm = cg_tol*cg_tol; Real rTr = r_vec.transpose()*r_vec; //precompute this Vector Apk(p_vec.size()); for (size_t cg_step = 0; (cg_step < maxiter) && (rTr > cg_norm); ++cg_step) { //compute step length; Apk = this->apply_K(VecMap(p_vec.data(), p_vec.size())); this->fft_engine->convolve(this->Ghats.data(), Apk.data(), Apk.data()); Real alpha = rTr/(p_vec.transpose()*Apk); delF0 += alpha * p_vec; r_vec += alpha * Apk; Real tmp_rTr = r_vec.transpose()*r_vec; Real beta = tmp_rTr/rTr; rTr = tmp_rTr; - std::cout << "at CG step " << cg_step << ": |r| = " << std::sqrt(rTr) - << ", cg_tol = " << cg_tol << std::endl; + + if (verbose) { + std::cout << "at CG step " << cg_step << ": |r| = " << std::sqrt(rTr) + << ", cg_tol = " << cg_tol << std::endl; + } p_vec = -r_vec + beta*p_vec; } if (rTr > cg_norm) { throw std::runtime_error("cg did not converge"); } } //----------------------------------------------------------------------------// template typename SystemBase::Vector SystemBase::apply_K(const VecMap& in) { const T2VecMap P(const_cast(in.data()), DimM, DimM, static_cast(this->size())); const std::array,2> - prod_dims {Eigen::IndexPair{2, 1}, - Eigen::IndexPair{3, 0}}; + prod_dims {Eigen::IndexPair{2, 0}, + Eigen::IndexPair{3, 1}}; Vector retval(P.size()); T2VecMap retmap(retval.data(), DimM, DimM, this->size()); for (size_t pix_id = 0; pix_id < this->size(); ++pix_id) { retmap.chip<2>(pix_id) = this->Ks.template chip<4>(pix_id).contract(P.chip<2>(pix_id), prod_dims); } return retval; } //----------------------------------------------------------------------------// template void SystemBase::build_ghats() { //! row-major storage for pixels Projection::FreqStruc freqs(this->sizes, this->nb_pixels); auto it = freqs.begin(); //drop the first frequency, as it is zero ++it; for (size_t pix_id = 1; pix_id < this->nb_pixel; ++pix_id, ++it) { auto xi = *it; xi/=xi.norm(); using T4shape = Eigen::Sizes; using T4 = Eigen::TensorFixedSize; Eigen::TensorMap Ghat_kk(&this->Ghats(0, 0, 0, 0, pix_id), DimS, DimS, DimS, DimS); for (Dim_t im = 0; im < DimS; ++im) { for (Dim_t j = 0; j < DimS; ++j) { for (Dim_t l = 0; l < DimS; ++l) { Ghat_kk(im, j, l, im) = xi(j)*xi(l); } } } } } //----------------------------------------------------------------------------// template void SystemBase::init_constitutive_laws() { for (const auto & mat_ptr: this->materials) { mat_ptr->initialize(); } } //----------------------------------------------------------------------------// template void SystemBase::eval_constitutive_laws() { for (const auto & mat_ptr: this->materials) { //using SecondArrayMap = typename MaterialBase::SecondArrayMap; //using FourthArrayMap = typename MaterialBase::FourthArrayMap; mat_ptr->compute_First_Piola_Kirchhoff (this->grad, this->stress, this->Ks, this->nb_pixels); } } //----------------------------------------------------------------------------// template void SystemBase::add_material(mat_Ptr material) { this->materials.push_back(material); } //----------------------------------------------------------------------------// template typename SystemBase::Index SystemBase::iterator::operator*() { size_t tmp_index = this->index; Index retval; for (Dim_t i = DimS-1; i > 0; --i) { retval[i] = tmp_index%this->parent.nb_pixels[i]; tmp_index/=this->parent.nb_pixels[i]; } retval[0] = tmp_index; return retval; } template class SystemBase<2, 2>; template class SystemBase<2, 3>; template class SystemBase<3, 3>; } // muSpectre diff --git a/src/system/system_base.hh b/src/system/system_base.hh index 0ed3ec9..d492986 100644 --- a/src/system/system_base.hh +++ b/src/system/system_base.hh @@ -1,155 +1,156 @@ /** * file system_base.hh * * @author Till Junge * * @date 09 May 2017 * * @brief Base class representing a periodic cell and it's solution * * @section LICENCE * * Copyright (C) 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 "common/common.hh" #include "materials/material.hh" #include "system/fft_engine.hh" #ifndef SYSTEM_BASE_H #define SYSTEM_BASE_H namespace muSpectre { namespace Projection { using vector = Eigen::Matrix; //----------------------------------------------------------------------------// //! computes fft frequecies. see function fftfreq in numpy's helper.py vector fft_freqs(Uint nb_samples, Real length); vector fft_freqs(Uint nb_samples); template class FreqStruc { public: FreqStruc(const std::array& sizes, const Index_t& nb_pixels); class iterator { public: iterator(const FreqStruc & parent, Index_t index); inline Eigen::Matrix operator*(); inline iterator& operator++ (); inline bool operator!=(const iterator & other) const {return this->index !=other.index;} protected: inline Index_t get_index(); const FreqStruc & parent; size_t index; }; inline iterator begin(); inline iterator end(); protected: const Index_t& nb_pixels; const std::array& sizes; //! scaled frequencies (k_i/L_i) std::array xi; }; } // Projection //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class SystemBase { public: using Index = Index_t; using T2 = Eigen::TensorFixedSize>; using T4Vector = Eigen::Tensor; using T2Vector = Eigen::Tensor; using T4VecMap = Eigen::TensorMap; using T2VecMap = Eigen::TensorMap; using Vector = Eigen::Matrix; using VecMap = Eigen::Map; using FFT_Ptr = std::shared_ptr>; using mat_Ptr = std::shared_ptr>; enum class Formulation{finite_strain, small_strain}; class iterator { public: iterator(const SystemBase & parent, bool begin=true) :parent(parent), index (begin ? 0 : parent.size()){} Index operator*(); iterator& operator++() {index++; return *this;}; bool operator!=(const iterator & other) const{return index!= other.index;}; protected: const SystemBase & parent; size_t index; }; SystemBase(std::array sizes, Index nb_pixels, Formulation = Formulation::finite_strain); virtual ~SystemBase() = default; void set_fft_engine(FFT_Ptr fft_eng); void add_material(mat_Ptr material); - void solve(const T2 & DeltaF, const Real newton_tol, const Real cg_tol, Uint maxiter=0); + void solve(const T2 & DeltaF, const Real newton_tol, const Real cg_tol, Uint maxiter=0, + bool newton_verbose = false, bool cg_verbose = false); const Index & get_nb_pixels() const {return this->nb_pixels;} inline size_t size() const {return this->nb_pixel;} iterator begin(){return iterator (*this);} iterator end(){return iterator (*this, false);} protected: void build_ghats(); virtual void init_constitutive_laws(); virtual void eval_constitutive_laws(); - void cg(Vector & delF0, const Vector& rhs, const Real cg_tol, const Uint maxiter); + void cg(Vector & delF0, const Vector& rhs, const Real cg_tol, const Uint maxiter, bool verbose); Vector apply_K (const VecMap & in); T4Vector get_scaled_freqs; std::array sizes; Index nb_pixels; size_t nb_pixel; Formulation form; FFT_Ptr fft_engine; std::vector materials; //! projection operator T4Vector Ghats; //! tangent stiffness tensors by pixel; T4Vector Ks; //! gradient or strain (depending on finite/small strain formulation) T2Vector grad; //! Cauchy or First Piola-Kirchhoff stress T2Vector stress; }; } // muSpectre #endif /* SYSTEM_BASE_H */ diff --git a/tests/test_fftw_comparison.cc b/tests/test_fftw_comparison.cc new file mode 100644 index 0000000..963d75d --- /dev/null +++ b/tests/test_fftw_comparison.cc @@ -0,0 +1,91 @@ +/** + * file test_fftw_c2c.cc + * + * @author Till Junge + * + * @date 11 May 2017 + * + * @brief test fftw redundant class + * + * @section LICENCE + * + * Copyright (C) 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 "system/fftw_engine_c2c.hh" +#include "system/fftw_engine_r2c.hh" +#include "tests.hh" +#include "materials/materials_toolbox.hh" + +namespace muSpectre { + template + struct fftw_comp_fixture{ + fftw_comp_fixture() + :fftwc(std_nb_pix, FFT_PlanFlags::estimate), + fftwr(std_nb_pix, FFT_PlanFlags::estimate){} + + const static typename FFTW_EngineC2C::Index std_nb_pix; + FFTW_EngineC2C fftwc; + FFTW_EngineR2C fftwr; + }; + + template<> const Index_t<1> fftw_comp_fixture<1,1>::std_nb_pix{3}; + template<> const Index_t<2> fftw_comp_fixture<2,2>::std_nb_pix{5,3}; + template<> const Index_t<2> fftw_comp_fixture<2,1>::std_nb_pix{5,3}; + template<> const Index_t<3> fftw_comp_fixture<3,3>::std_nb_pix{7,5,3}; + + using fix2 = fftw_comp_fixture<2, 2>; + BOOST_FIXTURE_TEST_CASE(fftw_comparison, fix2) { + + const size_t dim_m = 2; + const size_t tot_nb_pixels = fftwc.get_tot_nb_pixels(); + Real in_arr[tot_nb_pixels*dim_m*dim_m]; + for (size_t i = 0; i < tot_nb_pixels; ++i) { + for (size_t d = 0; d < dim_m*dim_m; ++d) { + in_arr[i*dim_m*dim_m+d] = i; + } + } + + Eigen::TensorG(dim_m, dim_m, dim_m, dim_m, tot_nb_pixels); + G.setZero(); + for (size_t i = 1; i < tot_nb_pixels; ++i) { + Eigen::TensorMap> Gmap(G.data()+dim_m*dim_m*dim_m*dim_m*i, dim_m, dim_m,dim_m,dim_m); + Gmap = Tensors::I4(); + } + + Eigen::Tensor resc(dim_m, dim_m, tot_nb_pixels); + Eigen::Tensor resr(dim_m, dim_m, tot_nb_pixels); + fftwc.convolve(G.data(), in_arr, resc.data()); + fftwr.convolve(G.data(), in_arr, resr.data()); + + Eigen::Mapcmap(resc.data(), resc.size()); + Eigen::Maprmap(resr.data(), resr.size()); + Real error = (rmap-cmap).norm()/rmap.norm(); + Real tol = 1e-14; + + BOOST_CHECK_LT(error, tol); + + } + +} // muSpectre diff --git a/tests/test_material_hyperelastic.cc b/tests/test_material_hyperelastic.cc index 0453096..5e1cdac 100644 --- a/tests/test_material_hyperelastic.cc +++ b/tests/test_material_hyperelastic.cc @@ -1,216 +1,216 @@ /** * file test_material_hyperelastic.cc * * @author Till Junge * * @date 01 May 2017 * * @brief Test the hyper-elastic material used in Section 4 of * de Geus et al. / Comput. Methods in Appl. Mech. Engrg. * 318, (2017), 412–430 https://doi.org/10.1016/j.cma.2016.12.032 * * @section LICENCE * * Copyright (C) 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_hyper_elastic.hh" #include #include #include #include "tests.hh" namespace muSpectre { // params from "material_hyper_elastic.ipynb"; const Real E = 1.9058231542461002; const Real nu = 0.2993067590987868; template struct mat_fixture: public MaterialHyperElastic{ mat_fixture(): MaterialHyperElastic(E, nu){ }; const static Dim_t DimS = s_dim; const static Dim_t DimM = m_dim; /** inputs with expected outputs computed with * "material_hyper_elastic.ipynb" in folder "helpers"*/ using parent = MaterialHyperElastic; using SecondTens = typename parent::SecondTens; using FourthVoigt = Eigen::Matrix; const static SecondTens F, P; const static FourthVoigt K; }; //----------------------------------------------------------------------------// template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::F = (mat_fixture<2, 2>::SecondTens() << 1.0, 0.1, 0.2, 1.0).finished(); template <> const mat_fixture<2, 2>::SecondTens mat_fixture<2, 2>::P = (mat_fixture<2, 2>::SecondTens()<< 0.07868216666666673, 0.22348781666666673, 0.2313560333333334, 0.07868216666666672).finished(); template <> const mat_fixture<2, 2>::FourthVoigt mat_fixture<2, 2>::K = (mat_fixture<2, 2>::FourthVoigt()<< 2.624580833333334, 1.1084346666666667, 0.40273666666666674, 0.5854533333333334, 1.1084346666666667, 2.6245808333333334, 0.40273666666666674, 0.5854533333333334, 0.5854533333333334, 0.5854533333333334, 0.7552753333333334, 0.8925028333333335, 0.40273666666666674, 0.40273666666666674, 0.7936838333333334, 0.7552753333333334).finished(); //----------------------------------------------------------------------------// template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::F = (mat_fixture<2, 3>::SecondTens()<< 1.0, 0.1, 0.0, 0.2, 1.0, 0.0, 0.0, 0.0, 1.0).finished(); template <> const mat_fixture<2, 3>::SecondTens mat_fixture<2, 3>::P = (mat_fixture<2, 3>::SecondTens()<< 0.07868216666666673, 0.22348781666666673, 0.0, 0.2313560333333334, 0.07868216666666672, 0.0, 0.0, 0.0, 0.027344166666666694).finished(); template <> const mat_fixture<2, 3>::FourthVoigt mat_fixture<2, 3>::K = (mat_fixture<2, 3>::FourthVoigt()<<2.624580833333334, 1.1084346666666667, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334, 1.1084346666666667, 2.6245808333333334, 1.0937666666666668, 0.0, 0.0, 0.40273666666666674, 0.0, 0.0, 0.5854533333333334, 1.0937666666666668, 1.0937666666666668, 2.5879108333333334, 0.0, 0.0, 0.10937666666666668, 0.0, 0.0, 0.21875333333333336, 0.0, 0.0, 0.0, 0.7334, 0.07334, 0.0, 0.7680781666666667, 0.22002000000000005, 0.0, 0.0, 0.0, 0.0, 0.14668, 0.7334, 0.0, 0.22002000000000005, 0.7900801666666668, 0.0, 0.5854533333333334, 0.5854533333333334, 0.21875333333333336, 0.0, 0.0, 0.7552753333333334, 0.0, 0.0, 0.8925028333333335, 0.0, 0.0, 0.0, 0.7900801666666668, 0.22002, 0.0, 0.7334, 0.14668, 0.0, 0.0, 0.0, 0.0, 0.22002, 0.7680781666666667, 0.0, 0.07334, 0.7334, 0.0, 0.40273666666666674, 0.40273666666666674, 0.10937666666666668, 0.0, 0.0, 0.7936838333333334, 0.0, 0.0, 0.7552753333333334).finished(); //----------------------------------------------------------------------------// template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::F = (mat_fixture<3, 3>::SecondTens()<< 1.0, 0.1, 0.3, 0.2, 1.0, 0.4, 0.5, 0.6, 1.0 ).finished(); template <> const mat_fixture<3, 3>::SecondTens mat_fixture<3, 3>::P = (mat_fixture<3, 3>::SecondTens()<< 0.9479714333333337, 0.7435627833333334, 0.92523635, 0.8402667666666668, 1.1591906333333337, 1.1568859333333334, 1.264590916666667, 1.4368351000000001, 1.4569510333333335).finished(); template <> const mat_fixture<3, 3>::FourthVoigt mat_fixture<3, 3>::K = (mat_fixture<3,3>::FourthVoigt() << 3.3442565, 1.1084346666666667, 1.2037766666666667, 0.4815106666666667, 1.193542, 0.6227566666666668, 0.69293, 1.5443073333333335, 0.6734613333333335, 1.1084346666666667, 3.4762685000000006, 1.2697826666666667, 1.4862686666666667, 0.35746600000000006, 0.4907446666666667, 1.90304, 0.6348913333333334, 0.8054733333333335, 1.2037766666666667, 1.2697826666666667, 3.6889545000000004, 1.5376066666666668, 1.178874, 0.2413886666666667, 1.851702, 1.5589753333333336, 0.3654333333333334, 0.69293, 1.90304, 1.851702, 0.9959040000000001, 0.270218, 0.7403540000000001, 2.6075758333333336, 0.9881900000000001, 0.49795200000000006, 1.5443073333333335, 0.6348913333333334, 1.5589753333333336, 0.3654333333333334, 0.8974650000000001, 0.4947283333333334, 0.9881900000000001, 2.3479155, 0.9894566666666668, 0.6734613333333335, 0.8054733333333335, 0.3654333333333334, 0.7915653333333335, 0.358986, 0.7552753333333334, 0.49795200000000006, 0.9894566666666668, 1.6635165000000003, 0.4815106666666667, 1.4862686666666667, 1.5376066666666668, 1.8534405000000003, 0.5272880000000001, 0.2637706666666667, 0.9959040000000001, 0.3654333333333334, 0.7915653333333335, 1.193542, 0.35746600000000006, 1.178874, 0.5272880000000001, 1.6521988333333337, 0.810217, 0.270218, 0.8974650000000001, 0.358986, 0.6227566666666668, 0.4907446666666667, 0.2413886666666667, 0.2637706666666667, 0.810217, 1.5940335000000003, 0.7403540000000001, 0.4947283333333334, 0.7552753333333334).finished(); typedef boost::mpl::list, mat_fixture<2, 3>, mat_fixture<3, 3>> test_types; BOOST_FIXTURE_TEST_CASE_TEMPLATE(material_hyper_elastic_test, F, test_types, F) { BOOST_CHECK_CLOSE(F::lambda, E*nu/((1+nu)*(1-2*nu)), tol); BOOST_CHECK_CLOSE(F::mu, E/2/(1+nu), tol); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_elastic_tensor, F, test_types, F) { Eigen::Matrix Cv; F::VoigtConv::fourth_to_voigt(F::C, Cv); auto & E = F::Young; auto & nu = F::Poisson; auto prefactor = E/(1+nu)/(1-2*nu); const auto diff{vsize(F::DimM) - F::DimM}; for (Dim_t i = 0; i < F::DimM; ++i) { BOOST_CHECK_CLOSE(Cv(i, i ), prefactor*(1-nu), 1e-12); for (Dim_t j = i + 1; j < F::DimM; ++j ) { BOOST_CHECK_CLOSE(Cv(i, j), prefactor*nu, 1e-12); BOOST_CHECK_CLOSE(Cv(j, i), prefactor*nu, 1e-12); } } for (Dim_t i = 0; i < diff; ++i) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, i+F::DimM), prefactor*(1-2*nu)/2, 1e-12); for (Dim_t j = i + 1; j < diff; ++j ) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, j+F::DimM), 0, 1e-12); } for (Dim_t j = 0; j < diff; ++j) { BOOST_CHECK_CLOSE(Cv(i+F::DimM, j ), 0, 1e-12); BOOST_CHECK_CLOSE(Cv(i, j+F::DimM), 0, 1e-12); } } } //----------------------------------------------------------------------------// BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constitutive_law, F, test_types, F) { //using SecondArray = typename MaterialHyperElastic::SecondArray; //using FourthArray = typename MaterialHyperElastic::FourthArray; using SecondTens = typename MaterialHyperElastic::SecondTens; using FourthTens = typename MaterialHyperElastic::FourthTens; //auto F_t = getMap(add_dim(F::F)); //auto K_t = Eigen::TensorMap(const_cast(F::K.data()), F::DimM, F::DimM, F::DimM, F::DimM, 1, 1); //std::cout << F_t; //std::cout << K_t; this->add_pixel(Index_t{0}); this->initialize(); std::array sizes; for (auto & s: sizes) { s=1; } auto F_t = F::get_second_array_map(F::F.data(), sizes); //std::cout << "hello: " << std::endl << F_t << std::endl; SecondTens P_; FourthTens K_; auto P_t = F::get_second_array_map(P_.data(), sizes); auto K_t = F::get_fourth_array_map(K_.data(), sizes); for (Dim_t i = 0; i < F::DimM; ++i) { for (Dim_t j = 0; j < F::DimM; ++j) { P_(i, j) = 10*(i+1)+j+1; } } F::compute_First_Piola_Kirchhoff_stress(F_t, P_t, K_t); - auto K_v = VoigtConversion::template fourth_to_voigt(K_); + auto K_v = VoigtConversion::template fourth_to_voigt(K_.shuffle(std::array{1,0,2,3})); Real P_error = (P_-F::P).norm(); if (P_error >= tol) { std::cout << "First Piola-Kirhoff stress wrong:" << std::endl << "P ref = " << std::endl << F::P << std::endl << "P comp = " << std::endl << P_ << std::endl << "P diff = " << std::endl << F::P-P_ << std::endl<< std::endl; } BOOST_CHECK_LT(P_error, tol); Real K_error = (K_v - F::K).norm(); if (K_error >= tol) { std::cout << "Tangent wrong: "<< std::endl << "K ref = " << std::endl << F::K < struct sys_fixture:public SystemBase { sys_fixture() :SystemBase(std_size, std_nb_pix){} const static std::array std_size; const static Index_t std_nb_pix; const static Dim_t dim_s=DimS, dim_m=DimM; }; template<> const std::array sys_fixture<2, 2>::std_size{5, 3}; template<> const std::array sys_fixture<2, 3>::std_size{5, 3}; template<> const std::array sys_fixture<3, 3>::std_size{7, 5, 3}; template<> const Index_t<2> sys_fixture<2, 2>::std_nb_pix{5, 3}; template<> const Index_t<2> sys_fixture<2, 3>::std_nb_pix{5, 3}; template<> const Index_t<3> sys_fixture<3, 3>::std_nb_pix{7, 5, 3}; #include "testdata_system_base.hh" //----------------------------------------------------------------------------// using test_systems = boost::mpl::list, sys_fixture<2, 3>, sys_fixture<3, 3>>; using dim_test_systems = boost::mpl::list, sys_fixture<3, 3>>; //----------------------------------------------------------------------------// BOOST_AUTO_TEST_CASE(fftfreq_test) { auto length = 15.2; auto n1{8}, n2{9}; Eigen::RowVectorXd ref1(n1); ref1 << 0.0, 0.06578947368421052, 0.13157894736842105, 0.19736842105263158, -0.2631578947368421, -0.19736842105263158, -0.13157894736842105, -0.06578947368421052; Eigen::RowVectorXd ref2(n2); ref2 << 0.0, 0.06578947368421052, 0.13157894736842105, 0.19736842105263158, 0.2631578947368421, -0.2631578947368421, -0.19736842105263158, -0.13157894736842105, -0.06578947368421052; auto challenge1(Projection::fft_freqs(n1, length)); auto challenge2(Projection::fft_freqs(n2, length)); Eigen::RowVectorXd er1 = ref1-challenge1; Real error1 = (ref1 - challenge1).norm(); Real error2 = (ref2 - challenge2).norm(); BOOST_CHECK_LT(error1, tol); BOOST_CHECK_LT(error2, tol); } //----------------------------------------------------------------------------// BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, F, test_systems, F) { } //----------------------------------------------------------------------------// BOOST_FIXTURE_TEST_CASE_TEMPLATE(Ghat_test, F, test_systems, F) { const int dim_s = F::dim_s; const int dim_m = F::dim_m; auto ghat_ref = Ghats(); //TODO; for unknown reason, the python ref gives different results when material and spatial dimension dffer. not sure who's wrong if(dim_s == dim_m) { for (size_t pix_id = 0; pix_id < F::nb_pixel; ++pix_id) { using Tshape = Eigen::Sizes; using Ttype = Eigen::TensorFixedSize; Ttype ghat_i = F::Ghats.template chip<4>(pix_id); Real error = (VoigtConversion::fourth_to_2d(ghat_i) -(*ghat_ref[pix_id])).norm(); if (error >= tol) { std::cout <<"For dim_m = " << dim_m << ", dim_s = " << dim_s << " reference is:\n" << VoigtConversion::fourth_to_2d(ghat_i) << "\nbut I got:\n" << *ghat_ref[pix_id] << std::endl; } BOOST_CHECK_LT(error, tol); } } } //----------------------------------------------------------------------------// BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_uniform_patchtest_hyperelastic, F, dim_test_systems, F) { const int dim_s = F::dim_s; const int dim_m = F::dim_m; typename F::FFT_Ptr fft_eng = std::make_shared> (F::std_nb_pix, FFT_PlanFlags::estimate); const Real E = 1.0030648180242636; const Real nu = 0.29930675909878679; typename F::mat_Ptr hyperelastic = std::make_shared> (E, nu); F::set_fft_engine(fft_eng); F::add_material(hyperelastic); Real newton_tol = 1e-4; Real cg_tol = 1e-4; size_t maxiter = 12; auto it = F::begin(); auto end = F::end(); for (; it != end; ++it) { hyperelastic->add_pixel(*it); } typename F::T2 DeltaF; DeltaF.setZero(); DeltaF(0, 1) = 1e-2; - F::solve(DeltaF, newton_tol, cg_tol, maxiter); + bool verbose = false; + F::solve(DeltaF, newton_tol, cg_tol, maxiter, verbose); } + //----------------------------------------------------------------------------// + BOOST_FIXTURE_TEST_CASE_TEMPLATE(dumb_uniform_patchtest_hyperelastic_c2c, F, dim_test_systems, F) { + const int dim_s = F::dim_s; + const int dim_m = F::dim_m; + typename F::FFT_Ptr fft_eng = + std::make_shared> + (F::std_nb_pix, FFT_PlanFlags::estimate); + + const Real E = 1.0030648180242636; + const Real nu = 0.29930675909878679; + + typename F::mat_Ptr hyperelastic = + std::make_shared> + (E, nu); + F::set_fft_engine(fft_eng); + F::add_material(hyperelastic); + Real newton_tol = 1e-4; + Real cg_tol = 1e-4; + size_t maxiter = 12; + + + auto it = F::begin(); + auto end = F::end(); + for (; it != end; ++it) { + hyperelastic->add_pixel(*it); + } + + typename F::T2 DeltaF; + DeltaF.setZero(); + DeltaF(0, 1) = 1e-2; + + bool verbose = false; + F::solve(DeltaF, newton_tol, cg_tol, maxiter, verbose); + + } + + + + //----------------------------------------------------------------------------// + BOOST_AUTO_TEST_CASE(non_uniform_test_hyperelastic_complex) { + const Int dim = 3; + const Int N = 5; + const std::array sizes{N, N, N}; + const std::array nb_pixels{N, N, N}; + using Sys_t = SystemBase; + SystemBase system(sizes, nb_pixels); + + typename Sys_t::FFT_Ptr fft_eng = + std::make_shared> + (nb_pixels, FFT_PlanFlags::estimate); + + + const Real E = 1.0030648180242636; + const Real nu = 0.29930675909878679; + + typename Sys_t::mat_Ptr hyperelastic_soft = + std::make_shared> + (E, nu); + + typename Sys_t::mat_Ptr hyperelastic_hard = + std::make_shared> + (10*E, nu); + + + system.set_fft_engine(fft_eng); + system.add_material(hyperelastic_soft); + system.add_material(hyperelastic_hard); + Real newton_tol = 1e-5; + Real cg_tol = 1e-8; + size_t maxiter = 35; + + + bool first = true; + for (const auto && pixel:system) { + if (first) { + first = false; + hyperelastic_hard->add_pixel(pixel); + } else { + hyperelastic_soft->add_pixel(pixel); + } + } + + typename Sys_t::T2 DeltaF; + DeltaF.setZero(); + DeltaF(0, 1) = 1.; + + bool verbose = true; + system.solve(DeltaF, newton_tol, cg_tol, maxiter, verbose); + //BOOST_CHECK_THROW(system.solve(DeltaF, newton_tol, cg_tol, maxiter, verbose), + // std::runtime_error); + } //----------------------------------------------------------------------------// BOOST_AUTO_TEST_CASE(non_uniform_test_hyperelastic) { const Int dim = 3; const Int N = 5; const std::array sizes{N, N, N}; const std::array nb_pixels{N, N, N}; using Sys_t = SystemBase; SystemBase system(sizes, nb_pixels); typename Sys_t::FFT_Ptr fft_eng = std::make_shared> (nb_pixels, FFT_PlanFlags::estimate); const Real E = 1.0030648180242636; const Real nu = 0.29930675909878679; typename Sys_t::mat_Ptr hyperelastic_soft = std::make_shared> (E, nu); typename Sys_t::mat_Ptr hyperelastic_hard = std::make_shared> (10*E, nu); system.set_fft_engine(fft_eng); system.add_material(hyperelastic_soft); system.add_material(hyperelastic_hard); Real newton_tol = 1e-5; Real cg_tol = 1e-8; size_t maxiter = 35; bool first = true; for (const auto && pixel:system) { if (first) { first = false; hyperelastic_hard->add_pixel(pixel); } else { hyperelastic_soft->add_pixel(pixel); } } typename Sys_t::T2 DeltaF; DeltaF.setZero(); DeltaF(0, 1) = 1.; - BOOST_CHECK_THROW(system.solve(DeltaF, newton_tol, cg_tol, maxiter), - std::runtime_error); + bool verbose = true; + system.solve(DeltaF, newton_tol, cg_tol, maxiter, verbose); + //BOOST_CHECK_THROW(system.solve(DeltaF, newton_tol, cg_tol, maxiter, verbose), + // std::runtime_error); } } // muSpectre