diff --git a/build/.forgit b/build/.forgit deleted file mode 100644 index e69de29..0000000 diff --git a/language_bindings/python/bind_py_cell.cc b/language_bindings/python/bind_py_cell.cc index 76d1c5b..1cbea18 100644 --- a/language_bindings/python/bind_py_cell.cc +++ b/language_bindings/python/bind_py_cell.cc @@ -1,233 +1,258 @@ /** * @file bind_py_cell.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief Python bindings for the cell factory function * * 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/ccoord_operations.hh" #include "cell/cell_factory.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include #include #include "pybind11/eigen.h" #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * cell factory for specific FFT engine */ #ifdef WITH_MPI template void add_parallel_cell_factory_helper(py::module & mod, const char *name) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def (name, [](Ccoord res, Rcoord lens, Formulation form, size_t comm) { return make_parallel_cell , FFTEngine> (std::move(res), std::move(lens), std::move(form), std::move(Communicator(MPI_Comm(comm)))); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain, "communicator"_a=size_t(MPI_COMM_SELF)); } #endif /** * the cell factory is only bound for default template parameters */ template void add_cell_factory_helper(py::module & mod) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; mod.def ("CellFactory", [](Ccoord res, Rcoord lens, Formulation form) { return make_cell(std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain); - mod.def + mod.def ("CellFactorySplit", [](Ccoord res, Rcoord lens, Formulation form) { - return make_cell_ptr>(std::move(res), std::move(lens), std::move(form)); + //return make_cell_ptr>(std::move(res), std::move(lens), std::move(form)); + return make_cell_split(std::move(res), std::move(lens), std::move(form)); }, "resolutions"_a, "lengths"_a=CcoordOps::get_cube(1.), "formulation"_a=Formulation::finite_strain); - /*mod.def - ("CellFactorySplit", - [](Ccoord res, Rcoord lens, Formulation form) { - CellSplit splitted_cell = make_cell_split(std::move(res), std::move(lens), std::move(form)); - }, - "resolutions"_a, - "lengths"_a=CcoordOps::get_cube(1.), - "formulation"_a=Formulation::finite_strain);*/ #ifdef WITH_FFTWMPI add_parallel_cell_factory_helper>( mod, "FFTWMPICellFactory"); #endif #ifdef WITH_PFFT add_parallel_cell_factory_helper>( mod, "PFFTCellFactory"); #endif + } void add_cell_factory(py::module & mod) { add_cell_factory_helper(mod); add_cell_factory_helper(mod); } /** * CellBase for which the material and spatial dimension are identical */ template void add_cell_base_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "CellBase" << dim << 'd'; const std::string name = name_stream.str(); using sys_t = CellBase; py::class_(mod, name.c_str()) .def("__len__", &sys_t::size) .def("__iter__", [](sys_t & s) { return py::make_iterator(s.begin(), s.end()); }) .def("initialise", &sys_t::initialise, "flags"_a=FFT_PlanFlags::estimate) .def("directional_stiffness", [](sys_t& cell, py::EigenDRef& v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; constexpr bool create_tangent{true}; auto & K = cell.get_tangent(create_tangent); auto & input = cell.get_managed_field(in_name); auto & output = cell.get_managed_field(out_name); input.eigen() = v; cell.directional_stiffness(K, input, output); return output.eigen(); }, "δF"_a) .def("project", [](sys_t& cell, py::EigenDRef& v) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } if (!cell.is_initialised()) { cell.initialise(); } const std::string in_name{"temp input for projection"}; auto & input = cell.get_managed_field(in_name); input.eigen() = v; cell.project(input); return input.eigen(); }, "field"_a) .def("get_strain",[](sys_t & s) { return Eigen::ArrayXXd(s.get_strain().eigen()); }) .def("get_stress",[](sys_t & s) { return Eigen::ArrayXXd(s.get_stress().eigen()); }) .def("size", &sys_t::size) .def("evaluate_stress_tangent", [](sys_t& cell, py::EigenDRef& v ) { if ((size_t(v.cols()) != cell.size() || size_t(v.rows()) != dim*dim)) { std::stringstream err{}; err << "need array of shape (" << dim*dim << ", " << cell.size() << ") but got (" << v.rows() << ", " << v.cols() << ")."; throw std::runtime_error(err.str()); } auto & strain{cell.get_strain()}; strain.eigen() = v; cell.evaluate_stress_tangent(strain); }, "strain"_a) .def("get_projection", &sys_t::get_projection) .def("get_subdomain_resolutions", &sys_t::get_subdomain_resolutions) .def("get_subdomain_locations", &sys_t::get_subdomain_locations) .def("get_domain_resolutions", &sys_t::get_domain_resolutions) .def("get_domain_lengths", &sys_t::get_domain_resolutions); } +template +void add_cell_split_helper(py::module & mod) { + std::stringstream name_stream{}; + name_stream<< "CellSplit" << dim << 'd'; + const auto name {name_stream.str()}; + using sys_split_t = CellSplit; + using sys_t = CellBase; + using Mat_t = MaterialBase; + py::class_(mod, name.c_str()) + .def("make_precipitate", + [](sys_split_t& cell, Mat_t & mat, std::vector> precipitate_vertices){ + cell.make_automatic_precipitate(precipitate_vertices, mat); + }, + "material"_a, "vertices"_a + ) + .def("complete_material_assignment", + [](sys_split_t& cell, Mat_t & mat){ + cell.complete_material_assignment(mat);}, + "material"_a); +} + void add_cell_base(py::module & mod) { py::class_(mod, "Cell"); add_cell_base_helper (mod); add_cell_base_helper(mod); } + +void add_cell_split(py::module & mod) { + add_cell_split_helper (mod); + add_cell_split_helper (mod); +} + + void add_cell(py::module & mod) { add_cell_factory(mod); auto cell{mod.def_submodule("cell")}; cell.doc() = "bindings for cells and cell factories"; - + cell.def("scale_by_2", [](py::EigenDRef& v) { v *= 2; }); add_cell_base(cell); + add_cell_split(cell); } + + diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh index 0a2be3c..384554b 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,524 +1,523 @@ /** * @file cell_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell cell with single * projection operator * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef CELL_BASE_H #define CELL_BASE_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include #include #include #include #include namespace muSpectre { /** * Cell adaptors implement the matrix-vector multiplication and * allow the system to be used like a sparse matrix in * conjugate-gradient-type solvers */ template class CellAdaptor; /** * Base class for cells that is not templated and therefore can be * in solvers that see cells as runtime-polymorphic objects. This * allows the use of standard * (i.e. spectral-method-implementation-agnostic) solvers, as for * instance the scipy solvers */ class Cell { public: //! sparse matrix emulation using Adaptor = CellAdaptor; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = Eigen::Matrix; //! dynamic matrix type for setting strains using Matrix_t = Eigen::Matrix; //! ref to constant vector using ConstVector_ref = Eigen::Map; //! output vector reference for solvers using Vector_ref = Eigen::Map; //! Default constructor Cell() = default; //! Copy constructor Cell(const Cell &other) = default; //! Move constructor Cell(Cell &&other) = default; //! Destructor virtual ~Cell() = default; //! Copy assignment operator Cell& operator=(const Cell &other) = default; //! Move assignment operator Cell& operator=(Cell &&other) = default; //! for handling double initialisations right bool is_initialised() const {return this->initialised;} //! returns the number of degrees of freedom in the cell virtual Dim_t get_nb_dof() const = 0; //! return the communicator object virtual const Communicator & get_communicator() const = 0; /** * formulation is hard set by the choice of the projection class */ virtual const Formulation & get_formulation() const = 0; /** * returns the material dimension of the problem */ virtual Dim_t get_material_dim() const = 0; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ virtual std::array get_strain_shape() const = 0; /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() = 0; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const = 0; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() = 0; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() = 0; /** * applies the projection operator in-place on the input vector */ virtual void apply_projection(Eigen::Ref vec) = 0; /** * freezes all the history variables of the materials */ virtual void save_history_variables() = 0; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) = 0; /** * set uniform strain (typically used to initialise problems */ virtual void set_uniform_strain(const Eigen::Ref &) = 0; //! get a sparse matrix view on the cell virtual Adaptor get_adaptor() = 0; protected: bool initialised{false}; //!< to handle double initialisation right private: }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase: public Cell { public: using Parent = Cell; using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< physical coordinates type //! global field collection using FieldCollection_t = GlobalFieldCollection; //! the collection is handled in a `std::unique_ptr` using Collection_ptr = std::unique_ptr; //! polymorphic base material type using Material_t = MaterialBase; //! materials handled through `std::unique_ptr`s using Material_ptr = std::unique_ptr; //! polymorphic base projection type using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; //! expected type for strain fields using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; //! iterator type over all cell pixel's using iterator = typename CcoordOps::Pixels::iterator; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = typename Parent::Vector_t; //! ref to constant vector using ConstVector_ref = typename Parent::ConstVector_ref; //! output vector reference for solvers using Vector_ref = typename Parent::Vector_ref; //! sparse matrix emulation using Adaptor = Parent::Adaptor; //! Default constructor CellBase() = delete; //! constructor using sizes and resolution CellBase(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::no); //! Copy constructor CellBase(const CellBase &other) = delete; //! Move constructor CellBase(CellBase &&other); //! Destructor virtual ~CellBase() = default; //! Copy assignment operator CellBase& operator=(const CellBase &other) = delete; //! Move assignment operator CellBase& operator=(CellBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this cell */ Material_t & add_material(Material_ptr mat); /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() override; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const override; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() override; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() override; - /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc. (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) override; //! return the template param DimM (required for polymorphic use of `Cell` Dim_t get_material_dim() const override final {return DimM;} /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const override final; /** * applies the projection operator in-place on the input vector */ void apply_projection(Eigen::Ref vec) override final; /** * set uniform strain (typically used to initialise problems */ void set_uniform_strain(const Eigen::Ref &) override; /** * evaluate all materials */ virtual FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF or G:K:δε) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * vectorized version for eigen solvers, no copy, but only works * when fields have ArrayStore=false */ Vector_ref directional_stiffness_vec(const Eigen::Ref & delF); /** * Evaluate directional stiffness into a temporary array and * return a copy. This is a costly and wasteful interface to * directional_stiffness and should only be used for debugging or * in the python interface */ Eigen::ArrayXXd directional_stiffness_with_copy (Eigen::Ref delF); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); //! returns a ref to the cell's strain field StrainField_t & get_strain(); //! returns a ref to the cell's stress field const StressField_t & get_stress() const; //! returns a ref to the cell's tangent stiffness field const TangentField_t & get_tangent(bool create = false); //! returns a ref to a temporary field managed by the cell StrainField_t & get_managed_field(std::string unique_name); /** * general initialisation; initialises the projection and * fft_engine (i.e. infrastructure) but not the materials. These * need to be initialised separately */ void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, this function * calls this update for each material in the cell */ void save_history_variables() override final; iterator begin(); //!< iterator to the first pixel iterator end(); //!< iterator past the last pixel //! number of pixels in the cell size_t size() const {return pixels.size();} //! return the subdomain resolutions of the cell const Ccoord & get_subdomain_resolutions() const { return this->subdomain_resolutions;} //! return the subdomain locations of the cell const Ccoord & get_subdomain_locations() const { return this->subdomain_locations;} //! return the domain resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->domain_resolutions;} //! return the sizes of the cell const Rcoord & get_domain_lengths() const {return this->domain_lengths;} //! return the real space coordinates of the center of pixel Rcoord get_pixel_coordinate(Ccoord & pixel) const; //! return the physical dimension of a pixel Rcoord get_pixel_lengths() const; //! check if the pixel is inside of the cell bool is_inside(Rcoord point); //! check if the point is inside of the cell bool is_inside(Ccoord pixel); //calculate the volume of each pixel: Real get_pixel_volume(); /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const override final { return this->projection->get_formulation();} /** * get a reference to the projection object. should only be * required for debugging */ Eigen::Map get_projection() { return this->projection->get_operator();} //! returns the spatial size constexpr static Dim_t get_sdim() {return DimS;}; //! return a sparse matrix adaptor to the cell Adaptor get_adaptor() override; //! returns the number of degrees of freedom in the cell Dim_t get_nb_dof() const override {return this->size()*ipow(DimS, 2);}; //! return the communicator object virtual const Communicator & get_communicator() const override { return this->projection->get_communicator(); } protected: //! make sure that every pixel is assigned to one and only one material virtual void check_material_coverage(); const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions const Ccoord & domain_resolutions; //!< the cell's domain resolutions CcoordOps::Pixels pixels; //!< helper to iterate over the pixels const Rcoord & domain_lengths; //!< the cell's lengths Collection_ptr fields; //!< handle for the global fields of the cell StrainField_t & F; //!< ref to strain field StressField_t & P; //!< ref to stress field //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; //! container of the materials present in the cell std::vector materials{}; Projection_ptr projection; //!< handle for the projection operator bool initialised{false}; //!< to handle double initialisation right SplittedCell is_cell_splitted ; private: }; /** * lightweight resource handle wrapping a `muSpectre::Cell` or * a subclass thereof into `Eigen::EigenBase`, so it can be * interpreted as a sparse matrix by Eigen solvers */ template class CellAdaptor: public Eigen::EigenBase> { public: using Scalar = double; //!< sparse matrix traits using RealScalar = double; //!< sparse matrix traits using StorageIndex = int; //!< sparse matrix traits enum { ColsAtCompileTime = Eigen::Dynamic, MaxColsAtCompileTime = Eigen::Dynamic, RowsAtCompileTime = Eigen::Dynamic, MaxRowsAtCompileTime = Eigen::Dynamic, IsRowMajor = false }; //! constructor CellAdaptor(Cell & cell):cell{cell}{} //!returns the number of logical rows Eigen::Index rows() const {return this->cell.get_nb_dof();} //!returns the number of logical columns Eigen::Index cols() const {return this->rows();} //! implementation of the evaluation template Eigen::Product operator*(const Eigen::MatrixBase& x) const { return Eigen::Product (*this, x.derived()); } Cell & cell; //!< ref to the cell }; } // muSpectre namespace Eigen { namespace internal { //! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` through a //! specialization of `Eigen::internal::generic_product_impl`: template struct generic_product_impl // GEMV stands for matrix-vector : generic_product_impl_base > { //! undocumented typedef typename Product::Scalar Scalar; //! undocumented template static void scaleAndAddTo(Dest& dst, const CellAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/) { // This method should implement "dst += alpha * lhs * rhs" inplace, // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, dst.noalias() += const_cast(lhs).cell. evaluate_projected_directional_stiffness(rhs); } }; } } #endif /* CELL_BASE_H */ diff --git a/src/cell/cell_factory.hh b/src/cell/cell_factory.hh index 182193c..e673b5b 100644 --- a/src/cell/cell_factory.hh +++ b/src/cell/cell_factory.hh @@ -1,177 +1,190 @@ /** * @file cell_factory.hh * * @author Till Junge * * @date 15 Dec 2017 * * @brief Cell factories to help create cells with ease * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef CELL_FACTORY_H #define CELL_FACTORY_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/projection_small_strain.hh" #include "fft/fftw_engine.hh" #ifdef WITH_MPI #include "common/communicator.hh" #include "fft/fftwmpi_engine.hh" #endif #include namespace muSpectre { /** * Create a unique ptr to a Projection operator (with appropriate * FFT_engine) to be used in a cell constructor */ template > inline std::unique_ptr> cell_input(Ccoord_t resolutions, Rcoord_t lengths, Formulation form ) { auto fft_ptr{ std::make_unique(resolutions, dof_for_formulation(form, DimS))}; switch (form) { case Formulation::finite_strain: { using Projection = ProjectionFiniteStrainFast; return std::make_unique(std::move(fft_ptr), lengths); break; } case Formulation::small_strain: { using Projection = ProjectionSmallStrain; return std::make_unique(std::move(fft_ptr), lengths); break; } default: { throw std::runtime_error("unknow formulation"); break; } } } /** * convenience function to create a cell (avoids having to build * and move the chain of unique_ptrs */ template , typename FFTEngine=FFTWEngine> inline Cell make_cell(Ccoord_t resolutions, Rcoord_t lengths, Formulation form) { auto && input = cell_input(resolutions, lengths, form); auto cell{Cell{std::move(input)}}; return cell; } template , typename FFTEngine=FFTWEngine> inline Cell make_cell_split(Ccoord_t resolutions, Rcoord_t lengths, Formulation form) { auto && input = cell_input(resolutions, lengths, form); auto cell{Cell{std::move(input)}}; return cell; } + // this function returns a unique pointer to the CellBase class of the cell + // all members of cell and its descending cell class such as CellSplit are available + template , + typename FFTEngine=FFTWEngine> + std::unique_ptr> make_cell_ptr (Ccoord_t resolutions, + Rcoord_t lengths, + Formulation form) { + auto && input = cell_input(resolutions, + lengths, + form); + return std::make_unique(std::move(input)); + } #ifdef WITH_MPI /** * Create a unique ptr to a parallel Projection operator (with appropriate * FFT_engine) to be used in a cell constructor */ template > inline std::unique_ptr> parallel_cell_input(Ccoord_t resolutions, Rcoord_t lengths, Formulation form, const Communicator & comm) { auto fft_ptr{std::make_unique(resolutions, dof_for_formulation(form, DimM), comm)}; switch (form) { case Formulation::finite_strain: { using Projection = ProjectionFiniteStrainFast; return std::make_unique(std::move(fft_ptr), lengths); break; } case Formulation::small_strain: { using Projection = ProjectionSmallStrain; return std::make_unique(std::move(fft_ptr), lengths); break; } default: { throw std::runtime_error("unknown formulation"); break; } } } /** * convenience function to create a parallel cell (avoids having to build * and move the chain of unique_ptrs */ template , typename FFTEngine=FFTWMPIEngine> inline Cell make_parallel_cell(Ccoord_t resolutions, Rcoord_t lengths, Formulation form, const Communicator & comm) { auto && input = parallel_cell_input(resolutions, lengths, form, comm); auto cell{Cell{std::move(input)}}; return cell; } #endif /* WITH_MPI */ } // muSpectre #endif /* CELL_FACTORY_H */ diff --git a/src/cell/cell_split.cc b/src/cell/cell_split.cc index ce2a3d1..94b545b 100644 --- a/src/cell/cell_split.cc +++ b/src/cell/cell_split.cc @@ -1,282 +1,273 @@ /** * @file cell_split.cc * * @author Ali Falsafi * * @date 19 Apr 2018 * * @brief Implementation for cell base class * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include "cell/cell_base.hh" #include "cell/cell_split.hh" +#include "common/intersection_octree.hh" #include #include #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template CellSplit::CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted ) :Parent(std::move(projection), is_cell_splitted) {} /* ---------------------------------------------------------------------- */ template std::vector CellSplit::get_assigned_ratios(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } return pixel_assigned_ratios; } /* ---------------------------------------------------------------------- */ template - std::vector CellSplit::get_unassigned_ratios_incompleted_pixels(){ + std::vector CellSplit::get_unassigned_ratios_incomplete_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); - std::vector pixel_assigned_ratios_incompleted_pixels; + std::vector pixel_assigned_ratios_incomplete_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && ratio:pixel_assigned_ratios){ if (ratio<1){ - pixel_assigned_ratios_incompleted_pixels.push_back(1.0-ratio); + pixel_assigned_ratios_incomplete_pixels.push_back(1.0-ratio); } } - return pixel_assigned_ratios_incompleted_pixels; + return pixel_assigned_ratios_incomplete_pixels; } /* ---------------------------------------------------------------------- */ template - auto CellSplit::make_incompleted_pixels() -> IncompletedPixels{ - return IncompletedPixels(*this); + auto CellSplit::make_incomplete_pixels() -> IncompletePixels{ + return IncompletePixels(*this); } /* ---------------------------------------------------------------------- */ template - std::vector CellSplit::get_index_incompleted_pixels(){ + std::vector CellSplit::get_index_incomplete_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); std::vector index_unassigned_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && tup: akantu::enumerate(this->pixels)){ auto && i{std::get<0> (tup)}; if (pixel_assigned_ratios[i] < 1){ index_unassigned_pixels.push_back(i); } } return index_unassigned_pixels; } /* ---------------------------------------------------------------------- */ template std::vector> CellSplit::get_unassigned_pixels(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratios(nb_pixels, 0.0); std::vector> unassigned_pixels; for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratios.at(index) += mat->get_assigned_ratio(pixel); } } for (auto && tup: akantu::enumerate(this->pixels)){ auto && index{std::get<0> (tup)}; auto && pix{std::get<1> (tup)}; if (pixel_assigned_ratios[index] < 1){ unassigned_pixels.push_back(pix); } } return unassigned_pixels; } /* ---------------------------------------------------------------------- */ - template - template - void CellSplit::complete_material_assignment(MaterialType& material){ - std::vector pixel_assigned_ratio(this->get_assigned_ratios()); - for (auto && tup: akantu::enumerate(*this)) { - auto && pixel = std::get<1>(tup); - auto iterator = std::get<0>(tup); - if (pixel_assigned_ratio[iterator] < 1.0){ - material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]); - } - } - } - template void CellSplit::check_material_coverage(){ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector pixel_assigned_ratio(nb_pixels, 0.0); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); pixel_assigned_ratio.at(index) += mat->get_assigned_ratio(pixel); } } std::vector> over_assigned_pixels; std::vector> under_assigned_pixels; for (size_t i = 0; i < nb_pixels; ++i) { if (pixel_assigned_ratio.at(i) > 1.0) { over_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, this->subdomain_locations, i)); }else if (pixel_assigned_ratio.at(i) < 1.0) { under_assigned_pixels.push_back(CcoordOps::get_ccoord(this->subdomain_resolutions, this->subdomain_locations, i)); } } if (over_assigned_pixels.size() != 0) { std::stringstream err {}; err << "Execesive material is assigned to the following pixels: "; for (auto & pixel: over_assigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } if (under_assigned_pixels.size() != 0) { std::stringstream err {}; err << "Insufficient material is assigned to the following pixels: "; for (auto & pixel: under_assigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } } /* ---------------------------------------------------------------------- */ //this piece of code handles the evaluation of stress an dtangent matrix //in case the cells have materials in which pixels are partially composed of //diffferent materials. template typename CellSplit::FullResponse_t CellSplit::evaluate_stress_tangent(StrainField_t & F){ return evaluate_split_stress_tangent(F); } /* ---------------------------------------------------------------------- */ template typename CellSplit::FullResponse_t CellSplit::evaluate_split_stress_tangent(StrainField_t & F){ if (this->initialised == false) { this->initialise(); } //! High level compatibility checks if (F.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); // Here we should first set P and K matrixes equal to zeros first because we want to add up contribution //of the partial influence of different materials assigend to each pixel. Therefore, this values should be // initiialised as zero filled tensors this->set_p_k_zero(); for (auto & mat: this->materials) { mat->compute_stresses_tangent(F, this->P, this->K.value(), this->get_formulation(), this->is_cell_splitted); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template void CellSplit::set_p_k_zero(){ this->P.set_zero(); this->K.value().get().set_zero(); } /* ---------------------------------------------------------------------- */ template - CellSplit::IncompletedPixels::IncompletedPixels(CellSplit& cell) - :cell(cell),incomplete_assigned_ratios(cell.get_unassigned_ratios_incompleted_pixels()), - index_incompleted_pixels(cell.get_index_incompleted_pixels()) + CellSplit::IncompletePixels::IncompletePixels(CellSplit& cell) + :cell(cell),incomplete_assigned_ratios(cell.get_unassigned_ratios_incomplete_pixels()), + index_incomplete_pixels(cell.get_index_incomplete_pixels()) {} /* ---------------------------------------------------------------------- */ template - CellSplit::IncompletedPixels::iterator::iterator - (const IncompletedPixels & pixels, bool begin) - : incompleted_pixels(pixels), - index{begin? 0:this->incompleted_pixels.index_incompleted_pixels.size()} + CellSplit::IncompletePixels::iterator::iterator + (const IncompletePixels & pixels, bool begin) + : incomplete_pixels(pixels), + index{begin? 0:this->incomplete_pixels.index_incomplete_pixels.size()} {} /* ---------------------------------------------------------------------- */ template - auto CellSplit::IncompletedPixels::iterator::operator++()->iterator&{ + auto CellSplit::IncompletePixels::iterator::operator++()->iterator&{ this->index++; return *this; } /* ---------------------------------------------------------------------- */ template - auto CellSplit::IncompletedPixels::iterator::operator!=(const iterator & other)->bool{ - // return (this->incompleted_pixels.index_incompleted_pixels[index] != - // other.incompleted_pixels.index_incompleted_pixels[index]); + auto CellSplit::IncompletePixels::iterator::operator!=(const iterator & other)->bool{ + // return (this->incomplete_pixels.index_incomplete_pixels[index] != + // other.incomplete_pixels.index_incomplete_pixels[index]); return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template - auto CellSplit::IncompletedPixels::iterator::operator*()->value_type const{ - auto ccoord = CcoordOps::get_ccoord(this->incompleted_pixels.cell.get_domain_resolutions(), - this->incompleted_pixels.cell.get_subdomain_locations(), - this->incompleted_pixels.index_incompleted_pixels[index]); - auto ratio = this->incompleted_pixels.incomplete_assigned_ratios[index]; + auto CellSplit::IncompletePixels::iterator::operator*()->value_type const{ + auto ccoord = CcoordOps::get_ccoord(this->incomplete_pixels.cell.get_domain_resolutions(), + this->incomplete_pixels.cell.get_subdomain_locations(), + this->incomplete_pixels.index_incomplete_pixels[index]); + auto ratio = this->incomplete_pixels.incomplete_assigned_ratios[index]; return std::make_tuple(ccoord, ratio); } + + /* ---------------------------------------------------------------------- */ + /* ---------------------------------------------------------------------- */ template class CellSplit; template class CellSplit; /* ---------------------------------------------------------------------- */ } //muspectre diff --git a/src/cell/cell_split.hh b/src/cell/cell_split.hh index efddc71..20820b6 100644 --- a/src/cell/cell_split.hh +++ b/src/cell/cell_split.hh @@ -1,163 +1,193 @@ /** * @file cell_split.hh * * @author Ali Falsafi * * @date 19 Apr 2018 * * @brief Base class representing a unit cell able to handle * split material assignments * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef CELL_SPLIT_H #define CELL_SPLIT_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include "cell/cell_base.hh" +#include "common/intersection_octree.hh" #include #include #include #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) /* This class handles the cells that has splitly assigned material to their pixels */ template class CellSplit: public CellBase { friend class CellBase; public: using Parent = CellBase; //!< base class //! global field collection using FieldCollection_t = GlobalFieldCollection; using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; //! Default constructor CellSplit() = delete; //! constructor using sizes and resolution CellSplit(Projection_ptr projection, SplittedCell is_cell_splitted = SplittedCell::yes); //! Copy constructor CellSplit(const CellSplit &other) = delete; //! Move constructor CellSplit(CellSplit &&other) = default; //! Destructor virtual ~CellSplit() = default; //! Copy assignment operator CellSplit& operator=(const CellSplit &other) = delete; //! Move assignment operator CellSplit& operator=(CellSplit &&other) = default; //completes the assignmnet of material with a specific material so all the under-assigned pixels would be assigned to a material - template - void complete_material_assignment(MaterialType& material); + void complete_material_assignment(MaterialBase& material); //returns the assigend ratios to each pixel std::vector get_assigned_ratios(); - std::vector get_unassigned_ratios_incompleted_pixels(); - std::vector get_index_incompleted_pixels(); + // + //template + void make_automatic_precipitate(std::vector> preciptiate_vertices, MaterialBase& material); + // + std::vector get_unassigned_ratios_incomplete_pixels(); + std::vector get_index_incomplete_pixels(); std::vector> get_unassigned_pixels(); - class IncompletedPixels{ + class IncompletePixels{ public: //! constructor - IncompletedPixels(CellSplit& cell); + IncompletePixels(CellSplit& cell); //! copy constructor - IncompletedPixels(const IncompletedPixels & other) = default; + IncompletePixels(const IncompletePixels & other) = default; //! move constructor - IncompletedPixels(IncompletedPixels & other) = default; + IncompletePixels(IncompletePixels & other) = default; //Deconstructor - virtual ~IncompletedPixels() = default; + virtual ~IncompletePixels() = default; //! iterator type over all incompletetedly assigned pixel's class iterator{ public: using value_type = std::tuple, Real>; //!< stl conformance using const_value_type = const value_type; //!< stl conformance using pointer = value_type*; //!< stl conformance using difference_type = std::ptrdiff_t; //!< stl conformance using iterator_category = std::forward_iterator_tag;//!value_type const; //! pre-increment auto operator++()->iterator&; //! inequality auto operator!=(const iterator & other)->bool; //! equality inline bool operator==(const iterator & other) const; private: - const IncompletedPixels& incompleted_pixels; + const IncompletePixels& incomplete_pixels; size_t index; }; inline iterator begin() const {return iterator(*this);} //! stl conformance inline iterator end() const {return iterator(*this, false);} //! stl conformance inline size_t size() const {return CcoordOps::get_size(this->cell.get_domain_resolutions());} private: CellSplit& cell; std::vector incomplete_assigned_ratios; - std::vector index_incompleted_pixels; + std::vector index_incomplete_pixels; }; - auto make_incompleted_pixels() -> IncompletedPixels; + auto make_incomplete_pixels() -> IncompletePixels; protected: void check_material_coverage() override final; void set_p_k_zero(); //full resppnse is consisted of the stresses and tangent matrix FullResponse_t evaluate_stress_tangent(StrainField_t & F) override final; FullResponse_t evaluate_split_stress_tangent(StrainField_t & F); private: }; -} + /* ---------------------------------------------------------------------- */ + template + void CellSplit::make_automatic_precipitate(std::vector> precipitate_vertices, MaterialBase& material){ + RootNode precipitate (*this, precipitate_vertices); + auto precipitate_intersects= precipitate.get_intersected_pixels(); + auto precipitate_intersection_ratios = precipitate.get_intersection_ratios(); + for (auto tup:akantu::zip(precipitate_intersects, precipitate_intersection_ratios)){ + auto pix = std::get<0> (tup); + auto ratio = std::get<1> (tup); + material.add_pixel_split(pix,ratio); + } + } + /* ---------------------------------------------------------------------- */ + template + void CellSplit::complete_material_assignment(MaterialBase& material){ + std::vector pixel_assigned_ratio(this->get_assigned_ratios()); + for (auto && tup: akantu::enumerate(*this)) { + auto && pixel = std::get<1>(tup); + auto iterator = std::get<0>(tup); + if (pixel_assigned_ratio[iterator] < 1.0){ + material.add_pixel_split(pixel, 1.0-pixel_assigned_ratio[iterator]); + } + } + } + +} //muspectre + #endif /* CELL_SPLIT_H */ diff --git a/src/common/field_collection_local.hh b/src/common/field_collection_local.hh index 1039eee..6369d34 100644 --- a/src/common/field_collection_local.hh +++ b/src/common/field_collection_local.hh @@ -1,184 +1,204 @@ /** * @file field_collection_local.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for local fields * * 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_COLLECTION_LOCAL_H #define FIELD_COLLECTION_LOCAL_H #include "common/field_collection_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ /** `LocalFieldCollection` derives from `FieldCollectionBase` and stores * local fields, i.e. fields that are only defined for a subset of all pixels * in the computational domain. The coordinates of these active pixels are * explicitly stored by this field collection. * `LocalFieldCollection::add_pixel` allows to add individual pixels to the * field collection. */ template class LocalFieldCollection: public FieldCollectionBase> { public: //! for compile time check constexpr static bool Global{false}; //! base class using Parent = FieldCollectionBase>; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Field_p = typename Parent::Field_p; //!< field pointer using ccoords_container = std::vector; //!< list of pixels //! iterator over managed pixels using iterator = typename ccoords_container::iterator; //! Default constructor LocalFieldCollection(); //! Copy constructor LocalFieldCollection(const LocalFieldCollection &other) = delete; //! Move constructor LocalFieldCollection(LocalFieldCollection &&other) = delete; //! Destructor virtual ~LocalFieldCollection() = default; //! Copy assignment operator LocalFieldCollection& operator=(const LocalFieldCollection &other) = delete; //! Move assignment operator LocalFieldCollection& operator=(LocalFieldCollection &&other) = delete; //! add a pixel/voxel to the field collection inline void add_pixel(const Ccoord & local_ccoord); + //! add a split pixel/voxel to the field collection + inline void add_pixel_split(const Ccoord & local_ccoord, Real ratio); + /** allocate memory, etc. at this point, the field collection knows how many entries it should have from the size of the coords containes (which grows by one every time add_pixel is called. The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' any field of a different size is wrong TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(); //! returns the linear index corresponding to cell coordinates template inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; //! iterator to first pixel inline iterator begin() {return this->ccoords.begin();} //! iterator past last pixel inline iterator end() {return this->ccoords.end();} protected: //! container of pixel coords for non-global collections ccoords_container ccoords{}; //! container of indices for non-global collections (slow!) std::map indices{}; + // + std::vector assigned_ratio{}; private: }; /* ---------------------------------------------------------------------- */ template LocalFieldCollection::LocalFieldCollection() :Parent(){} /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: add_pixel(const Ccoord & local_ccoord) { if (this->is_initialised) { throw FieldCollectionError ("once a field collection has been initialised, you can't add new " "pixels."); } this->indices[local_ccoord] = this->ccoords.size(); this->ccoords.push_back(local_ccoord); this->size_++; } + /* ---------------------------------------------------------------------- */ + template + void LocalFieldCollection:: + add_pixel_split(const Ccoord & local_ccoord, Real ratio) { + if (this->is_initialised) { + throw FieldCollectionError + ("once a field collection has been initialised, you can't add new " + "pixels."); + } + this->indices[local_ccoord] = this->ccoords.size(); + this->ccoords.push_back(local_ccoord); + this->assigned_ratio.push_back(ratio); + this->size_++; + + } /* ---------------------------------------------------------------------- */ template void LocalFieldCollection:: initialise() { if (this->is_initialised) { throw std::runtime_error("double initialisation"); } std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template template size_t LocalFieldCollection::get_index(CcoordRef && ccoord) const { static_assert(std::is_same< Ccoord, std::remove_const_t< std::remove_reference_t>>::value, "can only be called with values or references of Ccoord"); return this->indices.at(std::forward(ccoord)); } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename LocalFieldCollection::Ccoord LocalFieldCollection::get_ccoord(size_t index) const { return this->ccoords[std::move(index)]; } } // muSpectre #endif /* FIELD_COLLECTION_LOCAL_H */ diff --git a/src/common/intersection_octree.cc b/src/common/intersection_octree.cc index bbf4fb2..0a4d256 100644 --- a/src/common/intersection_octree.cc +++ b/src/common/intersection_octree.cc @@ -1,193 +1,192 @@ /** * @file intersection_octree.cc * * @author Ali Falsafi * * @date May 2018 * * @brief common operations on pixel addressing * * Copyright © 2018 Ali Falsafi * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/ccoord_operations.hh" #include "common/common.hh" #include "cell/cell_base.hh" #include "materials/material_base.hh" #include "common/intersection_octree.hh" #include "common/intersection_volume_calculator.hh" #include #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template Node::Node(const Rcoord &new_origin, const Ccoord &new_lengths, int depth, RootNode_t& root, bool is_root ) :root_node(root), origin(new_origin), Clengths(new_lengths), depth(depth), is_pixel((depth == root.max_depth)), children_no (((is_pixel)? 0 : pow(2, DimS))) { for (int i = 0 ; i < DimS ; i++){ this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; } if (not is_root){ this->check_node(); } } - + /* ---------------------------------------------------------------------- */ template void Node::check_node() { Real total_volume = 1.0, intersected_volume = 0.0, intersection_ratio = 0.0; constexpr Real machine_prescision = 1e-12; for (int i = 0; i < DimS; i++){ total_volume *= this->Rlengths[i]; } // this volume should be calculated by CGAL as the intersection volume of the precipitate and the Node intersected_volume = total_volume ; intersected_volume = intersection_volume_calculator(this->root_node.precipitate_vertices, this->origin, this->Rlengths); intersection_ratio = intersected_volume / total_volume; if (intersection_ratio > (1.0 - machine_prescision)) { // all the pixels in this node should be assigned to the precipitate material Real pix_num = pow(2, (this->root_node.max_depth - this->depth)); Ccoord origin_point, pixels_number; for (int i = 0; i < DimS; i++ ) { origin_point[i] = std::round(this->origin[i] / this->root_node.pixel_lengths[i]); pixels_number[i] = pix_num; } CcoordOps::Pixels pixels(pixels_number, origin_point); for (auto && pix:pixels){ this->root_node.intersected_pixels.push_back(pix); this->root_node.intersection_ratios.push_back(1.0); } }else if (intersection_ratio > machine_prescision) { this->split_node(intersection_ratio); } } - - + /* ---------------------------------------------------------------------- */ template void Node::split_node(Real intersection_ratio){ if (this->depth == this->root_node.max_depth) { Ccoord pixel; for (int i = 0 ; i < DimS ; i++){ pixel[i] = std::round(this->origin[i] / this->root_node.pixel_lengths[i]); } this->root_node.intersected_pixels.push_back(pixel); this->root_node.intersection_ratios.push_back(intersection_ratio); } else { this->divide_node(); } } - + /* ---------------------------------------------------------------------- */ template void Node::divide_node(){ Rcoord new_origin; Ccoord new_length; this->children.reserve(children_no); CcoordOps::Pixels sub_nodes (CcoordOps::get_cube(2)); for (auto && sub_node: sub_nodes ) { for (int i = 0 ; i < DimS ; i++){ new_length[i] = std::round(this->Clengths[i] *0.5); new_origin[i] = this->origin[i] + sub_node[i] * Rlengths[i] * 0.5; } this->children.emplace_back(new_origin, new_length, this->depth+1, this->root_node, false); } } - + /* ---------------------------------------------------------------------- */ template RootNode::RootNode(CellBase& cell, std::vector vert_precipitate) /*Calling parent constructing method simply by argumnets of (coordinates of origin, 2^depth_max in each direction, 0 ,*this)*/ :Node(Rcoord{}, CcoordOps::get_cube(pow(2, std::ceil(std::log2(cell.get_domain_resolutions().at(std::distance(std::max_element (cell.get_domain_resolutions().begin(), cell.get_domain_resolutions().end()) ,cell.get_domain_resolutions().begin()) ) ) ) ) ), 0, *this, true), cell(cell), cell_length (cell.get_domain_lengths()), pixel_lengths (cell.get_pixel_lengths()), cell_resolution (cell.get_domain_resolutions()), max_resolution (this->cell_resolution.at(std::distance(std::max_element(this->cell_resolution.begin(), this->cell_resolution.end()) ,this->cell_resolution.begin()))), max_depth (std::ceil(std::log2(this->max_resolution))), precipitate_vertices (vert_precipitate) { for (int i = 0 ; i < DimS ; i++){ this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; } this->check_node(); } - + /* ---------------------------------------------------------------------- */ template std::vector> RootNode::get_intersected_pixels(){ return this->intersected_pixels; } - + /* ---------------------------------------------------------------------- */ template std::vector RootNode::get_intersection_ratios(){ return this->intersection_ratios; } - + /* ---------------------------------------------------------------------- */ template void RootNode:: check_root_node() { for (int i = 0 ; i < DimS ; i++){ this->Rlengths[i] = this->Clengths[i] * this->root_node.pixel_lengths[i]; } this->check_node(); } - + /* ---------------------------------------------------------------------- */ template class RootNode; template class RootNode; template class Node; template class Node; } //muSpctre diff --git a/src/common/intersection_octree.hh b/src/common/intersection_octree.hh index d9a39fe..b9fac69 100644 --- a/src/common/intersection_octree.hh +++ b/src/common/intersection_octree.hh @@ -1,145 +1,145 @@ /** * @file intersection_octree.hh * * @author Ali Falsafi * * @date May 2018 * * @brief common operations on pixel addressing * * Copyright © 2018 Ali Falsafi * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef INTERSECTION_OCTREE_H #define INTERSECTION_OCTREE_H #include #include "common/ccoord_operations.hh" #include "common/common.hh" #include "cell/cell_base.hh" #include "materials/material_base.hh" #include "common/intersection_volume_calculator.hh" namespace muSpectre { template class RootNode; template class Node { public: using Rcoord = Rcoord_t; //!< physical coordinates type using Ccoord = Ccoord_t; //!< cell coordinates type using RootNode_t = RootNode; //! Default constructor Node() = delete; // Construct by origin, lenghts, and depth Node(const Rcoord &new_origin, const Ccoord &new_lenghts, int depth, RootNode_t& root, bool is_root ); // Construct by cell - Node(const CellBase & cell, RootNode_t& root); + // Node(const CellBase& cell, RootNode_t& root); //! Copy constructor Node(const Node &other) = delete; //! Move constructor Node(Node &&other) = default; //! Destructor ~Node() = default; //This function checks the status of the node and orders its devision into //smaller nodes or asssign material to it virtual void check_node(); //This function gives the ratio of the node which happens to be inside the precipitate //and assign materila to it if node is a pixel or divide it furhter if it is not void split_node(Real ratio); //this function constructs children of a node void divide_node(); protected: //////// RootNode_t& root_node; Rcoord origin, Rlengths{}; Ccoord Clengths{}; int depth; bool is_pixel; int children_no; std::vector children{}; }; template class RootNode: public Node { friend class Node; public: using Rcoord = Rcoord_t; //!< physical coordinates type using Ccoord = Ccoord_t; //!< cell coordinates type using Parent = Node; //!< base class //!Default Constructor RootNode() = delete ; //! Constructing a root node for a cell and a preticipate inside that cell RootNode(CellBase& cell, std::vector vert_precipitate); //! Copy constructor RootNode(const RootNode &other) = delete; //! Move constructor RootNode(RootNode &&other) = default; //! Destructor ~RootNode() = default; //returns the pixels which have intersection raio with the preipitate std::vector get_intersected_pixels (); //return the intersection ratios of corresponding to the pixels returned by get_intersected_pixels() std::vector get_intersection_ratios (); //checking rootnode: void check_root_node(); protected: CellBase& cell; Rcoord cell_length, pixel_lengths; Ccoord cell_resolution; int max_resolution ; int max_depth; std::vector precipitate_vertices{}; std::vector intersected_pixels{}; std::vector intersection_ratios{}; }; } //muSpectre #endif /* INTERSECTION_OCTREE_H */ diff --git a/src/materials/material_base.cc b/src/materials/material_base.cc index 7a956c1..f26f7de 100644 --- a/src/materials/material_base.cc +++ b/src/materials/material_base.cc @@ -1,86 +1,92 @@ /** * @file material_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @brief implementation of materi * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "materials/material_base.hh" namespace muSpectre { //----------------------------------------------------------------------------// template MaterialBase::MaterialBase(std::string name) :name(name),assigned_ratio{make_field("assigned ratio", this->internal_fields)}, assigned_ratio_mapped{assigned_ratio}{ static_assert((DimM == oneD)|| (DimM == twoD)|| (DimM == threeD), "only 1, 2, or threeD supported"); } /* ---------------------------------------------------------------------- */ template const std::string & MaterialBase::get_name() const { return this->name; } /* ---------------------------------------------------------------------- */ template void MaterialBase::add_pixel(const Ccoord &ccoord) { this->internal_fields.add_pixel(ccoord); } + /* ---------------------------------------------------------------------- */ + template < Dim_t DimS, Dim_t DimM> + void MaterialBase< DimS, DimM>:: + add_pixel_split(const Ccoord &local_ccoord, Real ratio){ + this->internal_fields.add_pixel_split(local_ccoord, ratio); + } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses(const Field_t & F, Field_t & P, Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses(StrainField_t::check_ref(F), StressField_t::check_ref(P), form, is_cell_splitted); } template Real MaterialBase::get_assigned_ratio(Ccoord pixel){ return this->assigned_ratio_mapped[pixel]; } /* ---------------------------------------------------------------------- */ template void MaterialBase::compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form, SplittedCell is_cell_splitted) { this->compute_stresses_tangent(StrainField_t::check_ref(F), StressField_t::check_ref(P), TangentField_t::check_ref(K), form, is_cell_splitted); } template class MaterialBase<2, 2>; template class MaterialBase<2, 3>; template class MaterialBase<3, 3>; } // muSpectre diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 64be10c..86bfb5f 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,165 +1,167 @@ /** * @file material_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials (constitutive models) * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIAL_BASE_H #define MATERIAL_BASE_H #include "common/common.hh" #include "common/field.hh" #include "common/field_collection.hh" #include namespace muSpectre { //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) and /** * @a DimM is the material dimension (i.e., the dimension of constitutive * law; even for e.g. two-dimensional problems the constitutive law could * live in three-dimensional space for e.g. plane strain or stress problems) */ template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for cell-wide fields, like stress, strain, etc using GFieldCollection_t = GlobalFieldCollection; //! field collection for internal variables, such as eigen-strains, //! plastic strains, damage variables, etc, but also for managing which //! pixels the material is responsible for using MFieldCollection_t = LocalFieldCollection; using iterator = typename MFieldCollection_t::iterator; //!< pixel iterator //! polymorphic base class for fields only to be used for debugging using Field_t = internal::FieldBase; //! Full type for stress fields using StressField_t = TensorField; //! Full type for strain fields using StrainField_t = StressField_t; //! Full type for tangent stiffness fields fields using TangentField_t = TensorField; using Ccoord = Ccoord_t; //!< cell coordinates type //! Default constructor MaterialBase() = delete; //! Construct by name MaterialBase(std::string name); //! Copy constructor MaterialBase(const MaterialBase &other) = delete; //! Move constructor MaterialBase(MaterialBase &&other) = delete; //! Destructor virtual ~MaterialBase() = default; //! Copy assignment operator MaterialBase& operator=(const MaterialBase &other) = delete; //! Move assignment operator MaterialBase& operator=(MaterialBase &&other) = delete; /** * take responsibility for a pixel identified by its cell coordinates * WARNING: this won't work for materials with additional info per pixel * (as, e.g. for eigenstrain), we need to pass more parameters. Materials * of this tye need to overload add_pixel */ virtual void add_pixel(const Ccoord & ccooord); + virtual void add_pixel_split (const Ccoord &local_ccoord, + Real ratio); //! allocate memory, etc, but also: wipe history variables! virtual void initialise() = 0; /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, the virtual * base implementation does nothing, but materials with history * variables need to implement this */ virtual void save_history_variables() {}; //! return the material's name const std::string & get_name() const; //! spatial dimension for static inheritance constexpr static Dim_t sdim() {return DimS;} //! material dimension for static inheritance constexpr static Dim_t mdim() {return DimM;} //! computes stress virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses, mostly for debugging and * testing. Has runtime-cost associated with compatibility-checking and * conversion of the Field_t arguments that can be avoided by using the * version with strongly typed field references */ void compute_stresses(const Field_t & F, Field_t & P, Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplittedCell is_cell_splitted) = 0; /** * Convenience function to compute stresses and tangent moduli, mostly for * debugging and testing. Has runtime-cost associated with * compatibility-checking and conversion of the Field_t arguments that can * be avoided by using the version with strongly typed field references */ void compute_stresses_tangent(const Field_t & F, Field_t & P, Field_t & K, Formulation form, SplittedCell is_cell_splitted = SplittedCell::no); // this function return the ratio of which the input pixel is consisted of this material Real get_assigned_ratio(Ccoord pixel); //! iterator to first pixel handled by this material inline iterator begin() {return this->internal_fields.begin();} //! iterator past the last pixel handled by this material inline iterator end() {return this->internal_fields.end();} //! number of pixels assigned to this material inline size_t size() const {return this->internal_fields.size();} protected: const std::string name; //!< material's name (for output and debugging) MFieldCollection_t internal_fields{};//!< storage for internal variables using ScalarField_t = ScalarField, Real> ; ScalarField_t & assigned_ratio; //!< field holding the assigning ratio of the material using ScalarFieldMap_t = ScalarFieldMap, Real, true>; ScalarFieldMap_t assigned_ratio_mapped; private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 5ca234c..932f2b8 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,903 +1,911 @@ /** * @file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef MATERIAL_MUSPECTRE_BASE_H #define MATERIAL_MUSPECTRE_BASE_H #include "common/common.hh" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.hh" #include #include #include #include namespace muSpectre { // Forward declaration for factory function template class CellBase; /** * material traits are used by `muSpectre::MaterialMuSpectre` to * break the circular dependence created by the curiously recurring * template parameter. These traits must define * - these `muSpectre::FieldMap`s: * - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a * constant second-order `muSpectre::TensorField` * - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a * writable secord-order `muSpectre::TensorField` * - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a * writable fourth-order `muSpectre::TensorField` * - `strain_measure`: the expected strain type (will be replaced by the * small-strain tensor ε * `muspectre::StrainMeasure::Infinitesimal` in small * strain computations) * - `stress_measure`: the measure of the returned stress. Is used by * `muspectre::MaterialMuSpectre` to transform it into * Cauchy stress (`muspectre::StressMeasure::Cauchy`) in * small-strain computations and into first * Piola-Kirchhoff stress `muspectre::StressMeasure::PK1` * in finite-strain computations * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing * internal variables */ template struct MaterialMuSpectre_traits { }; template class MaterialMuSpectre; /** * Base class for most convenient implementation of materials */ template class MaterialMuSpectre: public MaterialBase { public: /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; //!< base class //! global field collection using GFieldCollection_t = typename Parent::GFieldCollection_t; //! expected type for stress fields using StressField_t = typename Parent::StressField_t; //! expected type for strain fields using StrainField_t = typename Parent::StrainField_t; //! expected type for tangent stiffness fields using TangentField_t = typename Parent::TangentField_t; using Ccoord = Ccoord_t; //!< cell coordinates type //! traits for the CRTP subclass using traits = MaterialMuSpectre_traits; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre &other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre &&other) = delete; //! Destructor virtual ~MaterialMuSpectre() = default; //! Factory template static Material & make(CellBase & cell, ConstructorArgs &&... args); //! Copy assignment operator MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete; //! Move assignment operator MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete; //! allocate memory, etc virtual void initialise() override; // this fucntion is speicalized to assign partial material to a pixel //void add_pixel_split(const Ccoord & local_ccoord, Real ratio = 1.0); template - void add_pixel_split(const Ccoord_t & pixel, - Real ratio, - InternalArgs ... args); + inline void add_pixel_split(const Ccoord_t & pixel, + Real ratio, + InternalArgs ... args); + + inline void add_pixel_split(const Ccoord_t & pixel, + Real ratio); + - template void add_split_pixels_precipitate(std::vector> intersected_precipitate, - std::vector intersection_ratios, - InternalArgs ... args); + std::vector intersection_ratios); //! computes stress using Parent::compute_stresses; using Parent::compute_stresses_tangent; virtual void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form, SplittedCell is_cell_splitted) override; //! computes stress and tangent modulus virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplittedCell is_cell_splitted)override; protected: //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P) __attribute__ ((visibility ("default"))); //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K) __attribute__ ((visibility ("default"))); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate over //! that does or does not include tangent moduli template class iterable_proxy; /** * inheriting classes with internal variables need to overload this function */ typename traits::InternalVariables& get_internals() { return static_cast(*this).get_internals();} bool is_initialised{false}; //!< to handle double initialisation right private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre:: MaterialMuSpectre(std::string name) :Parent(name){ using stress_compatible = typename traits::StressMap_t:: template is_compatible; using strain_compatible = typename traits::StrainMap_t:: template is_compatible; using tangent_compatible = typename traits::TangentMap_t:: template is_compatible; static_assert((stress_compatible::value && stress_compatible::explain()), "The material's declared stress map is not compatible " "with the stress field. More info in previously shown " "assert."); static_assert((strain_compatible::value && strain_compatible::explain()), "The material's declared strain map is not compatible " "with the strain field. More info in previously shown " "assert."); static_assert((tangent_compatible::value && tangent_compatible::explain()), "The material's declared tangent map is not compatible " "with the tangent field. More info in previously shown " "assert."); } /* ---------------------------------------------------------------------- */ template template Material & MaterialMuSpectre:: make(CellBase & cell, ConstructorArgs && ... args) { auto mat = std::make_unique(args...); auto & mat_ref = *mat; cell.add_material(std::move(mat)); return mat_ref; } /* ---------------------------------------------------------------------- */ template - template void MaterialMuSpectre:: - add_pixel_split(const Ccoord &local_ccoord, Real ratio, InternalArgs ...args){ + add_pixel_split(const Ccoord &local_ccoord, Real ratio) { auto & this_mat = static_cast(*this); - this_mat.add_pixel(local_ccoord, args...); + this_mat.add_pixel(local_ccoord); this->assigned_ratio.push_back(ratio); } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: + add_pixel_split(const Ccoord_t & pixel, + Real ratio, + InternalArgs ... args){ + auto & this_mat = static_cast(*this); + this_mat.add_pixel(pixel, args...); + this->assigned_ratio.push_back(ratio); + } + /* ---------------------------------------------------------------------- */ + template + void MaterialMuSpectre:: add_split_pixels_precipitate(std::vector> intersected_precipitate, - std::vector intersection_ratios, - InternalArgs ... args){ + std::vector intersection_ratios){ //assign precipitate materials: - //typedef typename std::vector>::iterator Ccoord_iter; - //typedef std::vector::iterator Real_iter; for (auto && tup: akantu::zip(intersected_precipitate, intersection_ratios) ){ auto pix = std::get<0> (tup); auto ratio = std::get<1> (tup); - this->add_pixel_split(pix,ratio , args...); + this->add_pixel_split(pix,ratio); } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: initialise() { if (!this->is_initialised) { this->internal_fields.initialise(); this->is_initialised = true; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses(const StrainField_t &F, StressField_t &P, Formulation form, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P); break; } } break; } case Formulation::small_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P); break; } } break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre:: compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form, SplittedCell is_cell_splitted) { switch (form) { case Formulation::finite_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P, K); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P, K); break; } } break; } case Formulation::small_strain: { switch (is_cell_splitted){ case (SplittedCell::no): { this->compute_stresses_worker(F, P, K); break; } case (SplittedCell::yes): { this->compute_stresses_worker(F, P, K); break; } } break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] (Strains_t Strains, Stresses_t Stresses, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli // for the case that cells do not contain any splitt cells auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ Stresses = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){ auto stress_tgt_contributions = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); auto && stress = std::get<0>(Stresses); auto && tangent = std::get<1>(Stresses); stress += ratio * std::get<0>(stress_tgt_contributions); tangent += ratio * std::get<1>(stress_tgt_contributions); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t Stresses, const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress // and tangent moduli // for the case that cells do not contain splitt cells auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); Stresses = MatTB::PK1_stress (std::move(grad), std::move(std::get<0>(stress_tgt)), std::move( std::get<1>(stress_tgt))); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){ auto stress_tgt = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...);}, internal_variables); auto && stress_tgt_contributions = MatTB::PK1_stress (std::move(grad), std::move(std::get<0>(stress_tgt)), std::move(std::get<1>(stress_tgt))); auto && stress = std::get<0>(Stresses); auto && tangent = std::get<1>(Stresses); stress += ratio * std::get<0>(stress_tgt_contributions); tangent += ratio * std::get<1>(stress_tgt_contributions); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; iterable_proxy fields{*this, F, P, K}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for tangent reference, check iterator" "value_type"); static_assert(std::is_same(arglist))>>::value, "Type mismatch for ratio reference, expected a real number"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre:: compute_stresses_worker(const StrainField_t & F, StressField_t & P){ /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this] (Strains_t Strains, Stresses_t Stresses,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // return value contains a tuple of rvalue_refs to both stress and tangent moduli auto & sigma = std::get<0>(Stresses); //for the case that cell does not contain any splitted pixel auto non_split_pixel = [&strain, &this_mat,ratio, &sigma, &internal_variables](){ sigma = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); }; // for the case that cells contain splitt cells auto split_pixel = [&strain, &this_mat, &ratio, &sigma, &internal_variables](){ sigma += ratio * apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; auto constitutive_law_finite_strain = [this] (Strains_t Strains, Stresses_t && Stresses,const Real & ratio, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // TODO: Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. //for the case that cell does not contain any splitted pixel auto non_split_pixel = [ &strain, &this_mat, &ratio, &Stresses, &internal_variables, &grad] (){ auto stress = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); auto & P = get<0>(Stresses); P = MatTB::PK1_stress (grad, stress); }; // for the case that cells contain splitted cells auto split_pixel = [&strain, &this_mat, ratio, &Stresses, internal_variables, &grad](){ auto stress = apply([&strain, &this_mat] (auto && ... internals) { return this_mat.evaluate_stress(std::move(strain), internals...);}, internal_variables); auto && P = get<0>(Stresses); P += ratio * MatTB::PK1_stress (grad, stress); }; switch (is_cell_splitted){ case SplittedCell::no : { non_split_pixel(); break; } case SplittedCell::yes : { split_pixel(); break; } } }; iterable_proxy fields{*this, F, P}; for (auto && arglist: fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_type in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ //auto && stress_tgt = std::get<0>(tuples); //auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same(std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same(std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same(arglist))>>::value, "Type mismatch for ratio reference, expected a real number"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K) :material{mat}, strain_field{F}, stress_tup{P,K}, internals{material.get_internals()}{}; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) :material{mat}, strain_field{F}, stress_tup{P}, internals{material.get_internals()}{}; //! Expected type for strain fields using StrainMap_t = typename traits::StrainMap_t; //! Expected type for stress fields using StressMap_t = typename traits::StressMap_t; //! Expected type for tangent stiffness fields using TangentMap_t = typename traits::TangentMap_t; //! expected type for strain values using Strain_t = typename traits::StrainMap_t::reference; //! expected type for stress values using Stress_t = typename traits::StressMap_t::reference; //! expected type for tangent stiffness values using Tangent_t = typename traits::TangentMap_t::reference; //! tuple of intervnal variables, depends on the material using InternalVariables = typename traits::InternalVariables; //! tuple containing a stress and possibly a tangent stiffness field using StressFieldTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness field map using StressMapTup = std::conditional_t <(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness value ref using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! Copy constructor iterable_proxy(const iterable_proxy &other) = default; //! Move constructor iterable_proxy(iterable_proxy &&other) = default; //! Destructor virtual ~iterable_proxy() = default; //! Copy assignment operator iterable_proxy& operator=(const iterable_proxy &other) = default; //! Move assignment operator iterable_proxy& operator=(iterable_proxy &&other) = default; /** * dereferences into a tuple containing strains, and internal * variables, as well as maps to the stress and potentially * stiffness maps where to write the response of a pixel */ class iterator { public: //! type to refer to internal variables owned by a CRTP material using InternalReferences = MatTB::ReferenceTuple_t; //! return type to be unpacked per pixel my the constitutive law using value_type = std::tuple, Stress_tTup, Real, InternalReferences>; using iterator_category = std::forward_iterator_tag; //!< stl conformance //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ iterator(const iterable_proxy & it, bool begin = true) : it{it}, strain_map{it.strain_field}, stress_map {it.stress_tup}, index{begin ? 0:it.material.internal_fields.size()}{} //! Copy constructor iterator(const iterator &other) = default; //! Move constructor iterator(iterator &&other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator& operator=(const iterator &other) = default; //! Move assignment operator iterator& operator=(iterator &&other) = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const iterable_proxy & it; //!< ref to the proxy StrainMap_t strain_map; //!< map onto the global strain field //! map onto the global stress field and possibly tangent stiffness StressMapTup stress_map; size_t index; //!< index or pixel currently referred to private: }; //! returns iterator to first pixel if this material iterator begin() {return std::move(iterator(*this));} //! returns iterator past the last pixel in this material iterator end() {return std::move(iterator(*this, false));} protected: MaterialMuSpectre & material; //!< reference to the proxied material const StrainField_t & strain_field; //!< cell's global strain field //! references to the global stress field and perhaps tangent stiffness StressFieldTup stress_tup; //! references to the internal variables InternalVariables & internals; private: }; /* ---------------------------------------------------------------------- */ template template bool MaterialMuSpectre::iterable_proxy::iterator:: operator!=(const iterator & other) const { return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy:: iterator & MaterialMuSpectre::iterable_proxy::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre:: template iterable_proxy::iterator:: value_type MaterialMuSpectre::iterable_proxy::iterator:: operator*() { const Ccoord_t pixel{ this->it.material.internal_fields.get_ccoord(this->index)}; auto && strain = std::make_tuple(this->strain_map[pixel]); auto && ratio = this->it.material.get_assigned_ratio(pixel); auto && stresses = apply([&pixel] (auto && ... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...);}, this->stress_map); auto && internal = this->it.material.get_internals(); const auto id{this->index}; auto && internals = apply([id] (auto && ... internals_) { return InternalReferences{internals_[id]...};}, internal); return std::make_tuple(std::move(strain), std::move(stresses), std::move(ratio), std::move(internals)); } } // muSpectre #endif /* MATERIAL_MUSPECTRE_BASE_H */ diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index 7a8c2e1..5c9285d 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,172 +1,172 @@ #!/usr/bin/env python3 """ file python_binding_tests.py @author Till Junge @date 09 Jan 2018 @brief Unit tests for python bindings @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 unittest import numpy as np from python_test_imports import µ from python_fft_tests import FFT_Check from python_projection_tests import * from python_material_linear_elastic3_test import MaterialLinearElastic3_Check -from python_material_linear_elastic4_test import MaterialLinearElastic4_Check +#from python_material_linear_elastic4_test import MaterialLinearElastic4_Check class CellCheck(unittest.TestCase): def test_Construction(self): """ Simple check for cell constructors """ resolution = [5,7] lengths = [5.2, 8.3] formulation = µ.Formulation.small_strain try: sys = µ.Cell(resolution, lengths, formulation) mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material", 210e9, .33) except Exception as err: print(err) raise err class MaterialLinearElastic1_2dCheck(unittest.TestCase): def setUp(self): self.resolution = [5,7] self.lengths = [5.2, 8.3] self.formulation = µ.Formulation.small_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat = µ.material.MaterialLinearElastic1_2d.make( self.sys, "material", 210e9, .33) def test_add_material(self): self.mat.add_pixel([2,1]) class SolverCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.finite_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.hard = µ.material.MaterialLinearElastic1_2d.make( self.sys, "hard", 210e9, .33) self.soft = µ.material.MaterialLinearElastic1_2d.make( self.sys, "soft", 70e9, .33) def test_solve(self): for i, pixel in enumerate(self.sys): if i < 3: self.hard.add_pixel(pixel) else: self.soft.add_pixel(pixel) self.sys.initialise() tol = 1e-6 Del0 = np.array([[0, .1], [0, 0]]) maxiter = 100 verbose = 0 solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose) r = µ.solvers.de_geus(self.sys, Del0, solver,tol, verbose) #print(r) class EigenStrainCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.small_strain self.cell1 = µ.Cell(self.resolution, self.lengths, self.formulation) self.cell2 = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat1 = µ.material.MaterialLinearElastic1_2d.make( self.cell1, "simple", 210e9, .33) self.mat2 = µ.material.MaterialLinearElastic2_2d.make( self.cell2, "eigen", 210e9, .33) def test_solve(self): verbose_test = False if verbose_test: print("start test_solve") grad = np.array([[1.1, .2], [ .3, 1.5]]) gl_strain = -0.5*(grad.T.dot(grad) - np.eye(2)) gl_strain = -0.5*(grad.T + grad - 2*np.eye(2)) grad = -gl_strain if verbose_test: print("grad =\n{}\ngl_strain =\n{}".format(grad, gl_strain)) for i, pixel in enumerate(self.cell1): self.mat1.add_pixel(pixel) self.mat2.add_pixel(pixel, gl_strain) self.cell1.initialise() self.cell2.initialise() tol = 1e-6 Del0_1 = grad Del0_2 = np.zeros_like(grad) maxiter = 2 verbose = 0 def solve(cell, grad): solver=µ.solvers.SolverCG(cell, tol, maxiter, verbose) r = µ.solvers.newton_cg(cell, grad, solver, tol, tol, verbose) return r results = [solve(cell, del0) for (cell, del0) in zip((self.cell1, self.cell2), (Del0_1, Del0_2))] P1 = results[0].stress P2 = results[1].stress error = np.linalg.norm(P1-P2)/np.linalg.norm(.5*(P1+P2)) if verbose_test: print("cell 1, no eigenstrain") print("P1:\n{}".format(P1[:,0])) print("F1:\n{}".format(results[0].grad[:,0])) print("cell 2, with eigenstrain") print("P2:\n{}".format(P2[:,0])) print("F2:\n{}".format(results[1].grad[:,0])) print("end test_solve") self.assertLess(error, tol) if __name__ == '__main__': unittest.main() diff --git a/tests/split_test_intersection_error_induced.cc b/tests/split_test_intersection_error_induced.cc index 5b80747..b353919 100644 --- a/tests/split_test_intersection_error_induced.cc +++ b/tests/split_test_intersection_error_induced.cc @@ -1,210 +1,199 @@ /** * @file test_intersection_error_induced.cc * * @author Ali Faslfi * * @date 21 Jun 2018 * * @brief Tests for split cells and octree material assignment * * Copyright © 2017 Till Ali Faslafi * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "tests.hh" #include "solver/deprecated_solvers.hh" #include "solver/deprecated_solver_cg.hh" #include "solver/deprecated_solver_cg_eigen.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain_fast.hh" #include "materials/material_linear_elastic1.hh" #include "common/iterators.hh" #include "common/ccoord_operations.hh" #include "common/common.hh" #include "cell/cell_factory.hh" #include "cell/cell_split.hh" #include "common/intersection_octree.hh" #include #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(split_cell_octree_error_test); BOOST_AUTO_TEST_CASE(area_calculation_and_eq_solution_of_35vs7_pixels) { constexpr Dim_t dim{twoD}; using Rcoord = Rcoord_t; using Ccoord = Ccoord_t; using Mat_t = MaterialLinearElastic1; const Real contrast {10}; const Real Young_soft{1.0030648180242636}, Poisson_soft{0.29930675909878679}; const Real Young_hard{contrast * Young_soft}, Poisson_hard{0.29930675909878679}; const Real machine_precision = 1e-13; constexpr Real length {11}; constexpr int high_res{55},low_res{11}; constexpr Rcoord lengths_split{length, length}; constexpr Ccoord resolutions_split_high_res{high_res, high_res}; constexpr Ccoord resolutions_split_low_res{low_res, low_res}; - auto fft_ptr_split_high_res{std::make_unique>(resolutions_split_high_res, ipow(dim, 2))}; auto fft_ptr_split_low_res{std::make_unique>(resolutions_split_low_res, ipow(dim, 2))}; - auto proj_ptr_split_high_res{std::make_unique> (std::move(fft_ptr_split_high_res), lengths_split)}; auto proj_ptr_split_low_res{std::make_unique> (std::move(fft_ptr_split_low_res), lengths_split)}; CellSplit sys_split_high_res(std::move(proj_ptr_split_high_res), SplittedCell::yes); CellSplit sys_split_low_res(std::move(proj_ptr_split_low_res), SplittedCell::yes); auto& Material_hard_split_high_res = Mat_t::make(sys_split_high_res, "hard", Young_hard, Poisson_hard); auto& Material_soft_split_high_res = Mat_t::make(sys_split_high_res, "soft", Young_soft, Poisson_soft); auto& Material_hard_split_low_res = Mat_t::make(sys_split_low_res, "hard", Young_hard, Poisson_hard); auto& Material_soft_split_low_res = Mat_t::make(sys_split_low_res, "soft", Young_soft, Poisson_soft); constexpr Rcoord center{5.5, 5.5}; constexpr Rcoord halfside1{2.3, 2.3}; constexpr Rcoord halfside2{2.4, 2.4}; std::vector precipitate_vertices ; precipitate_vertices.push_back ({center[0] - halfside1[0], center[1] - halfside1[1]}); precipitate_vertices.push_back ({center[0] + halfside2[0], center[1] - halfside1[1]}); precipitate_vertices.push_back ({center[0] - halfside1[0], center[1] + halfside2[1]}); precipitate_vertices.push_back ({center[0] + halfside2[0], center[1] + halfside2[1]}); - //assign precipitate materials: - typedef std::vector::iterator Ccoord_iter_low_res, Ccoord_iter_high_res; - typedef std::vector::iterator Real_iter_low_res, Real_iter_high_res; - RootNode precipitate_high_res (sys_split_high_res, precipitate_vertices); - //Extracting the intersected pixels and their correspondent intersection ratios: auto precipitate_high_res_intersects = precipitate_high_res.get_intersected_pixels(); auto precipitate_high_res_intersection_ratios = precipitate_high_res.get_intersection_ratios(); - - //assign precipitate materials: - for (std::pair - it(precipitate_high_res_intersects.begin(), precipitate_high_res_intersection_ratios.begin()) ; - it.first != precipitate_high_res_intersects.end(); - ++it.first , ++it.second ){ - Material_hard_split_high_res.add_pixel_split(*(it.first), *(it.second)); - //Material_soft_split_high_res.add_pixel_split(*(it.first), 1.0-*(it.second)); + for (auto tup:akantu::zip(precipitate_high_res_intersects, precipitate_high_res_intersection_ratios)){ + auto pix = std::get<0> (tup); + auto ratio = std::get<1>(tup); + Material_hard_split_high_res.add_pixel_split(pix,ratio); + Material_soft_split_high_res.add_pixel_split(pix, 1.0-ratio); } std::vector assigned_ratio_high_res = sys_split_high_res.get_assigned_ratios(); - int iterator = 0 ; for (auto && tup: akantu::enumerate(sys_split_high_res)) { auto && pixel = std::get<1>(tup); + auto iterator = std::get<0>(tup); if (assigned_ratio_high_res[iterator] < 1.0){ Material_soft_split_high_res.add_pixel_split(pixel, 1.0-assigned_ratio_high_res[iterator]); } - iterator++; } - RootNode precipitate_low_res (sys_split_low_res, precipitate_vertices); + RootNode precipitate_low_res (sys_split_low_res, precipitate_vertices); //Extracting the intersected pixels and their correspondent intersection ratios: auto precipitate_low_res_intersects = precipitate_low_res.get_intersected_pixels(); auto precipitate_low_res_intersection_ratios = precipitate_low_res.get_intersection_ratios(); - - for (std::pair it(precipitate_low_res_intersects.begin(), precipitate_low_res_intersection_ratios.begin()) ; - it.first != precipitate_low_res_intersects.end(); ++it.first , ++it.second ){ - - Material_hard_split_low_res.add_pixel_split(*(it.first), *(it.second)); - Material_soft_split_low_res.add_pixel_split(*(it.first), 1.0-*(it.second)); + for (auto tup:akantu::zip(precipitate_low_res_intersects, precipitate_low_res_intersection_ratios)){ + auto pix = std::get<0> (tup); + auto ratio = std::get<1>(tup); + Material_hard_split_low_res.add_pixel_split(pix,ratio); + Material_soft_split_low_res.add_pixel_split(pix, 1.0-ratio); } std::vector assigned_ratio_low_res = sys_split_low_res.get_assigned_ratios(); - iterator = 0 ; for (auto && tup: akantu::enumerate(sys_split_low_res)) { auto && pixel = std::get<1>(tup); + auto iterator = std::get<0>(tup); if (assigned_ratio_low_res[iterator] != 1.0){ Material_soft_split_low_res.add_pixel_split(pixel, 1.0-assigned_ratio_low_res[iterator]); } - iterator++; } + Real area_preticipitate_low_res = 0.0 ,area_preticipitate_high_res = 0.0 ; for (auto&& precepetiate_area:precipitate_low_res_intersection_ratios){ area_preticipitate_low_res += precepetiate_area*sys_split_low_res.get_pixel_volume(); } for (auto&& precepetiate_area:precipitate_high_res_intersection_ratios){ area_preticipitate_high_res += precepetiate_area*sys_split_high_res.get_pixel_volume(); } + BOOST_CHECK_LE(abs(area_preticipitate_high_res - area_preticipitate_low_res), area_preticipitate_low_res * machine_precision); ///////////////////////////////////////////////////////////////////////////////////////////////// //Finding equilbrium strain fileds with Sppectral method: sys_split_low_res.initialise(); sys_split_high_res.initialise(); Grad_t delF0_1; delF0_1 << 1.0, 0.0, 0.0, 1.0; constexpr Real cg_tol_low_res{1e-8}, newton_tol_low_res{1e-5}; constexpr Uint maxiter_low_res{CcoordOps::get_size(resolutions_split_low_res)*ipow(dim, secondOrder)*10}; constexpr bool verbose_low_res{false}; constexpr Real cg_tol_high_res{1e-8}, newton_tol_high_res{1e-5}; constexpr Uint maxiter_high_res{CcoordOps::get_size(resolutions_split_high_res)*ipow(dim, secondOrder)*10}; constexpr bool verbose_high_res{false}; - GradIncrements grads_low_res; grads_low_res.push_back(delF0_1); - GradIncrements grads_high_res; grads_high_res.push_back(delF0_1); - DeprecatedSolverCG cg_low_res{sys_split_low_res, cg_tol_low_res, maxiter_low_res, bool(verbose_low_res)}; - Eigen::ArrayXXd res_low_res{deprecated_newton_cg(sys_split_low_res, grads_low_res, cg_low_res, newton_tol_low_res, verbose_low_res)[0].grad}; + Eigen::ArrayXXd res_low_res{deprecated_newton_cg(sys_split_low_res, delF0_1, cg_low_res, newton_tol_low_res, verbose_low_res).grad}; DeprecatedSolverCGEigen cg_high_res{sys_split_high_res, cg_tol_high_res, maxiter_high_res, bool(verbose_high_res)}; - Eigen::ArrayXXd res_high_res{deprecated_newton_cg(sys_split_high_res, grads_high_res, cg_high_res, newton_tol_high_res, verbose_high_res)[0].grad}; + Eigen::ArrayXXd res_high_res{deprecated_newton_cg(sys_split_high_res, delF0_1, cg_high_res, newton_tol_high_res, verbose_high_res).grad}; - std::vector> result_low_res{}, result_high_res_filtered{}, result_comparison{} ; + /*std::vector> result_low_res{}, result_high_res_filtered{}, result_comparison{} ; std::vector temp_vec_low_res{}, temp_vec_high_res{}, temp_vec_result_comparison{}; for (int i{0} ; i < high_res ; i++){ if (i%5 == 0){ for(int j{0} ; j < high_res ; j++){ if(j%5 == 0){ temp_vec_high_res.push_back(res_high_res.col(i* high_res +j)); } } result_high_res_filtered.push_back(temp_vec_high_res); + temp_vec_high_res.clear(); } } for (int i{0} ; i < low_res ; i++){ for(int j{0} ; j < low_res ; j++){ temp_vec_low_res.push_back(res_low_res.col(i* low_res +j)); } result_low_res.push_back(temp_vec_low_res); + temp_vec_low_res.clear(); } for (unsigned int i{0} ; i < result_low_res.size(); i++){ for (unsigned int j{0} ; j < result_low_res.size(); j++){ temp_vec_result_comparison.push_back(result_high_res_filtered[i][j]-result_low_res[i][j]); } result_comparison.push_back(temp_vec_result_comparison); - } + temp_vec_result_comparison.clear(); + }*/ BOOST_CHECK_LE(abs(0),cg_tol_low_res); } BOOST_AUTO_TEST_SUITE_END(); } diff --git a/tests/split_test_intersection_sym_check.cc b/tests/split_test_intersection_sym_check.cc index 332e446..fb2dcb9 100644 --- a/tests/split_test_intersection_sym_check.cc +++ b/tests/split_test_intersection_sym_check.cc @@ -1,190 +1,191 @@ /** * @file test_intersection_error_induced.cc * * @author Ali Faslfi * * @date 21 Jun 2018 * * @brief Tests for split cells and octree material assignment * * Copyright © 2017 Till Ali Faslafi * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "tests.hh" #include "solver/deprecated_solvers.hh" #include "solver/deprecated_solver_cg.hh" #include "solver/deprecated_solver_cg_eigen.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain_fast.hh" #include "materials/material_linear_elastic1.hh" #include "common/iterators.hh" #include "common/ccoord_operations.hh" #include "common/common.hh" #include "cell/cell_factory.hh" #include "cell/cell_split.hh" #include "common/intersection_octree.hh" #include #include #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(symmetric_test); BOOST_AUTO_TEST_CASE(checking_the_symmetry) { constexpr Dim_t dim{twoD}; using Rcoord = Rcoord_t; using Ccoord = Ccoord_t; //typedef std::vector::iterator Array_iter; using Mat_t = MaterialLinearElastic1; const Real contrast {10}; const Real Young_soft{1.0030648180242636}, Poisson_soft{0.29930675909878679}; const Real Young_hard{contrast * Young_soft}, Poisson_hard{0.29930675909878679}; constexpr Real length {25}; constexpr Rcoord center{12.5, 12.5}; constexpr Rcoord halfside_neg{-3.1, -3.1}; constexpr Rcoord halfside_pos{+4.00, +4.00}; constexpr Real cg_tol{1e-8}, newton_tol{1e-5}; constexpr bool verbose{false}; std::vector precipitate_vertices ; precipitate_vertices.push_back ({center[0] + halfside_neg[0], center[1] + halfside_neg[1]}); precipitate_vertices.push_back ({center[0] + halfside_neg[0], center[1] + halfside_pos[1]}); precipitate_vertices.push_back ({center[0] + halfside_pos[0], center[1] + halfside_neg[1]}); precipitate_vertices.push_back ({center[0] + halfside_pos[0], center[1] + halfside_pos[1]}); Grad_t delF0_1; delF0_1 << 0.002, 0.0, 0.001, 0.005; ////LOW RESOLUTION///////////////////// constexpr int resolution_l{25}; constexpr Rcoord lengths_split_l{ length, length}; constexpr Ccoord resolutions_split_l{resolution_l, resolution_l}; auto fft_ptr_split_l{std::make_unique>(resolutions_split_l, ipow(dim, 2))}; auto proj_ptr_split_l{std::make_unique>(std::move(fft_ptr_split_l), lengths_split_l)}; CellSplit sys_split_l(std::move(proj_ptr_split_l), SplittedCell::yes); - auto& Material_hard_split_l = Mat_t::make(sys_split_l, "hard", Young_hard, Poisson_hard); auto& Material_soft_split_l = Mat_t::make(sys_split_l, "soft", Young_soft, Poisson_soft); - + sys_split_l.make_automatic_precipitate(precipitate_vertices, Material_hard_split_l); + sys_split_l.complete_material_assignment(Material_soft_split_l); //assign precipitate materials: - RootNode precipitate_l (sys_split_l, precipitate_vertices); + /* RootNode precipitate_l (sys_split_l, precipitate_vertices); //Extracting the intersected pixels and their correspondent intersection ratios: auto precipitate_intersects_l = precipitate_l.get_intersected_pixels(); auto precipitate_intersection_ratios_l = precipitate_l.get_intersection_ratios(); Material_hard_split_l.add_split_pixels_precipitate(precipitate_intersects_l, - precipitate_intersection_ratios_l); + precipitate_intersection_ratios_l);*/ //std::vector assigned_ratio_l = sys_split_l.get_assigned_ratios(); - auto incompleted_pixels_l = sys_split_l.make_incompleted_pixels(); + /*auto incompleted_pixels_l = sys_split_l.make_incompleted_pixels(); for (auto && tup:incompleted_pixels_l){ auto && pix = std::get<0> (tup); auto && ratio = std::get<1> (tup); Material_soft_split_l.add_pixel_split(pix, ratio); } - /* for (auto && tup: akantu::enumerate(sys_split_l)) { + for (auto && tup: akantu::enumerate(sys_split_l)) { auto && pixel = std::get<1>(tup); auto iterator = std::get<0>(tup); if (assigned_ratio_l[iterator] < 1.0){ Material_soft_split_l.add_pixel_split(pixel, 1.0-assigned_ratio_l[iterator]); } }*/ ////HIGH RESOLUTION///////////////////// constexpr int resolution_h{125}; constexpr Rcoord lengths_split_h{ length, length}; constexpr Ccoord resolutions_split_h{resolution_h, resolution_h}; auto fft_ptr_split_h{std::make_unique>(resolutions_split_h, ipow(dim, 2))}; auto proj_ptr_split_h{std::make_unique>(std::move(fft_ptr_split_h), lengths_split_h)}; CellSplit sys_split_h(std::move(proj_ptr_split_h), SplittedCell::yes); auto& Material_hard_split_h = Mat_t::make(sys_split_h, "hard", Young_hard, Poisson_hard); auto& Material_soft_split_h = Mat_t::make(sys_split_h, "soft", Young_soft, Poisson_soft); + sys_split_h.make_automatic_precipitate(precipitate_vertices, Material_hard_split_h); + sys_split_h.complete_material_assignment(Material_soft_split_h); //assign precipitate materials: - RootNode precipitate_h (sys_split_h, precipitate_vertices); + /*RootNode precipitate_h (sys_split_h, precipitate_vertices); //Extracting the intersected pixels and their correspondent intersection ratios: auto precipitate_intersects_h= precipitate_h.get_intersected_pixels(); auto precipitate_intersection_ratios_h = precipitate_h.get_intersection_ratios(); Material_hard_split_h.add_split_pixels_precipitate(precipitate_intersects_h, precipitate_intersection_ratios_h); //std::vector assigned_ratio_h = sys_split_h.get_assigned_ratios(); - /* + auto incompleted_pixels_h = sys_split_h.make_incompleted_pixels(); for (auto && tup:incompleted_pixels_h){ auto && pix = std::get<0> (tup); auto && ratio = std::get<1> (tup); Material_soft_split_h.add_pixel_split(pix, ratio); - } - */ - sys_split_h.complete_material_assignment(Material_soft_split_h); + }*/ + /* for (auto && tup: akantu::enumerate(sys_split_h)) { auto && pixel = std::get<1>(tup); auto iterator = std::get<0>(tup); if (assigned_ratio_h[iterator] < 1.0){ Material_soft_split_h.add_pixel_split(pixel, 1.0-assigned_ratio_h[iterator]); } }*/ //Finding equilbrium strain fileds with Sppectral method: sys_split_l.initialise(); sys_split_h.initialise(); constexpr Uint maxiter_l{CcoordOps::get_size(resolutions_split_l)*ipow(dim, secondOrder)*10}; constexpr Uint maxiter_h{CcoordOps::get_size(resolutions_split_h)*ipow(dim, secondOrder)*10}; DeprecatedSolverCG cg_l{sys_split_l, cg_tol, maxiter_l, bool(verbose)}; auto output_l {deprecated_newton_cg(sys_split_l, delF0_1, cg_l, newton_tol, verbose)}; Eigen::ArrayXXd res_grad_l{output_l.grad}; Eigen::ArrayXXd res_P_l{output_l.stress}; DeprecatedSolverCG cg_h{sys_split_h, cg_tol, maxiter_h, bool(verbose)}; auto output_h {deprecated_newton_cg(sys_split_h, delF0_1, cg_h, newton_tol, verbose)}; Eigen::ArrayXXd res_grad_h{output_h.grad}; Eigen::ArrayXXd res_P_h{output_h.stress}; std::ofstream file_grad_l("grad_l.csv"); std::ofstream file_P_l("P_l.csv"); for (int j{0}; j < res_grad_l.cols(); j++){ file_grad_l << res_grad_l.col(j).transpose()<<"\n"; file_P_l << res_P_l.col(j).transpose()<<"\n"; } std::ofstream file_grad_h("grad_h.csv"); std::ofstream file_P_h("P_h.csv"); for (int j{0}; j < res_grad_h.cols(); j++){ file_grad_h << res_grad_h.col(j).transpose()<<"\n"; file_P_h << res_P_h.col(j).transpose()<<"\n"; } BOOST_CHECK_LE(abs(0),cg_tol); } BOOST_AUTO_TEST_SUITE_END(); }