diff --git a/bin/demonstrator1.cc b/bin/demonstrator1.cc index d2e4de6..1dd4ff2 100644 --- a/bin/demonstrator1.cc +++ b/bin/demonstrator1.cc @@ -1,139 +1,141 @@ /** * file demonstrator1.cc * * @author Till Junge * * @date 03 Jan 2018 * * @brief larger problem to show off * * @section LICENCE * * Copyright (C) 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 "external/cxxopts.hpp" #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 opt_ptr = std::unique_ptr; opt_ptr parse_args(int argc, char **argv) { opt_ptr options = std::make_unique(argv[0], "Tests MPI fft scalability"); try { options->add_options() ("0,N0", "number of rows", cxxopts::value(), "N0") ("h,help", "print help") ("positional", "Positional arguments: these are the arguments that are entered " "without an option", cxxopts::value>()); options->parse_positional(std::vector{"N0", "positional"}); options->parse(argc, argv); if (options->count("help")) { std::cout << options->help({"", "Group"}) << std::endl; exit(0); } if (options->count("N0") != 1 ) { throw cxxopts::OptionException("Parameter N0 missing"); } else if ((*options)["N0"].as()%2 != 1) { throw cxxopts::OptionException("N0 must be odd"); } else if (options->count("positional") > 0) { throw cxxopts::OptionException("There are too many positional arguments"); } } catch (const cxxopts::OptionException & e) { std::cout << "Error parsing options: " << e.what() << std::endl; exit(1); } return options; } using namespace muSpectre; int main(int argc, char *argv[]) { auto options{parse_args(argc, argv)}; auto & opt{*options}; const Dim_t size{opt["N0"].as()}; constexpr Real fsize{1.}; constexpr Dim_t dim{3}; const Dim_t nb_dofs{ipow(size, dim)*ipow(dim, 2)}; std::cout << "Number of dofs: " << nb_dofs << std::endl; + constexpr Formulation form{Formulation::finite_strain}; + const Rcoord_t lengths{CcoordOps::get_cube(fsize)}; const Ccoord_t resolutions{CcoordOps::get_cube(size)}; auto system{make_system(resolutions, lengths)}; constexpr Real E{1.0030648180242636}; constexpr Real nu{0.29930675909878679}; using Material_t = MaterialHyperElastic1; auto Material_soft{std::make_unique("soft", E, nu)}; auto Material_hard{std::make_unique("hard", 10*E, nu)}; int counter{0}; for (const auto && pixel:system) { int sum = 0; for (Dim_t i = 0; i < dim; ++i) { sum += pixel[i]*2 / resolutions[i]; } if (sum == 0) { Material_hard->add_pixel(pixel); counter ++; } else { Material_soft->add_pixel(pixel); } } std::cout << counter << " Pixel out of " << system.size() << " are in the hard material" << std::endl; system.add_material(std::move(Material_soft)); system.add_material(std::move(Material_hard)); system.initialise(FFT_PlanFlags::measure); constexpr Real newton_tol{1e-4}; constexpr Real cg_tol{1e-7}; const size_t maxiter = nb_dofs; 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}; - de_geus(system, grads, cg_tol, newton_tol, maxiter, verbose); + de_geus(system, grads, form, 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; return 0; } diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc index 8dfa0bf..b355061 100644 --- a/src/solver/solvers.cc +++ b/src/solver/solvers.cc @@ -1,271 +1,301 @@ /** * file solvers.cc * * @author Till Junge * * @date 20 Dec 2017 * * @brief implementation of solver functions * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "solver/solvers.hh" #include "solver/solver_cg.hh" #include "common/iterators.hh" #include #include namespace muSpectre { template typename SystemBase::StrainField_t & de_geus (SystemBase & sys, const GradIncrements & delFs, + Formulation form, 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{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "Algo 5.2 with newton_tol = " << newton_tol << ", cg_tol = " << cg_tol << " maxiter = " << maxiter << " and ΔF =" <(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 + // initialise F = I or ε = 0 auto & F{sys.get_strain()}; - F.get_map() = Matrices::I2(); + 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 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) || (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 } return F; } template typename SystemBase::StrainField_t & de_geus (SystemBase & sys, const GradIncrements& delF0, + Formulation form, 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 & de_geus (SystemBase & sys, const GradIncrements& delF0, + Formulation form, const Real cg_tol, const Real newton_tol, Uint maxiter, Dim_t verbose); /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & newton_cg (SystemBase & sys, const GradIncrements & delFs, + Formulation form, 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{}; if (verbose > 0) { //setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111) std::cout << "Algo 5.2 with newton_tol = " << newton_tol << ", cg_tol = " << cg_tol << " maxiter = " << maxiter << " and ΔF =" <(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 + // initialise F = I or ε = 0 auto & F{sys.get_strain()}; - F.get_map() = Matrices::I2(); + 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 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; ++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; //store history variables here } return F; } template typename SystemBase::StrainField_t & newton_cg (SystemBase & sys, const GradIncrements& delF0, + Formulation form, 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 & newton_cg (SystemBase & sys, const GradIncrements& delF0, + Formulation form, 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 df5ab90..1e3fee0 100644 --- a/src/solver/solvers.hh +++ b/src/solver/solvers.hh @@ -1,81 +1,85 @@ /** * file solvers.hh * * @author Till Junge * * @date 20 Dec 2017 * * @brief Free functions for solving * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef SOLVERS_H #define SOLVERS_H #include "solver/solver_cg.hh" namespace muSpectre { template using Grad_t = Matrices::Tens2_t; template using GradIncrements = std::vector, Eigen::aligned_allocator>>; /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & newton_cg (SystemBase & sys, const GradIncrements & delF0, + Formulation form, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0); /* ---------------------------------------------------------------------- */ template inline typename SystemBase::StrainField_t & newton_cg (SystemBase & sys, const Grad_t & delF0, + Formulation form, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0){ - return newton_cg(sys, GradIncrements{delF0}, + return newton_cg(sys, GradIncrements{delF0}, form, cg_tol, newton_tol, maxiter, verbose); } /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & de_geus (SystemBase & sys, const GradIncrements & delF0, + Formulation form, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0); /* ---------------------------------------------------------------------- */ template inline typename SystemBase::StrainField_t & de_geus (SystemBase & sys, const Grad_t & delF0, + Formulation form, const Real cg_tol, const Real newton_tol, const Uint maxiter=0, Dim_t verbose = 0){ - return de_geus(sys, GradIncrements{delF0}, + return de_geus(sys, GradIncrements{delF0}, form, cg_tol, newton_tol, maxiter, verbose); } } // muSpectre #endif /* SOLVERS_H */ diff --git a/src/system/system_base.cc b/src/system/system_base.cc index a901ae4..97dd2c3 100644 --- a/src/system/system_base.cc +++ b/src/system/system_base.cc @@ -1,199 +1,202 @@ /** * file system_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief Implementation for system base class * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "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_) + SystemBase::SystemBase(Projection_ptr projection_, + Formulation form) :resolutions{projection_->get_resolutions()}, pixels(resolutions), lengths{projection_->get_lengths()}, fields{}, F{make_field("Gradient", this->fields)}, P{make_field("Piola-Kirchhoff-1", this->fields)}, - projection{std::move(projection_)} + projection{std::move(projection_)}, + form{form} { } /* ---------------------------------------------------------------------- */ template void SystemBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); } /* ---------------------------------------------------------------------- */ template typename SystemBase::FullResponse_t SystemBase::evaluate_stress_tangent(StrainField_t & grad) { if (this->is_initialised == false) { this->initialise(); } //! High level compatibility checks if (grad.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } if (!this->K) { this->K = make_field("Tangent Stiffness", this->fields); } for (auto & mat: this->materials) { - mat->compute_stresses_tangent(grad, this->P, this->K.value(), form); + 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::StressField_t & SystemBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename SystemBase::StrainField_t & SystemBase::get_strain() { if (this->is_initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename SystemBase::StressField_t & SystemBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ 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->is_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 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 0919846..8be827d 100644 --- a/src/system/system_base.hh +++ b/src/system/system_base.hh @@ -1,155 +1,160 @@ /** * file system_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell system with single * projection operator * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ +#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 #include #include #include - -#ifndef SYSTEM_BASE_H -#define SYSTEM_BASE_H - namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class SystemBase { public: - constexpr static Formulation form{Formulation::finite_strain}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; using FieldCollection_t = FieldCollection; - using Material_ptr = std::unique_ptr>; + using Material_t = MaterialBase; + using Material_ptr = std::unique_ptr; using Projection_ptr = std::unique_ptr>; - using StrainField_t = TensorField; - using StressField_t = TensorField; - using TangentField_t = TensorField; - using FullResponse_t = std::tuple; + using StrainField_t = + TensorField; + using StressField_t = + TensorField; + using TangentField_t = + TensorField; + using FullResponse_t = + std::tuple; using iterator = typename CcoordOps::Pixels::iterator; //! Default constructor SystemBase() = delete; //! constructor using sizes and resolution - SystemBase(Projection_ptr projection); + SystemBase(Projection_ptr projection, + Formulation form=Formulation::finite_strain); //! Copy constructor SystemBase(const SystemBase &other) = delete; //! Move constructor SystemBase(SystemBase &&other) noexcept = default; //! Destructor virtual ~SystemBase() noexcept = 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 */ void add_material(Material_ptr mat); /** * evaluates all materials */ FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * 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; /** * 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;} 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; FieldCollection_t fields; StrainField_t & F; StressField_t & P; - //! Tangent field might not even be required; so this is a - //! pointer instead of a ref + //! 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; + Formulation form; bool is_initialised{false}; private: }; } // muSpectre #endif /* SYSTEM_BASE_H */ diff --git a/src/system/system_factory.hh b/src/system/system_factory.hh index 285523a..46f8e25 100644 --- a/src/system/system_factory.hh +++ b/src/system/system_factory.hh @@ -1,80 +1,82 @@ /** * file system_factory.hh * * @author Till Junge * * @date 15 Dec 2017 * * @brief System factories to help create systems with ease * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef SYSTEM_FACTORY_H #define SYSTEM_FACTORY_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "system/system_base.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/projection_finite_strain.hh" #include "fft/fftw_engine.hh" #include namespace muSpectre { /** * Create a unique ptr to a Projection operator (with appropriate * FFT_engine) to be used in a system constructor */ template , typename Projection=ProjectionFiniteStrainFast> inline std::unique_ptr system_input(Ccoord_t resolutions, Rcoord_t lengths=CcoordOps::get_cube(1.)) { auto fft_ptr{std::make_unique(resolutions, lengths)}; return std::make_unique(std::move(fft_ptr)); } /** * convenience function to create a system (avoids having to build * and move the chain of unique_ptrs */ - template , typename FFT_Engine=FFTW_Engine, typename Projection=ProjectionFiniteStrainFast> inline System make_system(Ccoord_t resolutions, - Rcoord_t lengths=CcoordOps::get_cube(1.)) { + Rcoord_t lengths=CcoordOps::get_cube(1.), + Formulation form=Formulation::finite_strain) { return System{ std::move(system_input - (resolutions, lengths))}; + (resolutions, lengths)), + form}; } } // muSpectre #endif /* SYSTEM_FACTORY_H */ diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc index 2f51d02..448205c 100644 --- a/tests/test_solver_newton_cg.cc +++ b/tests/test_solver_newton_cg.cc @@ -1,90 +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 (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "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 = std::make_unique("hard", 10*Young, Poisson); auto Material_soft = std::make_unique("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.add_material(std::move(Material_hard)); sys.add_material(std::move(Material_soft)); 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}; + constexpr Formulation form{Formulation::finite_strain}; 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, form, cg_tol, newton_tol, maxiter, verbose).eigen()}; - Eigen::ArrayXXd res2{newton_cg(sys, grads, cg_tol, newton_tol, maxiter, verbose).eigen()}; + Eigen::ArrayXXd res2{newton_cg(sys, grads, form, cg_tol, newton_tol, maxiter, verbose).eigen()}; 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{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; + 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 bool verbose{false}; + + auto & result = newton_cg(sys, delEps0, form, 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; + } + + /** + * 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); + } else { + BOOST_CHECK_LE((Eps_soft-result.get_map()[pixel]).norm(), tol); + } + } + + + } + BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_system_base.cc b/tests/test_system_base.cc index b714898..fc7eff0 100644 --- a/tests/test_system_base.cc +++ b/tests/test_system_base.cc @@ -1,172 +1,195 @@ /** * file test_system_base.cc * * @author Till Junge * * @date 14 Dec 2017 * * @brief Tests for the basic system class * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include "tests.hh" #include "common/common.hh" #include "common/iterators.hh" #include "common/field_map.hh" #include "common/test_goodies.hh" #include "system/system_factory.hh" #include "materials/material_hyper_elastic1.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(system_base); template struct Sizes { }; template<> struct Sizes { constexpr static Dim_t sdim{twoD}; constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8};} }; template<> struct Sizes { constexpr static Dim_t sdim{threeD}; constexpr static Ccoord_t get_resolution() { return Ccoord_t{3, 5, 7};} constexpr static Rcoord_t get_lengths() { return Rcoord_t{3.4, 5.8, 6.7};} }; - template + template struct SystemBaseFixture: SystemBase { constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; + constexpr static Formulation formulation{form}; SystemBaseFixture() :SystemBase{ std::move(system_input(Sizes::get_resolution(), - Sizes::get_lengths()))} {} + Sizes::get_lengths())), + form} {} }; - using fixlist = boost::mpl::list, - SystemBaseFixture>; + using fixlist = boost::mpl::list, + SystemBaseFixture, + SystemBaseFixture, + SystemBaseFixture>; BOOST_AUTO_TEST_CASE(manual_construction) { constexpr Dim_t dim{twoD}; Ccoord_t resolutions{3, 3}; Rcoord_t lengths{2.3, 2.7}; auto fft_ptr{std::make_unique>(resolutions, lengths)}; auto proj_ptr{std::make_unique>(std::move(fft_ptr))}; SystemBase sys{std::move(proj_ptr)}; auto sys2{make_system(resolutions, lengths)}; auto sys2b{std::move(sys2)}; BOOST_CHECK_EQUAL(sys2b.size(), sys.size()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_THROW(fix::check_material_coverage(), std::runtime_error); BOOST_CHECK_THROW(fix::initialise(), std::runtime_error); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(add_material_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; using Material_t = MaterialHyperElastic1; auto Material_hard = std::make_unique("hard", 210e9, .33); BOOST_CHECK_NO_THROW(fix::add_material(std::move(Material_hard))); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(simple_evaluation_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; + constexpr Formulation form{fix::formulation}; using Mat_t = MaterialHyperElastic1; const Real Young{210e9}, Poisson{.33}; const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))}; const Real mu{Young/(2*(1+Poisson))}; auto Material_hard = std::make_unique("hard", Young, Poisson); for (auto && pixel: *this) { Material_hard->add_pixel(pixel); } fix::add_material(std::move(Material_hard)); auto & F = fix::get_strain(); auto F_map = F.get_map(); + // finite strain formulation expects the deformation gradient F, + // while small strain expects infinitesimal strain ε for (auto grad: F_map) { - grad = grad.Identity(); + switch (form) { + case Formulation::finite_strain: { + grad = grad.Identity(); + break; + } + case Formulation::small_strain: { + grad = grad.Zero(); + break; + } + default: + BOOST_CHECK(false); + break; + } } fix::initialise_materials(); auto res_tup{fix::evaluate_stress_tangent(F)}; auto stress{std::get<0>(res_tup).get_map()}; auto tangent{std::get<1>(res_tup).get_map()}; auto tup = testGoodies::objective_hooke_explicit (lambda, mu, Matrices::I2()); auto P_ref = std::get<0>(tup); for (auto mat: stress) { Real norm = (mat - P_ref).norm(); BOOST_CHECK_EQUAL(norm, 0.); } auto tan_ref = std::get<1>(tup); for (const auto tan: tangent) { Real norm = (tan - tan_ref).norm(); BOOST_CHECK_EQUAL(norm, 0.); } } BOOST_FIXTURE_TEST_CASE_TEMPLATE(evaluation_test, fix, fixlist, fix) { constexpr Dim_t dim{fix::sdim}; using Mat_t = MaterialHyperElastic1; auto Material_hard = std::make_unique("hard", 210e9, .33); auto Material_soft = std::make_unique("soft", 70e9, .3); for (auto && cnt_pixel: akantu::enumerate(*this)) { auto counter = std::get<0>(cnt_pixel); auto && pixel = std::get<1>(cnt_pixel); if (counter < 5) { Material_hard->add_pixel(pixel); } else { Material_soft->add_pixel(pixel); } } fix::add_material(std::move(Material_hard)); fix::add_material(std::move(Material_soft)); auto & F = fix::get_strain(); fix::initialise_materials(); fix::evaluate_stress_tangent(F); fix::evaluate_stress_tangent(F); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre