diff --git a/CMakeLists.txt b/CMakeLists.txt index 2670f99..d40ee97 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,140 +1,141 @@ # ============================================================================= # 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.7.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) add_compile_options(-Wall -Wextra -Weffc++) enable_testing() set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) include(muspectreTools) if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") # 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() find_package(Boost COMPONENTS unit_test_framework REQUIRED) add_external_package(Eigen3 VERSION 3.3 CONFIG) add_external_package(pybind11 VERSION 2.2 CONFIG) #find_package(Eigen3 3.3 REQUIRED NO_MODULE) include_directories( ${CMAKE_SOURCE_DIR}/tests ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR} ) #find_package(MPI REQUIRED) find_package(FFTW REQUIRED) #find_package(FFTWMPI REQUIRED) #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) add_test(main_test_suite main_test_suite --report_level=detailed --build_info=TRUE) #add_test(main_test_suite main_test_suite --build_info=TRUE) # compile the library add_compile_options( -Werror) add_subdirectory( ${CMAKE_SOURCE_DIR}/src/ ) add_subdirectory( ${CMAKE_SOURCE_DIR}/language_bindings/ ) #compile executables file( GLOB binaries "${CMAKE_SOURCE_DIR}/bin/*.cc") 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}) # compile 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}) # copy python test #build tests 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}) add_test(python_binding_test python_binding_tests.py) diff --git a/bin/small_case.py b/bin/small_case.py index 745bb3c..eaae229 100755 --- a/bin/small_case.py +++ b/bin/small_case.py @@ -1,92 +1,92 @@ #!/usr/bin/env python3 """ file small_case.py @author Till Junge @date 12 Jan 2018 @brief small case for debugging @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. """ import sys import os import numpy as np sys.path.append(os.path.join(os.getcwd(), "language_bindings/python")) import pyMuSpectre as µ -resolution = [555, 555] +resolution = [251, 251] center = np.array([r//2 for r in resolution]) incl = resolution[0]//5 lengths = [7., 5.] formulation = µ.Formulation.small_strain rve = µ.SystemFactory(resolution, lengths, formulation) hard = µ.material.MaterialHooke2d.make( rve, "hard", 210e9, .33) soft = µ.material.MaterialHooke2d.make( rve, "soft", 70e9, .33) for i, pixel in enumerate(rve): - #if np.linalg.norm(center - np.array(pixel)) * * @date 09 Jan 2018 * * @brief python bindings for the muSpectre solvers * * @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. */ #include "common/common.hh" #include "solver/solvers.hh" #include #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * Solvers instanciated for systems with equal spatial and material dimension */ + + template void add_newton_cg_helper(py::module & mod) { const char name []{"newton_cg"}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; using grad_vec = GradIncrements; mod.def(name, [](sys & s, const grad & g, Real ct, Real nt, Uint max, Dim_t verb) -> OptimizeResult { return newton_cg(s, g, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); mod.def(name, [](sys & s, const grad_vec & g, Real ct, Real nt, Uint max, Dim_t verb) -> std::vector { return newton_cg(s, g, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); } template void add_de_geus_helper(py::module & mod) { const char name []{"de_geus"}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; using grad_vec = GradIncrements; mod.def(name, [](sys & s, const grad & g, Real ct, Real nt, Uint max, Dim_t verb) -> OptimizeResult { return de_geus(s, g, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); mod.def(name, [](sys & s, const grad_vec & g, Real ct, Real nt, Uint max, Dim_t verb) -> std::vector { return de_geus(s, g, ct, nt, max, verb); }, "system"_a, "ΔF₀"_a, "cg_tol"_a, "newton_tol"_a, "maxiter"_a=0, "verbose"_a=0); } template void add_solver_helper(py::module & mod) { add_newton_cg_helper(mod); add_de_geus_helper (mod); } void add_solvers(py::module & mod) { auto solvers{mod.def_submodule("solvers")}; solvers.doc() = "bindings for solvers"; py::class_(mod, "OptimizeResult") .def_readwrite("grad", &OptimizeResult::grad) .def_readwrite("stress", &OptimizeResult::stress) .def_readwrite("success", &OptimizeResult::success) .def_readwrite("status", &OptimizeResult::status) .def_readwrite("message", &OptimizeResult::message) .def_readwrite("nb_it", &OptimizeResult::nb_it) .def_readwrite("nb_fev", &OptimizeResult::nb_fev); add_solver_helper(solvers); add_solver_helper(solvers); } diff --git a/src/common/field.hh b/src/common/field.hh index 8781b84..205c53e 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,797 +1,807 @@ /** * file field.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief header-only implementation of a field for field collections * * @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. */ #ifndef FIELD_H #define FIELD_H #include "common/T4_map_proxy.hh" #include #include #include #include #include #include #include #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ class FieldCollectionError: public std::runtime_error { public: explicit FieldCollectionError(const std::string& what) :std::runtime_error(what){} explicit FieldCollectionError(const char * what) :std::runtime_error(what){} }; class FieldError: public FieldCollectionError { using Parent = FieldCollectionError; public: explicit FieldError(const std::string& what) :Parent(what){} explicit FieldError(const char * what) :Parent(what){} }; class FieldInterpretationError: public FieldError { public: explicit FieldInterpretationError(const std::string & what) :FieldError(what){} explicit FieldInterpretationError(const char * what) :FieldError(what){} }; namespace internal { /* ---------------------------------------------------------------------- */ template class FieldBase { protected: //! constructor //! unique name (whithin Collection) //! number of components //! collection to which this field belongs (eg, material, system) FieldBase(std::string unique_name, size_t nb_components, FieldCollection & collection); public: using collection_t = FieldCollection; //! Copy constructor FieldBase(const FieldBase &other) = delete; //! Move constructor FieldBase(FieldBase &&other) = delete; //! Destructor virtual ~FieldBase() = default; //! Copy assignment operator FieldBase& operator=(const FieldBase &other) = delete; //! Move assignment operator FieldBase& operator=(FieldBase &&other) = delete; /* ---------------------------------------------------------------------- */ //!Identifying accessors //! return field name inline const std::string & get_name() const; //! return field type //inline const Field_t & get_type() const; //! return my collection (for iterating) inline const FieldCollection & get_collection() const; //! return my collection (for iterating) inline const size_t & get_nb_components() const; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const = 0; virtual size_t size() const = 0; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() = 0; //! give access to collections friend FieldCollection; friend typename FieldCollection::Parent; protected: /* ---------------------------------------------------------------------- */ //! allocate memory etc virtual void resize(size_t size) = 0; const std::string name; const size_t nb_components; const FieldCollection & collection; private: }; /** * dummy intermediate class to provide a run-time polymorphic * typed field. Mainly for binding python */ template class TypedFieldBase: public FieldBase { public: using Parent = FieldBase; using collection_t = typename Parent::collection_t; using Scalar = T; using Base = Parent; using EigenRep = Eigen::Array; using EigenMap = Eigen::Map; + using EigenVec = Eigen::Map; //! Default constructor TypedFieldBase() = delete; TypedFieldBase(std::string unique_name, size_t nb_components, FieldCollection& collection); //! Copy constructor TypedFieldBase(const TypedFieldBase &other) = delete; //! Move constructor TypedFieldBase(TypedFieldBase &&other) = delete; //! Destructor virtual ~TypedFieldBase() = default; //! Copy assignment operator TypedFieldBase& operator=(const TypedFieldBase &other) = delete; //! Move assignment operator TypedFieldBase& operator=(TypedFieldBase &&other) = delete; //! return type_id of stored type virtual const std::type_info & get_stored_typeid() const override final; virtual size_t size() const override = 0; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) virtual void set_zero() override = 0; virtual T* data() = 0; virtual const T* data() const = 0; EigenMap eigen(); + EigenVec eigenvec(); protected: private: }; /* ---------------------------------------------------------------------- */ //! declaraton for friending template class FieldMap; /* ---------------------------------------------------------------------- */ template class TypedSizedFieldBase: public TypedFieldBase { friend class FieldMap; friend class FieldMap; public: constexpr static auto nb_components{NbComponents}; using Parent = TypedFieldBase; using collection_t = typename Parent::collection_t; using Scalar = T; using Base = typename Parent::Base; //using storage_type = Eigen::Array; using StoredType = Eigen::Array; using StorageType = std::conditional_t >, std::vector>>; using EigenRep = Eigen::Array; using EigenMap = std::conditional_t< ArrayStore, Eigen::Map>, Eigen::Map>; using ConstEigenMap = std::conditional_t< ArrayStore, Eigen::Map>, Eigen::Map>; TypedSizedFieldBase(std::string unique_name, FieldCollection& collection); virtual ~TypedSizedFieldBase() = default; //! initialise field to zero (do more complicated initialisations through //! fully typed maps) inline void set_zero() override final; template inline void push_back(const std::enable_if_t & value); template inline std::enable_if_t push_back(const StoredType & value); //! Number of stored arrays (i.e. total number of stored //! scalars/NbComponents) size_t size() const override final; static TypedSizedFieldBase & check_ref(Base & other); static const TypedSizedFieldBase & check_ref(const Base & parent); inline T* data() override final {return this->get_ptr_to_entry(0);} inline const T* data() const override final {return this->get_ptr_to_entry(0);} inline EigenMap eigen(); inline ConstEigenMap eigen() const; inline typename Parent::EigenMap dyn_eigen() {return Parent::eigen();} template inline Real inner_product(const TypedSizedFieldBase & other) const; protected: template inline std::enable_if_t get_ptr_to_entry(const size_t&& index); template inline std::enable_if_t get_ptr_to_entry(const size_t&& index) const; template inline T* get_ptr_to_entry(std::enable_if_t index); template inline const T* get_ptr_to_entry(std::enable_if_t index) const; inline virtual void resize(size_t size) override final; StorageType values{}; }; } // internal /* ---------------------------------------------------------------------- */ template class TensorField: public internal::TypedSizedFieldBase { public: using Parent = internal::TypedSizedFieldBase; using Base = typename Parent::Base; using Field_p = typename FieldCollection::Field_p; using component_type = T; using Scalar = typename Parent::Scalar; //! Copy constructor TensorField(const TensorField &other) = delete; //! Move constructor TensorField(TensorField &&other) = delete; //! Destructor virtual ~TensorField() = default; //! Copy assignment operator TensorField& operator=(const TensorField &other) = delete; //! Move assignment operator TensorField& operator=(TensorField &&other) = delete; //! accessors inline Dim_t get_order() const; inline Dim_t get_dim() const; //! factory function template friend FieldType& make_field(std::string unique_name, CollectionType & collection, Args&&... args); static TensorField & check_ref(Base & other) { return static_cast(Parent::check_ref(other));} static const TensorField & check_ref(const Base & other) { return static_cast(Parent::check_ref(other));} /** * Pure convenience functions to get a MatrixFieldMap of * appropriate dimensions mapped to this field. You can also * create other types of maps, as long as they have the right * fundamental type (T) and the correct size (nbComponents). */ decltype(auto) get_map(); decltype(auto) get_map() const; protected: //! constructor protected! TensorField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ template class MatrixField: public internal::TypedSizedFieldBase { public: using Parent = internal::TypedSizedFieldBase; using Base = typename Parent::Base; using Field_p = std::unique_ptr>; using component_type = T; //! Copy constructor MatrixField(const MatrixField &other) = delete; //! Move constructor MatrixField(MatrixField &&other) = delete; //! Destructor virtual ~MatrixField() = default; //! Copy assignment operator MatrixField& operator=(const MatrixField &other) = delete; //! Move assignment operator MatrixField& operator=(MatrixField &&other) = delete; //! accessors inline Dim_t get_nb_row() const; inline Dim_t get_nb_col() const; //! factory function template friend FieldType& make_field(std::string unique_name, CollectionType & collection, Args&&... args); static MatrixField & check_ref(Base & other) { return static_cast(Parent::check_ref(other));} static const MatrixField & check_ref(const Base & other) { return static_cast(Parent::check_ref(other));} decltype(auto) get_map(); decltype(auto) get_map() const; protected: //! constructor protected! MatrixField(std::string unique_name, FieldCollection & collection); private: }; /* ---------------------------------------------------------------------- */ //! convenience alias ( template using ScalarField = MatrixField; /* ---------------------------------------------------------------------- */ // Implementations namespace internal { /* ---------------------------------------------------------------------- */ template FieldBase::FieldBase(std::string unique_name, size_t nb_components_, FieldCollection & collection_) :name(unique_name), nb_components(nb_components_), collection(collection_) {} /* ---------------------------------------------------------------------- */ template inline const std::string & FieldBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template inline const FieldCollection & FieldBase:: get_collection() const { return this->collection; } /* ---------------------------------------------------------------------- */ template inline const size_t & FieldBase:: get_nb_components() const { return this->nb_components; } /* ---------------------------------------------------------------------- */ template TypedFieldBase:: TypedFieldBase(std::string unique_name, size_t nb_components, FieldCollection & collection) :Parent(unique_name, nb_components, collection) {} /* ---------------------------------------------------------------------- */ //! return type_id of stored type template const std::type_info & TypedFieldBase:: get_stored_typeid() const { return typeid(T); } /* ---------------------------------------------------------------------- */ template typename TypedFieldBase::EigenMap TypedFieldBase:: eigen() { return EigenMap(this->data(), this->get_nb_components(), this->size()); } + /* ---------------------------------------------------------------------- */ + template + typename TypedFieldBase::EigenVec + TypedFieldBase:: + eigenvec() { + return EigenVec(this->data(), this->get_nb_components() * this->size()); + } + /* ---------------------------------------------------------------------- */ template TypedSizedFieldBase:: TypedSizedFieldBase(std::string unique_name, FieldCollection & collection) :Parent(unique_name, NbComponents, collection){ static_assert ((std::is_arithmetic::value || std::is_same::value), "Use TypedSizedFieldBase for integer, real or complex scalars for T"); static_assert(NbComponents > 0, "Only fields with more than 0 components"); } /* ---------------------------------------------------------------------- */ template void TypedSizedFieldBase:: set_zero() { std::fill(this->values.begin(), this->values.end(), T{}); } /* ---------------------------------------------------------------------- */ template size_t TypedSizedFieldBase:: size() const { if (ArrayStore) { return this->values.size(); } else { return this->values.size()/NbComponents; } } /* ---------------------------------------------------------------------- */ template TypedSizedFieldBase & TypedSizedFieldBase:: check_ref(Base & other) { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::string err ="Cannot create a Reference of requested type " +( "for field '" + other.get_name() + "' of type '" + other.get_stored_typeid().name() + "'"); throw std::runtime_error (err); } //check size compatibility if (NbComponents != other.get_nb_components()) { throw std::runtime_error ("Cannot create a Reference to a field with " + std::to_string(NbComponents) + " components " + "for field '" + other.get_name() + "' with " + std::to_string(other.get_nb_components()) + " components"); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template const TypedSizedFieldBase & TypedSizedFieldBase:: check_ref(const Base & other) { if (typeid(T).hash_code() != other.get_stored_typeid().hash_code()) { std::stringstream err_str{}; err_str << "Cannot create a Reference of requested type " << "for field '" << other.get_name() << "' of type '" << other.get_stored_typeid().name() << "'"; throw std::runtime_error (err_str.str()); } //check size compatibility if (NbComponents != other.get_nb_components()) { throw std::runtime_error ("Cannot create a Reference to a field with " + std::to_string(NbComponents) + " components " + "for field '" + other.get_name() + "' with " + std::to_string(other.get_nb_components()) + " components"); } return static_cast(other); } /* ---------------------------------------------------------------------- */ template typename TypedSizedFieldBase::EigenMap TypedSizedFieldBase:: eigen() { return EigenMap(this->data(), NbComponents, this->size()); } /* ---------------------------------------------------------------------- */ template typename TypedSizedFieldBase::ConstEigenMap TypedSizedFieldBase:: eigen() const { return ConstEigenMap(this->data(), NbComponents, this->size()); } /* ---------------------------------------------------------------------- */ template template Real TypedSizedFieldBase:: inner_product(const TypedSizedFieldBase & other) const { return (this->eigen() * other.eigen()).sum(); } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedSizedFieldBase:: get_ptr_to_entry(const size_t&& index) { static_assert (isArray == ArrayStore, "SFINAE"); return &this->values[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template template T* TypedSizedFieldBase:: get_ptr_to_entry(std::enable_if_t index) { static_assert (noArray != ArrayStore, "SFINAE"); return &this->values[NbComponents*std::move(index)]; } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedSizedFieldBase:: get_ptr_to_entry(const size_t&& index) const { static_assert (isArray == ArrayStore, "SFINAE"); return &this->values[std::move(index)](0, 0); } /* ---------------------------------------------------------------------- */ template template const T* TypedSizedFieldBase:: get_ptr_to_entry(std::enable_if_t index) const { static_assert (noArray != ArrayStore, "SFINAE"); return &this->values[NbComponents*std::move(index)]; } /* ---------------------------------------------------------------------- */ template void TypedSizedFieldBase:: resize(size_t size) { if (ArrayStore) { this->values.resize(size); } else { this->values.resize(size*NbComponents); } } /* ---------------------------------------------------------------------- */ template template void TypedSizedFieldBase:: push_back(const std::enable_if_t & value) { static_assert(isArrayStore == ArrayStore, "SFINAE"); this->values.push_back(value); } /* ---------------------------------------------------------------------- */ template template std::enable_if_t TypedSizedFieldBase:: push_back(const StoredType & value) { static_assert(componentStore != ArrayStore, "SFINAE"); for (Dim_t i = 0; i < NbComponents; ++i) { this->values.push_back(value(i)); } } } // internal /* ---------------------------------------------------------------------- */ //! Factory function, guarantees that only fields get created that are //! properly registered and linked to a collection. template FieldType & make_field(std::string unique_name, FieldCollection & collection, Args&&... args) { auto ptr = std::unique_ptr{ new FieldType(unique_name, collection, args...)}; auto& retref = *ptr; collection.register_field(std::move(ptr)); return retref; } /* ---------------------------------------------------------------------- */ template TensorField:: TensorField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_order() const { return order; } /* ---------------------------------------------------------------------- */ template Dim_t TensorField:: get_dim() const { return dim; } /* ---------------------------------------------------------------------- */ template MatrixField:: MatrixField(std::string unique_name, FieldCollection & collection) :Parent(unique_name, collection) {} /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_col() const { return NbCol; } /* ---------------------------------------------------------------------- */ template Dim_t MatrixField:: get_nb_row() const { return NbRow; } } // muSpectre #include "common/field_map.hh" namespace muSpectre { namespace internal { /* ---------------------------------------------------------------------- */ template struct tensor_map_type { }; template struct tensor_map_type { using type = MatrixFieldMap; }; template struct tensor_map_type { using type = MatrixFieldMap; }; template struct tensor_map_type { using type = T4MatrixFieldMap; }; /* ---------------------------------------------------------------------- */ template struct matrix_map_type { using type = MatrixFieldMap; }; template struct matrix_map_type { using type = ScalarFieldMap; }; } // internal /* ---------------------------------------------------------------------- */ template decltype(auto) TensorField:: get_map() { constexpr bool map_constness{false}; using RawMap_t = typename internal::tensor_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template decltype(auto) TensorField:: get_map() const { constexpr bool map_constness{true}; using RawMap_t = typename internal::tensor_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template decltype(auto) MatrixField:: get_map() { constexpr bool map_constness{false}; using RawMap_t = typename internal::matrix_map_type::type; return RawMap_t(*this); } /* ---------------------------------------------------------------------- */ template decltype(auto) MatrixField:: get_map() const { constexpr bool map_constness{true}; using RawMap_t = typename internal::matrix_map_type::type; return RawMap_t(*this); } } // muSpectre #endif /* FIELD_H */ diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc index e1f6991..59e9b4f 100644 --- a/src/solver/solvers.cc +++ b/src/solver/solvers.cc @@ -1,383 +1,397 @@ /** * file solvers.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief implementation of solver functions * * @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 "solver/solvers.hh" #include "solver/solver_cg.hh" #include "common/iterators.hh" +#include + #include #include namespace muSpectre { template std::vector de_geus (SystemBase & sys, const GradIncrements & delFs, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose) { using Field_t = typename MaterialBase::StrainField_t; auto solver_fields{std::make_unique>()}; solver_fields->initialise(sys.get_resolutions()); // Corresponds to symbol δF or δε auto & incrF{make_field("δF", *solver_fields)}; // Corresponds to symbol ΔF or Δε auto & DeltaF{make_field("ΔF", *solver_fields)}; // field to store the rhs for cg calculations auto & rhs{make_field("rhs", *solver_fields)}; - SolverCG cg(sys.get_resolutions(), - cg_tol, maxiter, verbose-1>0); - cg.initialise(); + // SolverCG cg(sys.get_resolutions(), + // cg_tol, maxiter, verbose-1>0); + // cg.initialise(); + Eigen::ConjugateGradient::Adaptor, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg{}; + cg.setTolerance(cg_tol); + cg.setMaxIterations(maxiter); + auto adaptor = sys.get_adaptor(); + cg.compute(adaptor); if (maxiter == 0) { maxiter = sys.size()*DimM*DimM*10; } size_t count_width{}; const auto form{sys.get_formulation()}; std::string strain_symb{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "de Geus for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "F"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << cg_tol << " maxiter = " << maxiter << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(maxiter))+1; } // initialise F = I or ε = 0 auto & F{sys.get_strain()}; switch (form) { case Formulation::finite_strain: { F.get_map() = Matrices::I2(); break; } case Formulation::small_strain: { F.get_map() = Matrices::I2().Zero(); for (const auto & delF: delFs) { if (!check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown formulation"); break; } // initialise return value std::vector ret_val{}; // initialise materials constexpr bool need_tangent{true}; sys.initialise_materials(need_tangent); Grad_t previous_grad{Grad_t::Zero()}; + Uint cg_counter{0}; for (const auto & delF: delFs) { //incremental loop Real incrNorm{2*newton_tol}, gradNorm{1}; auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () { return incrNorm/gradNorm <= newton_tol; }; Uint newt_iter{0}; for (; (newt_iter < maxiter) && (!convergence_test() || (newt_iter==1)); ++newt_iter) { // obtain material response auto res_tup{sys.evaluate_stress_tangent(F)}; auto & P{std::get<0>(res_tup)}; auto & K{std::get<1>(res_tup)}; auto tangent_effect = [&sys, &K] (const Field_t & dF, Field_t & dP) { sys.directional_stiffness(K, dF, dP); }; if (newt_iter == 0) { DeltaF.get_map() = -(delF-previous_grad); // neg sign because rhs tangent_effect(DeltaF, rhs); - cg.solve(tangent_effect, rhs, incrF); + // cg.solve(tangent_effect, rhs, incrF); + incrF.eigenvec() = cg.solve(rhs.eigenvec()); F.eigen() -= DeltaF.eigen(); } else { rhs.eigen() = -P.eigen(); sys.project(rhs); - cg.solve(tangent_effect, rhs, incrF); + //cg.solve(tangent_effect, rhs, incrF); + incrF.eigenvec() = cg.solve(rhs.eigenvec()); + } + if (cg.info() != Eigen::Success) { + throw SolverError("cg did not converge"); } + cg_counter += cg.iterations(); F.eigen() += incrF.eigen(); incrNorm = incrF.eigen().matrix().norm(); gradNorm = F.eigen().matrix().norm(); if (verbose>0) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << F.get_map().mean() << std::endl; } } } // update previous gradient previous_grad = delF; ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(), convergence_test(), Int(convergence_test()), "message not yet implemented", - newt_iter, cg.get_counter()}); + newt_iter, cg_counter}); //!store history variables here } return ret_val; } template std::vector de_geus (SystemBase & sys, const GradIncrements& delF0, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); // template typename SystemBase::StrainField_t & // de_geus (SystemBase & sys, const GradIncrements& delF0, // const Real cg_tol, const Real newton_tol, Uint maxiter, // Dim_t verbose); template std::vector de_geus (SystemBase & sys, const GradIncrements& delF0, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); /* ---------------------------------------------------------------------- */ template std::vector newton_cg (SystemBase & sys, const GradIncrements & delFs, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose) { using Field_t = typename MaterialBase::StrainField_t; auto solver_fields{std::make_unique>()}; solver_fields->initialise(sys.get_resolutions()); // Corresponds to symbol δF or δε auto & incrF{make_field("δF", *solver_fields)}; // field to store the rhs for cg calculations auto & rhs{make_field("rhs", *solver_fields)}; SolverCG cg(sys.get_resolutions(), cg_tol, maxiter, verbose-1>0); cg.initialise(); if (maxiter == 0) { maxiter = sys.size()*DimM*DimM*10; } size_t count_width{}; const auto form{sys.get_formulation()}; std::string strain_symb{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "Newton-CG for "; switch (form) { case Formulation::small_strain: { strain_symb = "ε"; std::cout << "small"; break; } case Formulation::finite_strain: { strain_symb = "Fy"; std::cout << "finite"; break; } default: throw SolverError("unknown formulation"); break; } std::cout << " strain with" << std::endl << "newton_tol = " << newton_tol << ", cg_tol = " << cg_tol << " maxiter = " << maxiter << " and Δ" << strain_symb << " =" <(tup)}; auto && grad{std::get<1>(tup)}; std::cout << "Step " << counter + 1 << ":" << std::endl << grad << std::endl; } count_width = size_t(std::log10(maxiter))+1; } // initialise F = I or ε = 0 auto & F{sys.get_strain()}; switch (sys.get_formulation()) { case Formulation::finite_strain: { F.get_map() = Matrices::I2(); break; } case Formulation::small_strain: { F.get_map() = Matrices::I2().Zero(); for (const auto & delF: delFs) { if (!check_symmetry(delF)) { throw SolverError("all Δε must be symmetric!"); } } break; } default: throw SolverError("Unknown formulation"); break; } // initialise return value std::vector ret_val{}; // initialise materials constexpr bool need_tangent{true}; sys.initialise_materials(need_tangent); Grad_t previous_grad{Grad_t::Zero()}; for (const auto & delF: delFs) { //incremental loop // apply macroscopic strain increment for (auto && grad: F.get_map()) { grad += delF - previous_grad; } Real incrNorm{2*newton_tol}, gradNorm{1}; auto convergence_test = [&incrNorm, &gradNorm, &newton_tol] () { return incrNorm/gradNorm <= newton_tol; }; Uint newt_iter{0}; for (; newt_iter < maxiter && !convergence_test(); ++newt_iter) { // obtain material response auto res_tup{sys.evaluate_stress_tangent(F)}; auto & P{std::get<0>(res_tup)}; auto & K{std::get<1>(res_tup)}; auto tangent_effect = [&sys, &K] (const Field_t & dF, Field_t & dP) { sys.directional_stiffness(K, dF, dP); }; rhs.eigen() = -P.eigen(); sys.project(rhs); cg.solve(tangent_effect, rhs, incrF); F.eigen() += incrF.eigen(); incrNorm = incrF.eigen().matrix().norm(); gradNorm = F.eigen().matrix().norm(); if (verbose > 0) { std::cout << "at Newton step " << std::setw(count_width) << newt_iter << ", |δ" << strain_symb << "|/|Δ" << strain_symb << "| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose-1>1) { std::cout << "<" << strain_symb << "> =" << std::endl << F.get_map().mean() << std::endl; } } } // update previous gradient previous_grad = delF; ret_val.push_back(OptimizeResult{F.eigen(), sys.get_stress().eigen(), convergence_test(), Int(convergence_test()), "message not yet implemented", newt_iter, cg.get_counter()}); //store history variables here } return ret_val; } template std::vector newton_cg (SystemBase & sys, const GradIncrements& delF0, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); // template typename SystemBase::StrainField_t & // newton_cg (SystemBase & sys, const GradIncrements& delF0, // const Real cg_tol, const Real newton_tol, Uint maxiter, // Dim_t verbose); template std::vector newton_cg (SystemBase & sys, const GradIncrements& delF0, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); /* ---------------------------------------------------------------------- */ bool check_symmetry(const Eigen::Ref& eps, Real rel_tol){ return rel_tol >= (eps-eps.transpose()).matrix().norm()/eps.matrix().norm(); } } // muSpectre diff --git a/src/system/system_base.cc b/src/system/system_base.cc index 66968c6..9ae8c8c 100644 --- a/src/system/system_base.cc +++ b/src/system/system_base.cc @@ -1,249 +1,283 @@ /** * file system_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief Implementation for system base class * * @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 "system/system_base.hh" #include "common/ccoord_operations.hh" #include "common/iterators.hh" #include "common/tensor_algebra.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template SystemBase::SystemBase(Projection_ptr projection_) :resolutions{projection_->get_resolutions()}, pixels(resolutions), lengths{projection_->get_lengths()}, fields{std::make_unique()}, F{make_field("Gradient", *this->fields)}, P{make_field("Piola-Kirchhoff-1", *this->fields)}, projection{std::move(projection_)}, form{projection->get_formulation()} { } /* ---------------------------------------------------------------------- */ template typename SystemBase::Material_t & SystemBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template typename SystemBase::FullResponse_t SystemBase::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->form); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template typename SystemBase::StressField_t & SystemBase::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 SystemBase::SolvVectorOut + SystemBase::directional_stiffness_vec(const SolvVectorIn &delF) { + if (!this->K) { + throw std::runtime_error + ("corrently only implemented for cases where a stiffness matrix " + "exists"); + } + if (delF.size() != this->nb_dof()) { + std::stringstream err{}; + err << "input should be of size ndof = ¶(" << this->resolutions <<") × " + << DimS << "² = "<< this->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); + SolvVectorOut(in_tempref.data(), this->nb_dof()) = delF; + + this->directional_stiffness(this->K.value(), in_tempref, out_tempref); + return SolvVectorOut(out_tempref.data(), this->nb_dof()); + + } + /* ---------------------------------------------------------------------- */ template Eigen::ArrayXXd SystemBase:: directional_stiffness_with_copy (Eigen::Ref delF) { if (!this->K) { throw std::runtime_error ("corrently 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 SystemBase::StressField_t & SystemBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & SystemBase::get_strain() { if (this->initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename SystemBase::StressField_t & SystemBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ template const typename SystemBase::TangentField_t & SystemBase::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 SystemBase::StrainField_t & SystemBase::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 SystemBase::initialise(FFT_PlanFlags flags) { // check that all pixels have been assigned exactly one material this->check_material_coverage(); // resize all global fields (strain, stress, etc) this->fields->initialise(this->resolutions); // initialise the projection and compute the fft plan this->projection->initialise(flags); this->initialised = true; } /* ---------------------------------------------------------------------- */ template void SystemBase::initialise_materials(bool stiffness) { for (auto && mat: this->materials) { mat->initialise(stiffness); } } /* ---------------------------------------------------------------------- */ template typename SystemBase::iterator SystemBase::begin() { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename SystemBase::iterator SystemBase::end() { return this->pixels.end(); } + /* ---------------------------------------------------------------------- */ + template + SystemAdaptor> + SystemBase::get_adaptor() { + return Adaptor(*this); + } + /* ---------------------------------------------------------------------- */ template void SystemBase::check_material_coverage() { auto nb_pixels = CcoordOps::get_size(this->resolutions); std::vector*> assignments(nb_pixels, nullptr); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->resolutions, 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->resolutions, 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 class SystemBase; template class SystemBase; } // muSpectre diff --git a/src/system/system_base.hh b/src/system/system_base.hh index 79a1ede..7451533 100644 --- a/src/system/system_base.hh +++ b/src/system/system_base.hh @@ -1,184 +1,251 @@ /** * file system_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell system with single * projection operator * * @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. */ #ifndef SYSTEM_BASE_H #define SYSTEM_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 "system/system_traits.hh" #include #include #include #include namespace muSpectre { - + template + class SystemAdaptor; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class SystemBase { public: using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; using FieldCollection_t = FieldCollection; using Collection_ptr = std::unique_ptr; using Material_t = MaterialBase; using Material_ptr = std::unique_ptr; using Projection_t = ProjectionBase; using Projection_ptr = std::unique_ptr; using StrainField_t = TensorField; using StressField_t = TensorField; using TangentField_t = TensorField; using FullResponse_t = std::tuple; using iterator = typename CcoordOps::Pixels::iterator; + using SolvVectorIn = Eigen::Ref; + using SolvVectorOut = Eigen::Map; + using Adaptor = SystemAdaptor; //! Default constructor SystemBase() = delete; //! constructor using sizes and resolution SystemBase(Projection_ptr projection); //! Copy constructor SystemBase(const SystemBase &other) = delete; //! Move constructor SystemBase(SystemBase &&other) = default; //! Destructor virtual ~SystemBase() = default; //! Copy assignment operator SystemBase& operator=(const SystemBase &other) = delete; //! Move assignment operator SystemBase& operator=(SystemBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this system */ Material_t & add_material(Material_ptr mat); /** * evaluates all materials */ 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 + */ + SolvVectorOut directional_stiffness_vec(const SolvVectorIn & 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); StrainField_t & get_strain(); const StressField_t & get_stress() const; const TangentField_t & get_tangent(bool create = false); 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); /** * initialise materials (including resetting any history variables) */ void initialise_materials(bool stiffness=false); iterator begin(); iterator end(); size_t size() const {return pixels.size();} const Ccoord & get_resolutions() const {return this->resolutions;} const Rcoord & get_lengths() const {return this->lengths;} /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const { return this->projection->get_formulation();} bool is_initialised() const {return this->initialised;} Eigen::Map get_projection() { return this->projection->get_operator();} + + constexpr static Dim_t get_sdim() {return DimS;}; + + Adaptor get_adaptor(); + Dim_t nb_dof() const {return this->size()*ipow(DimS, 2);}; + + protected: //! make sure that every pixel is assigned to one and only one material void check_material_coverage(); const Ccoord & resolutions; CcoordOps::Pixels pixels; const Rcoord & lengths; Collection_ptr fields; StrainField_t & F; StressField_t & P; //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; std::vector materials{}; Projection_ptr projection; bool initialised{false}; const Formulation form; private: }; + + template + class SystemAdaptor: public Eigen::EigenBase> { + + public: + using Scalar = double; + using RealScalar = double; + using StorageIndex = int; + enum { + ColsAtCompileTime = Eigen::Dynamic, + MaxColsAtCompileTime = Eigen::Dynamic, + RowsAtCompileTime = Eigen::Dynamic, + MaxRowsAtCompileTime = Eigen::Dynamic, + IsRowMajor = false + }; + + SystemAdaptor(System & sys):sys{sys}{} + Eigen::Index rows() const {return this->sys.nb_dof();} + Eigen::Index cols() const {return this->rows();} + + template + Eigen::Product + operator*(const Eigen::MatrixBase& x) const { + return Eigen::Product + (*this, x.derived()); + } + System & sys; + }; + } // muSpectre +// Implementation of SystemAdaptor * Eigen::DenseVector though a specialization of internal::generic_product_impl: +namespace Eigen { + namespace internal { + template + struct generic_product_impl // GEMV stands for matrix-vector + : generic_product_impl_base > + { + typedef typename Product::Scalar Scalar; + template + static void scaleAndAddTo(Dest& dst, const SystemAdaptor& 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).sys.directional_stiffness_vec(rhs); + } + }; + } +} + + #endif /* SYSTEM_BASE_H */ diff --git a/src/system/system_traits.hh b/src/system/system_traits.hh new file mode 100644 index 0000000..bc7ae0f --- /dev/null +++ b/src/system/system_traits.hh @@ -0,0 +1,50 @@ +/** + * file system_traits.hh + * + * @author Till Junge + * + * @date 19 Jan 2018 + * + * @brief Provides traits for Eigen solvers to be able to use Systems + * + * @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. + */ + +#include "common/common.hh" + +#include + +namespace muSpectre { + + template + class SystemAdaptor; + +} // muSpectre + +namespace Eigen { + namespace internal { + using Dim_t = muSpectre::Dim_t; + using Real = muSpectre::Real; + template + struct traits> : + public Eigen::internal::traits > + {}; + } // internal +} // Eigen diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index bb2aba4..20a6485 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,175 +1,175 @@ /** * file test_solver_newton_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief Tests for the standard Newton-Raphson + Conjugate Gradient solver * * @section LICENCE * * 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 "tests.hh" #include "solver/solvers.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain_fast.hh" #include "materials/material_hyper_elastic1.hh" #include "common/iterators.hh" #include "common/ccoord_operations.hh" #include "system/system_factory.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(newton_cg_tests); BOOST_AUTO_TEST_CASE(manual_construction_test) { // constexpr Dim_t dim{twoD}; constexpr Dim_t dim{threeD}; // constexpr Ccoord_t resolutions{3, 3}; // constexpr Rcoord_t lengths{2.3, 2.7}; constexpr Ccoord_t resolutions{5, 5, 5}; constexpr Rcoord_t lengths{5, 5, 5}; auto fft_ptr{std::make_unique>(resolutions, lengths)}; auto proj_ptr{std::make_unique>(std::move(fft_ptr))}; SystemBase sys(std::move(proj_ptr)); using Mat_t = MaterialHyperElastic1; //const Real Young{210e9}, Poisson{.33}; const Real Young{1.0030648180242636}, Poisson{0.29930675909878679}; // const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))}; // const Real mu{Young/(2*(1+Poisson))}; auto& Material_hard = Mat_t::make(sys, "hard", 10*Young, Poisson); auto& Material_soft = Mat_t::make(sys, "soft", Young, Poisson); for (auto && tup: akantu::enumerate(sys)) { auto && pixel = std::get<1>(tup); if (std::get<0>(tup) == 0) { Material_hard.add_pixel(pixel); } else { Material_soft.add_pixel(pixel); } } sys.initialise(); Grad_t delF0; delF0 << 0, 1., 0, 0, 0, 0, 0, 0, 0; constexpr Real cg_tol{1e-8}, newton_tol{1e-5}; constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10}; constexpr bool verbose{false}; GradIncrements grads; grads.push_back(delF0); Eigen::ArrayXXd res1{de_geus(sys, grads, cg_tol, newton_tol, maxiter, verbose)[0].grad}; Eigen::ArrayXXd res2{newton_cg(sys, grads, cg_tol, newton_tol, maxiter, verbose)[0].grad}; BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol); } - // BOOST_AUTO_TEST_CASE(small_strain_patch_test) { - // constexpr Dim_t dim{threeD}; - // using Ccoord = Ccoord_t; - // using Rcoord = Rcoord_t; - // constexpr Ccoord resolutions{CcoordOps::get_cube(11)}; - // constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; - // constexpr Formulation form{Formulation::small_strain}; - - // // number of layers in the hard material - // constexpr Uint nb_lays{1}; - // constexpr Real contrast{2}; - // static_assert(nb_lays < resolutions[0], - // "the number or layers in the hard material must be smaller " - // "than the total number of layers in dimension 0"); - - // auto sys{make_system(resolutions, lengths, form)}; - - // using Mat_t = MaterialHyperElastic1; - // constexpr Real Young{2.}, Poisson{.33}; - // auto material_hard{std::make_unique("hard", contrast*Young, Poisson)}; - // auto material_soft{std::make_unique("soft", Young, Poisson)}; - - // for (const auto & pixel: sys) { - // if (pixel[0] < Dim_t(nb_lays)) { - // material_hard->add_pixel(pixel); - // } else { - // material_soft->add_pixel(pixel); - // } - // } - - // sys.add_material(std::move(material_hard)); - // sys.add_material(std::move(material_soft)); - // sys.initialise(); - - // Grad_t delEps0{Grad_t::Zero()}; - // constexpr Real eps0 = 1.; - // delEps0(0, 1) = delEps0(1, 0) = eps0; - - // constexpr Real cg_tol{1e-8}, newton_tol{1e-5}; - // constexpr Uint maxiter{dim*10}; - // constexpr Dim_t verbose{3}; - - // auto result = de_geus(sys, delEps0, cg_tol, newton_tol, maxiter, verbose); - // if (verbose) { - // std::cout << "result:" << std::endl << result.grad << std::endl; - // std::cout << "mean strain = " << std::endl - // << sys.get_strain().get_map().mean() << std::endl; - // } - - // /** - // * verification of resultant strains: subscript ₕ for hard and ₛ - // * for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast - // * - // * Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ - // * => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ - // * - // * σ is constant across all layers - // * σₕ = σₛ - // * => Eₕ εₕ = Eₛ εₛ - // * => εₕ = 1/k εₛ - // * => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ - // */ - // constexpr Real factor{1/contrast * Real(nb_lays)/resolutions[0] - // + 1.-nb_lays/Real(resolutions[0])}; - // constexpr Real eps_soft{eps0/factor}; - // constexpr Real eps_hard{eps_soft/contrast}; - // if (verbose) { - // std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl; - // std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl; - // } - // Grad_t Eps_hard; Eps_hard << eps_hard, 0, 0, 0; - // Grad_t Eps_soft; Eps_soft << eps_soft, 0, 0, 0; - - // //TODO small strain projection incorrect - // // for (const auto & pixel: sys) { - // // if (pixel[0] < Dim_t(nb_lays)) { - // // BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol); - // // } else { - // // BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol); - // // } - // // } - - - // } + BOOST_AUTO_TEST_CASE(small_strain_patch_test) { + constexpr Dim_t dim{twoD}; + using Ccoord = Ccoord_t; + using Rcoord = Rcoord_t; + constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; + constexpr Rcoord lengths{CcoordOps::get_cube(1.)}; + constexpr Formulation form{Formulation::small_strain}; + + // number of layers in the hard material + constexpr Uint nb_lays{1}; + constexpr Real contrast{2}; + static_assert(nb_lays < resolutions[0], + "the number or layers in the hard material must be smaller " + "than the total number of layers in dimension 0"); + + auto sys{make_system(resolutions, lengths, form)}; + + using Mat_t = MaterialHyperElastic1; + constexpr Real Young{2.}, Poisson{.33}; + auto material_hard{std::make_unique("hard", contrast*Young, Poisson)}; + auto material_soft{std::make_unique("soft", Young, Poisson)}; + + for (const auto & pixel: sys) { + if (pixel[0] < Dim_t(nb_lays)) { + material_hard->add_pixel(pixel); + } else { + material_soft->add_pixel(pixel); + } + } + + sys.add_material(std::move(material_hard)); + sys.add_material(std::move(material_soft)); + sys.initialise(); + + Grad_t delEps0{Grad_t::Zero()}; + constexpr Real eps0 = 1.; + delEps0(0, 1) = delEps0(1, 0) = eps0; + + constexpr Real cg_tol{1e-8}, newton_tol{1e-5}; + constexpr Uint maxiter{dim*10}; + constexpr Dim_t verbose{3}; + + auto result = de_geus(sys, delEps0, cg_tol, newton_tol, maxiter, verbose); + if (verbose) { + std::cout << "result:" << std::endl << result.grad << std::endl; + std::cout << "mean strain = " << std::endl + << sys.get_strain().get_map().mean() << std::endl; + } + + /** + * verification of resultant strains: subscript ₕ for hard and ₛ + * for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast + * + * Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ + * => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ + * + * σ is constant across all layers + * σₕ = σₛ + * => Eₕ εₕ = Eₛ εₛ + * => εₕ = 1/k εₛ + * => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ + */ + constexpr Real factor{1/contrast * Real(nb_lays)/resolutions[0] + + 1.-nb_lays/Real(resolutions[0])}; + constexpr Real eps_soft{eps0/factor}; + constexpr Real eps_hard{eps_soft/contrast}; + if (verbose) { + std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl; + std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl; + } + Grad_t Eps_hard; Eps_hard << eps_hard, 0, 0, 0; + Grad_t Eps_soft; Eps_soft << eps_soft, 0, 0, 0; + + //TODO small strain projection incorrect + // for (const auto & pixel: sys) { + // if (pixel[0] < Dim_t(nb_lays)) { + // BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol); + // } else { + // BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol); + // } + // } + + + } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre