diff --git a/bin/demonstrator2.cc b/bin/demonstrator2.cc index d83e24c..ad70863 100644 --- a/bin/demonstrator2.cc +++ b/bin/demonstrator2.cc @@ -1,92 +1,91 @@ /** * file demonstrator1.cc * * @author Till Junge * * @date 03 Jan 2018 * * @brief larger problem to show off * * @section LICENCE * * 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 #include #include #include "common/common.hh" #include "common/ccoord_operations.hh" #include "system/system_factory.hh" #include "materials/material_hyper_elastic1.hh" #include "solver/solvers.hh" #include "solver/solver_cg.hh" using namespace muSpectre; int main() { banner("demonstrator1", 2018, "Till Junge "); constexpr Dim_t dim{2}; constexpr Formulation form{Formulation::finite_strain}; const Rcoord_t lengths{5.2, 8.3}; const Ccoord_t resolutions{5, 7}; auto system{make_system(resolutions, lengths, form)}; constexpr Real E{1.0030648180242636}; constexpr Real nu{0.29930675909878679}; using Material_t = MaterialHyperElastic1; auto & soft{Material_t::make(system, "soft", E, nu)}; auto & hard{Material_t::make(system, "hard", 10*E, nu)}; int counter{0}; for (const auto && pixel:system) { if (counter < 3) { hard.add_pixel(pixel); counter++; } else { soft.add_pixel(pixel); } } std::cout << counter << " Pixel out of " << system.size() << " are in the hard material" << std::endl; system.initialise(); constexpr Real newton_tol{1e-4}; constexpr Real cg_tol{1e-7}; const size_t maxiter = 100; Grad_t DeltaF{Grad_t::Zero()}; DeltaF(0, 1) = .1; Dim_t verbose {1}; auto start = std::chrono::high_resolution_clock::now(); - GradIncrements grads{DeltaF}; - auto &res = de_geus(system, grads, cg_tol, newton_tol, maxiter, verbose); + auto res = de_geus(system, DeltaF, cg_tol, newton_tol, maxiter, verbose); std::chrono::duration dur = std::chrono::high_resolution_clock::now() - start; std::cout << "Resolution time = " << dur.count() << "s" << std::endl; - std::cout << res.eigen().transpose() << std::endl; + std::cout << res.grad.transpose() << std::endl; return 0; } diff --git a/bin/small_case.cc b/bin/small_case.cc index 34d0dd1..3fa57f7 100644 --- a/bin/small_case.cc +++ b/bin/small_case.cc @@ -1,83 +1,83 @@ /** * file small_case.cc * * @author Till Junge * * @date 12 Jan 2018 * * @brief small case for debugging * * @section LICENCE * * 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 "common/iterators.hh" #include "system/system_factory.hh" #include "materials/material_hyper_elastic1.hh" #include "solver/solvers.hh" #include using namespace muSpectre; int main() { constexpr Dim_t dim{twoD}; Ccoord_t resolution{3, 3}; Rcoord_t lengths{CcoordOps::get_cube(3.)};//{5.2e-9, 8.3e-9, 8.3e-9}; Formulation form{Formulation::finite_strain}; auto rve{make_system(resolution, lengths, form)}; auto & hard{MaterialHyperElastic1::make (rve, "hard", 210e9, .33)}; auto & soft{MaterialHyperElastic1::make (rve, "soft", 70e9, .33)}; for (auto && tup: akantu::enumerate(rve)) { auto & i = std::get<0>(tup); auto & pixel = std::get<1>(tup); if (i < 3) { hard.add_pixel(pixel); } else { soft.add_pixel(pixel); } } rve.initialise(); Real tol{1e-6}; Grad_t Del0{}; Del0 << 0, .1, 0, 0; Dim_t maxiter{31}; Dim_t verbose{3}; - auto & res = de_geus(rve, Del0, tol, tol, maxiter, verbose); - std::cout << res.eigen().transpose() << std::endl; + auto res = de_geus(rve, Del0, tol, tol, maxiter, verbose); + std::cout << res.grad.transpose() << std::endl; return 0; } diff --git a/bin/small_case.py b/bin/small_case.py index 9ff2693..f114c95 100755 --- a/bin/small_case.py +++ b/bin/small_case.py @@ -1,71 +1,73 @@ #!/usr/bin/env python3 """ file small_case.py @author Till Junge @date 12 Jan 2018 @brief small case for debugging @section LICENCE 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(), "src/language_bindings/python")) +sys.path.append(os.path.join(os.getcwd(), "language_bindings/python")) import pyMuSpectre as µ resolution = [3, 3] lengths = [3., 3.] formulation = µ.Formulation.finite_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 i < 3: hard.add_pixel(pixel) else: soft.add_pixel(pixel) print("{}, {}".format(i, tuple(pixel))) rve.initialise() tol = 1e-6 Del0 = np.array([[0, .1], [0, 0]]) maxiter = 31 verbose = 2 -r = µ.solvers.de_geus2d(rve, Del0, tol, tol, maxiter, verbose) -print(r) #print(r.eigen().T) +r = µ.solvers.de_geus(rve, Del0, tol, tol, maxiter, verbose) +print(r.grad.T) +print(r.stress.T) +print(r) diff --git a/language_bindings/python/bind_py_solvers.cc b/language_bindings/python/bind_py_solvers.cc index e7eee98..421001f 100644 --- a/language_bindings/python/bind_py_solvers.cc +++ b/language_bindings/python/bind_py_solvers.cc @@ -1,131 +1,134 @@ /** * file bind_py_solver.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for the muSpectre solvers * * @section LICENCE * * 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 -#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) { - std::stringstream name_stream {}; - name_stream << "newton_cg" << sdim << "d"; + const char name []{"newton_cg"}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; - using grad_vec = Grad_t; + using grad_vec = GradIncrements; - mod.def(name_stream.str().c_str(), + mod.def(name, [](sys & s, const grad & g, Real ct, Real nt, - Uint max, Dim_t verb) -> Dim_t {//typename sys::StrainField_t & { - newton_cg(s, g, ct, nt, max, verb); - return -1; //! temporary dummy return + 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_stream.str().c_str(), + mod.def(name, [](sys & s, const grad_vec & g, Real ct, Real nt, - Uint max, Dim_t verb) -> Dim_t {//typename sys::StrainField_t & { - newton_cg(s, g, ct, nt, max, verb); - return -1; + 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) { - std::stringstream name_stream {}; - name_stream << "de_geus" << sdim << "d"; - + const char name []{"de_geus"}; constexpr Dim_t mdim{sdim}; using sys = SystemBase; using grad = Grad_t; - using grad_vec = Grad_t; + using grad_vec = GradIncrements; - mod.def(name_stream.str().c_str(), + mod.def(name, [](sys & s, const grad & g, Real ct, Real nt, - Uint max, Dim_t verb) -> Dim_t {//typename sys::StrainField_t & { - de_geus(s, g, ct, nt, max, verb); - return -1; + 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_stream.str().c_str(), + mod.def(name, [](sys & s, const grad_vec & g, Real ct, Real nt, - Uint max, Dim_t verb) -> Dim_t {//typename sys::StrainField_t & { - de_geus(s, g, ct, nt, max, verb); - return -1; + 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 bf7daba..755d17e 100644 --- a/src/common/field.hh +++ b/src/common/field.hh @@ -1,778 +1,796 @@ /** * file field.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief header-only implementation of a field for field collections * * @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. */ #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; //! 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(); 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() {return this->get_ptr_to_entry(0);} - inline const T* data() const {return this->get_ptr_to_entry(0);} + 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; 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 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/solver_cg.cc b/src/solver/solver_cg.cc index 9fc2200..6760a02 100644 --- a/src/solver/solver_cg.cc +++ b/src/solver/solver_cg.cc @@ -1,113 +1,125 @@ /** * file solver_cg.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief Implementation of cg 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 "solver/solver_cg.hh" #include "solver/solver_error.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template SolverCG::SolverCG(Ccoord resolutions, Real tol, Uint maxiter, bool verbose) :Parent(resolutions), r_k{make_field("residual r_k", this->collection)}, p_k{make_field("search direction r_k", this->collection)}, Ap_k{make_field("Effect of tangent A*p_k", this->collection)}, - tol{tol}, maxiter{maxiter}, verbose{verbose} + tol{tol}, maxiter{maxiter}, verbose{verbose}, counter{0} {} /* ---------------------------------------------------------------------- */ template void SolverCG::solve(const Fun_t& tangent_effect, const Field_t & rhs, Field_t & x_f) { // Following implementation of algorithm 5.2 in Nocedal's Numerical Optimization (p. 112) auto r = this->r_k.eigen(); auto p = this->p_k.eigen(); auto Ap = this->Ap_k.eigen(); auto x = x_f.eigen(); // initialisation of algo tangent_effect(x_f, this->r_k); r -= rhs.eigen(); p = -r; this->converged = false; Real rdr = (r*r).sum(); Real tol2 = ipow(this->tol,2); size_t count_width{}; // for output formatting in verbose case if (this->verbose) { count_width = size_t(std::log10(this->maxiter))+1; } - for (Uint i = 0; i < this->maxiter && rdr > tol2; ++i) { + for (Uint i = 0; i < this->maxiter && rdr > tol2; ++i, ++this->counter) { tangent_effect(this->p_k, this->Ap_k); Real alpha = rdr/(p*Ap).sum(); x += alpha * p; r += alpha * Ap; Real new_rdr = (r*r).sum(); Real beta = new_rdr/rdr; rdr = new_rdr; if (this->verbose) { std::cout << " at CG step " << std::setw(count_width) << i << ": |r| = " << std::setw(17) << std::sqrt(rdr) << ", cg_tol = " << this->tol << std::endl; } p = -r+beta*p; } if (rdr < tol2) { this->converged=true; } else { throw ConvergenceError("Conjugate gradient has not converged"); } } + /* ---------------------------------------------------------------------- */ + template + void SolverCG::reset_counter() { + this->counter = 0; + } + + /* ---------------------------------------------------------------------- */ + template + Uint SolverCG::get_counter() const { + return this->counter; + } + /* ---------------------------------------------------------------------- */ template typename SolverCG::Tg_req_t SolverCG::get_tangent_req() const { return tangent_requirement; } template class SolverCG; //template class SolverCG; template class SolverCG; } // muSpectre diff --git a/src/solver/solver_cg.hh b/src/solver/solver_cg.hh index 60264cb..26fdc0f 100644 --- a/src/solver/solver_cg.hh +++ b/src/solver/solver_cg.hh @@ -1,96 +1,104 @@ /** * file solver_cg.hh * * @author Till Junge * * @date 20 Dec 2017 * * @brief class for a simple implementation of a conjugate gradient * solver. This follows algorithm 5.2 in Nocedal's Numerical * Optimization (p 112) * * @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. */ #ifndef SOLVER_CG_H #define SOLVER_CG_H #include "solver/solver_base.hh" #include "common/field.hh" #include namespace muSpectre { template class SolverCG: public SolverBase { public: using Parent = SolverBase; using Ccoord = typename Parent::Ccoord; using Tg_req_t = typename Parent::TangentRequirement; // cg only needs to handle fields that look like strain and stress using Field_t = TensorField< typename Parent::Collection_t, Real, secondOrder, DimM>; using Fun_t = std::function; constexpr static Tg_req_t tangent_requirement{Tg_req_t::NeedEffect}; //! Default constructor SolverCG() = delete; //! Constructor with domain resolutions, etc, SolverCG(Ccoord resolutions, Real tol, Uint maxiter=0, bool verbose =false); //! Copy constructor SolverCG(const SolverCG &other) = delete; //! Move constructor SolverCG(SolverCG &&other) = default; //! Destructor virtual ~SolverCG() = default; //! Copy assignment operator SolverCG& operator=(const SolverCG &other) = delete; //! Move assignment operator SolverCG& operator=(SolverCG &&other) = default; //! actual solver void solve(const Fun_t & tangent_effect, const Field_t & rhs, Field_t & x); + //! reset the iteration counter to zero + void reset_counter(); + + //! get the count of how many solve steps have been executed since + //! construction of most recent counter reset + Uint get_counter() const; + protected: Tg_req_t get_tangent_req() const override final; Field_t & r_k; // residual Field_t & p_k; // search direction Field_t & Ap_k; // effect of tangent on search direction const Real tol; const Uint maxiter; const bool verbose; + Uint counter{}; private: }; } // muSpectre #endif /* SOLVER_CG_H */ diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc index 593ee0a..58fc12a 100644 --- a/src/solver/solvers.cc +++ b/src/solver/solvers.cc @@ -1,328 +1,354 @@ /** * file solvers.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief implementation of solver functions * * @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 "solver/solvers.hh" #include "solver/solver_cg.hh" #include "common/iterators.hh" #include #include namespace muSpectre { template - typename SystemBase::StrainField_t & + 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(); if (maxiter == 0) { maxiter = sys.size()*DimM*DimM*10; } size_t count_width{}; const auto form{sys.get_formulation()}; 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: { std::cout << "small"; break; } case Formulation::finite_strain: { 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 ΔF =" <(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(); 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 Real incrNorm{2*newton_tol}, gradNorm{1}; - for (Uint newt_iter{0}; - (newt_iter < maxiter) && ((incrNorm/gradNorm > newton_tol) || + 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 & delF, Field_t & delP) { sys.directional_stiffness(K, delF, delP); }; if (newt_iter == 0) { DeltaF.get_map() = -(delF-previous_grad); // neg sign because rhs tangent_effect(DeltaF, rhs); cg.solve(tangent_effect, rhs, incrF); F.eigen() -= DeltaF.eigen(); } else { 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 << ", |δF|/|ΔF| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose>1) { std::cout << " =" << std::endl << F.get_map().mean() << std::endl; } } } // update previous gradient previous_grad = delF; - //store history variables here + 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 F; + return ret_val; } - template typename SystemBase::StrainField_t & + 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 typename SystemBase::StrainField_t & + 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 & + 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()}; 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: { std::cout << "small"; break; } case Formulation::finite_strain: { 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 ΔF =" <(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(); 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}; - for (Uint newt_iter{0}; - newt_iter < maxiter && incrNorm/gradNorm> newton_tol; + 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 fun = [&sys, &K] (const Field_t & delF, Field_t & delP) { sys.directional_stiffness(K, delF, delP); }; rhs.eigen() = -P.eigen(); sys.project(rhs); cg.solve(fun, 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 << ", |δF|/|ΔF| = " << std::setw(17) << incrNorm/gradNorm << ", tol = " << newton_tol << std::endl; if (verbose>1) { std::cout << " =" << 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 F; + return ret_val; } - template typename SystemBase::StrainField_t & + 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 typename SystemBase::StrainField_t & + template std::vector newton_cg (SystemBase & sys, const GradIncrements& delF0, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); } // muSpectre diff --git a/src/solver/solvers.hh b/src/solver/solvers.hh index 3faeaef..044038b 100644 --- a/src/solver/solvers.hh +++ b/src/solver/solvers.hh @@ -1,81 +1,108 @@ /** * file solvers.hh * * @author Till Junge * * @date 20 Dec 2017 * * @brief Free functions for solving * * @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. */ #ifndef SOLVERS_H #define SOLVERS_H #include "solver/solver_cg.hh" +#include + +#include +#include + namespace muSpectre { + /** + * emulates scipy.optimize.OptimizeResult + */ + struct OptimizeResult + { + //! Strain ε or Gradient F at solution + Eigen::ArrayXXd grad; + //! Cauchy stress σ or first Piola-Kirchhoff stress P at solution + Eigen::ArrayXXd stress; + //! whether or not the solver exited successfully + bool success; + //! Termination status of the optimizer. Its value depends on the + //! underlying solver. Refer to message for details. + Int status; + //! Description of the cause of the termination. + std::string message; + //! number of iterations + Uint nb_it; + //! number of system evaluations + Uint nb_fev; + }; + template using Grad_t = Matrices::Tens2_t; template using GradIncrements = std::vector, Eigen::aligned_allocator>>; /* ---------------------------------------------------------------------- */ template - typename SystemBase::StrainField_t & + std::vector newton_cg (SystemBase & sys, const GradIncrements & delF0, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0); /* ---------------------------------------------------------------------- */ template - inline typename SystemBase::StrainField_t & + inline OptimizeResult newton_cg (SystemBase & sys, const Grad_t & delF0, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0){ return newton_cg(sys, GradIncrements{delF0}, - cg_tol, newton_tol, maxiter, verbose); + cg_tol, newton_tol, maxiter, verbose)[0]; } /* ---------------------------------------------------------------------- */ template - typename SystemBase::StrainField_t & + std::vector de_geus (SystemBase & sys, const GradIncrements & delF0, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0); /* ---------------------------------------------------------------------- */ template - inline typename SystemBase::StrainField_t & + OptimizeResult de_geus (SystemBase & sys, const Grad_t & delF0, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0){ return de_geus(sys, GradIncrements{delF0}, - cg_tol, newton_tol, maxiter, verbose); + cg_tol, newton_tol, maxiter, verbose)[0]; } } // muSpectre #endif /* SOLVERS_H */ diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index 83c2f95..777657e 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,173 +1,174 @@ /** * 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).eigen()}; + 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).eigen()}; + 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{twoD}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; constexpr Ccoord resolutions{3, 3}; constexpr Rcoord lengths{3, 3}; 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{2e9}, 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; constexpr Real eps0 = 1.; delEps0 << eps0, 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 Dim_t verbose{2}; - auto & result = newton_cg(sys, delEps0, cg_tol, newton_tol, maxiter, verbose); + auto result = newton_cg(sys, delEps0, cg_tol, newton_tol, maxiter, verbose); if (verbose) { - std::cout << "result:" << std::endl << result.eigen() << std::endl; - std::cout << "mean strain = " << std::endl << result.get_map().mean() << std::endl; + 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; for (const auto & pixel: sys) { if (pixel[0] < Dim_t(nb_lays)) { - BOOST_CHECK_LE((Eps_hard-result.get_map()[pixel]).norm(), tol); + BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol); } else { - BOOST_CHECK_LE((Eps_soft-result.get_map()[pixel]).norm(), tol); + BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol); } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre