diff --git a/src/common/field_collection.hh b/src/common/field_collection.hh index ce5f92b..82132f9 100644 --- a/src/common/field_collection.hh +++ b/src/common/field_collection.hh @@ -1,50 +1,35 @@ /** * @file field_collection.hh * * @author Till Junge * * @date 07 Sep 2017 * * @brief Provides pixel-iterable containers for scalar and tensorial fields, * addressable by field name * * 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_H #define FIELD_COLLECTION_H #include "common/field_collection_global.hh" #include "common/field_collection_local.hh" -namespace muSpectre { - - /** - * Convenience aliases to unify the local and global - * Field_collection interface - */ - template - using FieldCollection = std::conditional_t< - Global, - GlobalFieldCollection, - LocalFieldCollection>; - - -} // muSpectre - #endif /* FIELD_COLLECTION_H */ diff --git a/src/fft/fft_engine_base.hh b/src/fft/fft_engine_base.hh index 7d64f89..44b47b3 100644 --- a/src/fft/fft_engine_base.hh +++ b/src/fft/fft_engine_base.hh @@ -1,145 +1,145 @@ /** * @file fft_engine_base.hh * * @author Till Junge * * @date 01 Dec 2017 * * @brief Interface for FFT engines * * 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 FFT_ENGINE_BASE_H #define FFT_ENGINE_BASE_H #include "common/common.hh" #include "common/field_collection.hh" namespace muSpectre { /** * Virtual base class for FFT engines. To be implemented by all * FFT_engine implementations. */ template class FFT_Engine_base { public: constexpr static Dim_t sdim{DimS}; //!< spatial dimension of the cell constexpr static Dim_t mdim{DimM}; //!< material dimension of the cell //! cell coordinates type using Ccoord = Ccoord_t; //! spatial coordinates type using Rcoord = std::array; //! global FieldCollection - using GFieldCollection_t = FieldCollection; + using GFieldCollection_t = GlobalFieldCollection; //! local FieldCollection (for Fourier-space pixels) - using LFieldCollection_t = FieldCollection; + using LFieldCollection_t = LocalFieldCollection; //! Field type on which to apply the projection using Field_t = TensorField; /** * Field type holding a Fourier-space representation of a * real-valued second-order tensor field */ using Workspace_t = TensorField; /** * iterator over Fourier-space discretisation point */ using iterator = typename LFieldCollection_t::iterator; //! Default constructor FFT_Engine_base() = delete; //! Constructor with system resolutions FFT_Engine_base(Ccoord resolutions, Rcoord lengths); //! Copy constructor FFT_Engine_base(const FFT_Engine_base &other) = delete; //! Move constructor FFT_Engine_base(FFT_Engine_base &&other) = default; //! Destructor virtual ~FFT_Engine_base() = default; //! Copy assignment operator FFT_Engine_base& operator=(const FFT_Engine_base &other) = delete; //! Move assignment operator FFT_Engine_base& operator=(FFT_Engine_base &&other) = default; //! compute the plan, etc virtual void initialise(FFT_PlanFlags /*plan_flags*/); //! forward transform (dummy for interface) virtual Workspace_t & fft(Field_t & /*field*/) = 0; //! inverse transform (dummy for interface) virtual void ifft(Field_t & /*field*/) const = 0; /** * iterators over only those pixels that exist in frequency space * (i.e. about half of all pixels, see rfft) */ //! returns an iterator to the first pixel in Fourier space inline iterator begin() {return this->work_space_container.begin();} //! returns an iterator past to the last pixel in Fourier space inline iterator end() {return this->work_space_container.end();} //! nb of pixels (mostly for debugging) size_t size() const; //! nb of pixels in Fourier space size_t workspace_size() const; //! returns the resolutions of the cell const Ccoord & get_resolutions() const {return this->resolutions;} //! returns the physical sizes of the cell const Rcoord & get_lengths() const {return this->lengths;} //! only required for testing and debugging LFieldCollection_t & get_field_collection() { return this->work_space_container;} //! only required for testing and debugging Workspace_t& get_work_space() {return this->work;} //! factor by which to multiply projection before inverse transform (this is //! typically 1/nb_pixels for so-called unnormalized transforms (see, //! e.g. http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data //! or https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html //! . Rather than scaling the inverse transform (which would cost one more //! loop), FFT engines provide this value so it can be used in the //! projection operator (where no additional loop is required) inline Real normalisation() const {return norm_factor;}; protected: /** * Field collection in which to store fields associated with * Fourier-space points */ LFieldCollection_t work_space_container{}; const Ccoord resolutions; //!< resolutions of the cell const Rcoord lengths; //!< physical sizes of the cell Workspace_t & work; //!< field to store the Fourier transform of P const Real norm_factor; //!< normalisation coefficient of fourier transform private: }; } // muSpectre #endif /* FFT_ENGINE_BASE_H */ diff --git a/src/fft/projection_default.hh b/src/fft/projection_default.hh index f45bfa4..7b02966 100644 --- a/src/fft/projection_default.hh +++ b/src/fft/projection_default.hh @@ -1,98 +1,98 @@ /** * @file projection_default.hh * * @author Till Junge * * @date 14 Jan 2018 * * @brief virtual base class for default projection implementation, where the * projection operator is stored as a full fourth-order tensor per * k-space point (as opposed to 'smart' faster implementations, such as * ProjectionFiniteStrainFast * * Copyright (C) 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef PROJECTION_DEFAULT_H #define PROJECTION_DEFAULT_H #include "fft/projection_base.hh" namespace muSpectre { /** * base class to inherit from if one implements a projection * operator that is stored in form of a fourth-order tensor of real * values per k-grid point */ template class ProjectionDefault: public ProjectionBase { public: using Parent = ProjectionBase; //!< base class //! polymorphic FFT pointer type using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type //! global field collection - using GFieldCollection_t = FieldCollection; + using GFieldCollection_t = GlobalFieldCollection; //! local field collection for Fourier-space fields - using LFieldCollection_t = FieldCollection; + using LFieldCollection_t = LocalFieldCollection; //! Real space second order tensor fields (to be projected) using Field_t = TensorField; //! Fourier-space field containing the projection operator itself using Proj_t = TensorField; //! iterable form of the operator using Proj_map = T4MatrixFieldMap; //! vectorized version of the Fourier-space second-order tensor field using Vector_map = MatrixFieldMap; //! Default constructor ProjectionDefault() = delete; //! Constructor with system sizes and formulation ProjectionDefault(FFT_Engine_ptr engine, Formulation form); //! Copy constructor ProjectionDefault(const ProjectionDefault &other) = delete; //! Move constructor ProjectionDefault(ProjectionDefault &&other) = default; //! Destructor virtual ~ProjectionDefault() = default; //! Copy assignment operator ProjectionDefault& operator=(const ProjectionDefault &other) = delete; //! Move assignment operator ProjectionDefault& operator=(ProjectionDefault &&other) = delete; //! apply the projection operator to a field void apply_projection(Field_t & field) override final; Eigen::Map get_operator() override final; protected: Proj_t & Gfield; //!< field holding the operator Proj_map Ghat; //!< iterable version of operator private: }; } // muSpectre #endif /* PROJECTION_DEFAULT_H */ diff --git a/src/fft/projection_finite_strain.hh b/src/fft/projection_finite_strain.hh index 05ce41a..7dc7a4b 100644 --- a/src/fft/projection_finite_strain.hh +++ b/src/fft/projection_finite_strain.hh @@ -1,90 +1,90 @@ /** * @file projection_finite_strain.hh * * @author Till Junge * * @date 05 Dec 2017 * * @brief Class for standard finite-strain gradient projections see de Geus et * al. (https://doi.org/10.1016/j.cma.2016.12.032) for derivation * * 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 PROJECTION_FINITE_STRAIN_H #define PROJECTION_FINITE_STRAIN_H #include "fft/projection_default.hh" #include "common/common.hh" #include "common/field_collection.hh" #include "common/field_map.hh" namespace muSpectre { /** * Implements the finite strain gradient projection operator as * defined in de Geus et * al. (https://doi.org/10.1016/j.cma.2016.12.032) for derivation */ template class ProjectionFiniteStrain: public ProjectionDefault { public: using Parent = ProjectionDefault; //!< base class //! polymorphic pointer to FFT engines using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type //! local field collection (for Fourier-space representations) - using LFieldCollection_t = FieldCollection; + using LFieldCollection_t = LocalFieldCollection; //! iterable operator using Proj_map = T4MatrixFieldMap; //! iterable vectorised version of the Fourier-space tensor field using Vector_map = MatrixFieldMap; //! Default constructor ProjectionFiniteStrain() = delete; //! Constructor with fft_engine ProjectionFiniteStrain(FFT_Engine_ptr engine); //! Copy constructor ProjectionFiniteStrain(const ProjectionFiniteStrain &other) = delete; //! Move constructor ProjectionFiniteStrain(ProjectionFiniteStrain &&other) = default; //! Destructor virtual ~ProjectionFiniteStrain() = default; //! Copy assignment operator ProjectionFiniteStrain& operator=(const ProjectionFiniteStrain &other) = delete; //! Move assignment operator ProjectionFiniteStrain& operator=(ProjectionFiniteStrain &&other) = default; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final; protected: private: }; } // muSpectre #endif /* PROJECTION_FINITE_STRAIN_H */ diff --git a/src/fft/projection_finite_strain_fast.hh b/src/fft/projection_finite_strain_fast.hh index 4575211..332cf9c 100644 --- a/src/fft/projection_finite_strain_fast.hh +++ b/src/fft/projection_finite_strain_fast.hh @@ -1,105 +1,105 @@ /** * @file projection_finite_strain_fast.hh * * @author Till Junge * * @date 12 Dec 2017 * * @brief Faster alternative to ProjectionFinitestrain * * 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 PROJECTION_FINITE_STRAIN_FAST_H #define PROJECTION_FINITE_STRAIN_FAST_H #include "fft/projection_base.hh" #include "common/common.hh" #include "common/field_collection.hh" #include "common/field_map.hh" namespace muSpectre { /** * replaces `muSpectre::ProjectionFiniteStrain` with a faster and * less memory-hungry alternative formulation. Use this if you don't * have a very good reason not to (and tell me (author) about it, * I'd be interested to hear it). */ template class ProjectionFiniteStrainFast: public ProjectionBase { public: using Parent = ProjectionBase; //!< base class //! polymorphic pointer to FFT engines using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type //! global field collection (for real-space representations) - using GFieldCollection_t = FieldCollection; + using GFieldCollection_t = GlobalFieldCollection; //! local field collection (for Fourier-space representations) - using LFieldCollection_t = FieldCollection; + using LFieldCollection_t = LocalFieldCollection; //! Real space second order tensor fields (to be projected) using Field_t = TensorField; //! Fourier-space field containing the projection operator itself using Proj_t = TensorField; //! iterable form of the operator using Proj_map = MatrixFieldMap; //! iterable Fourier-space second-order tensor field using Grad_map = MatrixFieldMap; //! Default constructor ProjectionFiniteStrainFast() = delete; //! Constructor with fft_engine ProjectionFiniteStrainFast(FFT_Engine_ptr engine); //! Copy constructor ProjectionFiniteStrainFast(const ProjectionFiniteStrainFast &other) = delete; //! Move constructor ProjectionFiniteStrainFast(ProjectionFiniteStrainFast &&other) = default; //! Destructor virtual ~ProjectionFiniteStrainFast() = default; //! Copy assignment operator ProjectionFiniteStrainFast& operator=(const ProjectionFiniteStrainFast &other) = delete; //! Move assignment operator ProjectionFiniteStrainFast& operator=(ProjectionFiniteStrainFast &&other) = default; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final; //! apply the projection operator to a field void apply_projection(Field_t & field) override final; Eigen::Map get_operator() override final; protected: Proj_t & xiField; //!< field of normalised wave vectors Proj_map xis; //!< iterable normalised wave vectors private: }; } // muSpectre #endif /* PROJECTION_FINITE_STRAIN_FAST_H */ diff --git a/src/fft/projection_small_strain.hh b/src/fft/projection_small_strain.hh index fd7ca7b..d872d61 100644 --- a/src/fft/projection_small_strain.hh +++ b/src/fft/projection_small_strain.hh @@ -1,93 +1,93 @@ /** * @file projection_small_strain.cc * * @author Till Junge * * @date 14 Jan 2018 * * @brief Small strain projection operator as defined in Appendix A1 of * DOI: 10.1002/nme.5481 ("A finite element perspective on nonlinear * FFT-based micromechanical simulations", Int. J. Numer. Meth. Engng * 2017; 111 :903–926) * * 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. */ #ifndef PROJECTION_SMALL_STRAIN_H #define PROJECTION_SMALL_STRAIN_H #include "fft/projection_default.hh" namespace muSpectre { /** * Implements the small strain projection operator as defined in * Appendix A1 of DOI: 10.1002/nme.5481 ("A finite element * perspective on nonlinear FFT-based micromechanical * simulations", Int. J. Numer. Meth. Engng 2017; 111 * :903–926) */ template class ProjectionSmallStrain: public ProjectionDefault { public: using Parent = ProjectionDefault; //!< base class //! polymorphic pointer to FFT engines using FFT_Engine_ptr = typename Parent::FFT_Engine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type //! local field collection (for Fourier-space representations) - using LFieldCollection_t = FieldCollection; + using LFieldCollection_t = LocalFieldCollection; //! Fourier-space field containing the projection operator itself using Proj_t = TensorField; //! iterable operator using Proj_map = T4MatrixFieldMap; //! iterable vectorised version of the Fourier-space tensor field using Vector_map = MatrixFieldMap; //! Default constructor ProjectionSmallStrain() = delete; //! Constructor with fft_engine ProjectionSmallStrain(FFT_Engine_ptr engine); //! Copy constructor ProjectionSmallStrain(const ProjectionSmallStrain &other) = delete; //! Move constructor ProjectionSmallStrain(ProjectionSmallStrain &&other) = default; //! Destructor virtual ~ProjectionSmallStrain() = default; //! Copy assignment operator ProjectionSmallStrain& operator=(const ProjectionSmallStrain &other) = delete; //! Move assignment operator ProjectionSmallStrain& operator=(ProjectionSmallStrain &&other) = delete; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final; protected: private: }; } // muSpectre #endif /* PROJECTION_SMALL_STRAIN_H */ diff --git a/src/materials/material_base.hh b/src/materials/material_base.hh index 47697fb..516c4c4 100644 --- a/src/materials/material_base.hh +++ b/src/materials/material_base.hh @@ -1,148 +1,148 @@ /** * @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) template class MaterialBase { public: //! typedefs for data handled by this interface //! global field collection for system-wide fields, like stress, strain, etc - using GFieldCollection_t = FieldCollection; + 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 = FieldCollection; + 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); //! allocate memory, etc, but also: wipe history variables! virtual void initialise(bool stiffness = false) = 0; //! return the materil'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) = 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); //! computes stress and tangent moduli virtual void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) = 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); //! 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 private: }; } // muSpectre #endif /* MATERIAL_BASE_H */ diff --git a/src/system/system_base.hh b/src/system/system_base.hh index 90c2f37..7d8ff0f 100644 --- a/src/system/system_base.hh +++ b/src/system/system_base.hh @@ -1,293 +1,293 @@ /** * @file system_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell system 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 SYSTEM_BASE_H #define SYSTEM_BASE_H #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "system/system_traits.hh" #include #include #include #include namespace muSpectre { template class SystemAdaptor; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class SystemBase { public: using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< physical coordinates type //! global field collection - using FieldCollection_t = FieldCollection; + 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; //! input vector for solvers using SolvVectorIn = Eigen::Ref; //! output vector for solvers using SolvVectorOut = Eigen::Map; //! sparse matrix emulation using Adaptor = SystemAdaptor; //! Default constructor SystemBase() = delete; //! constructor using sizes and resolution SystemBase(Projection_ptr projection); //! Copy constructor SystemBase(const SystemBase &other) = delete; //! Move constructor SystemBase(SystemBase &&other) = default; //! Destructor virtual ~SystemBase() = default; //! Copy assignment operator SystemBase& operator=(const SystemBase &other) = delete; //! Move assignment operator SystemBase& operator=(SystemBase &&other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this system */ Material_t & add_material(Material_ptr mat); /** * evaluates all materials */ FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF or G:K:δε) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * vectorized version for eigen solvers, no copy, but only works * when fields have ArrayStore=false */ SolvVectorOut directional_stiffness_vec(const SolvVectorIn & delF); /** * Evaluate directional stiffness into a temporary array and * return a copy. This is a costly and wasteful interface to * directional_stiffness and should only be used for debugging or * in the python interface */ Eigen::ArrayXXd directional_stiffness_with_copy (Eigen::Ref delF); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); //! 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 system StrainField_t & get_managed_field(std::string unique_name); /** * general initialisation; initialises the projection and * fft_engine (i.e. infrastructure) but not the materials. These * need to be initialised separately */ void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); /** * initialise materials (including resetting any history variables) */ void initialise_materials(bool stiffness=false); iterator begin(); //!< iterator 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 resolutions of the cell const Ccoord & get_resolutions() const {return this->resolutions;} //! return the sizes of the cell const Rcoord & get_lengths() const {return this->lengths;} /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const { return this->projection->get_formulation();} //! for handling double initialisations right bool is_initialised() const {return this->initialised;} /** * 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(); //! returns the number of degrees of freedom in the cell Dim_t nb_dof() const {return this->size()*ipow(DimS, 2);}; protected: //! make sure that every pixel is assigned to one and only one material void check_material_coverage(); const Ccoord & resolutions; //!< the cell's resolutions CcoordOps::Pixels pixels; //!< helper to iterate over the pixels const Rcoord & 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 const Formulation form; //!< formulation for solution private: }; /** * lightweight resource handle wrapping a `muSpectre::SystemBase` or * a subclass thereof into `Eigen::EigenBase`, so it can be * interpreted as a sparse matrix by Eigen solvers */ template class SystemAdaptor: 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 SystemAdaptor(System & sys):sys{sys}{} //!returns the number of logical rows Eigen::Index rows() const {return this->sys.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()); } System & sys; //!< ref to the cell }; } // muSpectre namespace Eigen { namespace internal { //! Implementation of `muSpectre::SystemAdaptor` * `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 SystemAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/) { // This method should implement "dst += alpha * lhs * rhs" inplace, // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, dst.noalias() += const_cast(lhs).sys.directional_stiffness_vec(rhs); } }; } } #endif /* SYSTEM_BASE_H */ diff --git a/tests/test_fftw_engine.cc b/tests/test_fftw_engine.cc index a0bab4a..d208c37 100644 --- a/tests/test_fftw_engine.cc +++ b/tests/test_fftw_engine.cc @@ -1,129 +1,129 @@ /** * @file test_fftw_engine.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief tests for the fftw fft engine implementation * * 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 #include "tests.hh" #include "fft/fftw_engine.hh" #include "common/ccoord_operations.hh" #include "common/field_collection.hh" #include "common/field_map.hh" #include "common/iterators.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(fftw_engine); /* ---------------------------------------------------------------------- */ template struct FFTW_fixture { constexpr static Dim_t box_resolution{resolution}; constexpr static Real box_length{4.5}; constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; constexpr static Ccoord_t res() { return CcoordOps::get_cube(box_resolution); } FFTW_fixture() :engine(res(), CcoordOps::get_cube(box_length)){} FFTW_Engine engine; }; struct FFTW_fixture_python_segfault{ constexpr static Dim_t dim{twoD}; constexpr static Dim_t sdim{twoD}; constexpr static Dim_t mdim{twoD}; constexpr static Ccoord_t res() {return {6, 4};} FFTW_fixture_python_segfault():engine{res(), {3., 3}} {} FFTW_Engine engine; }; using fixlist = boost::mpl::list, FFTW_fixture< twoD, threeD, 3>, FFTW_fixture, FFTW_fixture< twoD, twoD, 4>, FFTW_fixture< twoD, threeD, 4>, FFTW_fixture, FFTW_fixture_python_segfault>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Constructor_test, Fix, fixlist, Fix) { BOOST_CHECK_NO_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate)); BOOST_CHECK_EQUAL(Fix::engine.size(), CcoordOps::get_size(Fix::res())); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(fft_test, Fix, fixlist, Fix) { Fix::engine.initialise(FFT_PlanFlags::estimate); constexpr Dim_t order{2}; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t fc; auto & input{make_field>("input", fc)}; auto & ref {make_field>("reference", fc)}; auto & result{make_field>("result", fc)}; fc.initialise(Fix::res()); using map_t = MatrixFieldMap; map_t inmap{input}; auto refmap{map_t{ref}}; auto resultmap{map_t{result}}; size_t cntr{0}; for (auto tup: akantu::zip(inmap, refmap)) { cntr++; auto & in_{std::get<0>(tup)}; auto & ref_{std::get<1>(tup)}; in_.setRandom(); ref_ = in_; } auto & complex_field = Fix::engine.fft(input); - using cmap_t = MatrixFieldMap, Complex, Fix::mdim, Fix::mdim>; + using cmap_t = MatrixFieldMap, Complex, Fix::mdim, Fix::mdim>; cmap_t complex_map(complex_field); Real error = complex_map[0].imag().norm(); BOOST_CHECK_LT(error, tol); /* make sure, the engine has not modified input (which is unfortunately const-casted internally, hence this test) */ for (auto && tup: akantu::zip(inmap, refmap)) { Real error{(std::get<0>(tup) - std::get<1>(tup)).norm()}; BOOST_CHECK_LT(error, tol); } /* make sure that the ifft of fft returns the original*/ Fix::engine.ifft(result); for (auto && tup: akantu::zip(resultmap, refmap)) { Real error{(std::get<0>(tup)*Fix::engine.normalisation() - std::get<1>(tup)).norm()}; BOOST_CHECK_LT(error, tol); if (error > tol) { std::cout << std::get<0>(tup).array()/std::get<1>(tup).array() << std::endl << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_field_collections.cc b/tests/test_field_collections.cc index f4be509..c1e9017 100644 --- a/tests/test_field_collections.cc +++ b/tests/test_field_collections.cc @@ -1,221 +1,221 @@ /** * @file test_field_collections.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test the FieldCollection classes which provide fast optimized iterators * over run-time typed 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. */ #include "test_field_collections_header.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_collection_tests); BOOST_AUTO_TEST_CASE(simple) { const Dim_t sdim = 2; const Dim_t mdim = 2; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t fc; BOOST_CHECK_EQUAL(FC_t::spatial_dim(), sdim); BOOST_CHECK_EQUAL(fc.get_spatial_dim(), sdim); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(Simple_construction_test, F, test_collections, F) { BOOST_CHECK_EQUAL(F::FC_t::spatial_dim(), F::sdim()); BOOST_CHECK_EQUAL(F::fc.get_spatial_dim(), F::sdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(get_field2_test, F, test_collections, F) { const auto order{2}; using FC_t = typename F::FC_t; using TF_t = TensorField; auto && myfield = make_field("TensorField real 2", F::fc); using TensorMap = TensorFieldMap; using MatrixMap = MatrixFieldMap; using ArrayMap = ArrayFieldMap; TensorMap TFM(myfield); MatrixMap MFM(myfield); ArrayMap AFM(myfield); BOOST_CHECK_EQUAL(TFM.info_string(), "Tensor(d, "+ std::to_string(order) + "_o, " + std::to_string(F::mdim()) + "_d)"); BOOST_CHECK_EQUAL(MFM.info_string(), "Matrix(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); BOOST_CHECK_EQUAL(AFM.info_string(), "Array(d, "+ std::to_string(F::mdim()) + "x" + std::to_string(F::mdim()) + ")"); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(multi_field_test, F, mult_collections, F) { using FC_t = typename F::FC_t; // possible maptypes for Real tensor fields using T_type = Real; using T_TFM1_t = TensorFieldMap; using T_TFM2_t = TensorFieldMap; //! dangerous using T4_Map_t = T4MatrixFieldMap; // impossible maptypes for Real tensor fields using T_SFM_t = ScalarFieldMap; using T_MFM_t = MatrixFieldMap; using T_AFM_t = ArrayFieldMap; using T_MFMw1_t = MatrixFieldMap; using T_MFMw2_t = MatrixFieldMap; using T_MFMw3_t = MatrixFieldMap; const std::string T_name{"Tensorfield Real o4"}; const std::string T_name_w{"TensorField Real o4 wrongname"}; BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(T_TFM1_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T_TFM2_t(F::fc.at(T_name))); BOOST_CHECK_NO_THROW(T4_Map_t(F::fc.at(T_name))); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(T_name_w)), std::out_of_range); BOOST_CHECK_THROW(T_MFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_AFM_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw1_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw2_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_MFMw3_t(F::fc.at(T_name)), FieldInterpretationError); BOOST_CHECK_THROW(T_SFM_t(F::fc.at(T_name_w)), std::out_of_range); // possible maptypes for integer scalar fields using S_type = Int; using S_SFM_t = ScalarFieldMap; using S_TFM1_t = TensorFieldMap; using S_TFM2_t = TensorFieldMap; using S_MFM_t = MatrixFieldMap; using S_AFM_t = ArrayFieldMap; using S4_Map_t = T4MatrixFieldMap; // impossible maptypes for integer scalar fields using S_MFMw1_t = MatrixFieldMap; using S_MFMw2_t = MatrixFieldMap; using S_MFMw3_t = MatrixFieldMap; const std::string S_name{"integer Scalar"}; const std::string S_name_w{"integer Scalar wrongname"}; BOOST_CHECK_NO_THROW(S_SFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM1_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_TFM2_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_MFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S_AFM_t(F::fc.at(S_name))); BOOST_CHECK_NO_THROW(S4_Map_t(F::fc.at(S_name))); BOOST_CHECK_THROW(S_MFMw1_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(T4_Map_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw2_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_MFMw3_t(F::fc.at(S_name)), FieldInterpretationError); BOOST_CHECK_THROW(S_SFM_t(F::fc.at(S_name_w)), std::out_of_range); // possible maptypes for complex matrix fields using M_type = Complex; using M_MFM_t = MatrixFieldMap; using M_AFM_t = ArrayFieldMap; // impossible maptypes for complex matrix fields using M_SFM_t = ScalarFieldMap; using M_MFMw1_t = MatrixFieldMap; using M_MFMw2_t = MatrixFieldMap; using M_MFMw3_t = MatrixFieldMap; const std::string M_name{"Matrixfield Complex sdim x mdim"}; const std::string M_name_w{"Matrixfield Complex sdim x mdim wrongname"}; BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_NO_THROW(M_MFM_t(F::fc.at(M_name))); BOOST_CHECK_NO_THROW(M_AFM_t(F::fc.at(M_name))); BOOST_CHECK_THROW(M_MFMw1_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw2_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_MFMw3_t(F::fc.at(M_name)), FieldInterpretationError); BOOST_CHECK_THROW(M_SFM_t(F::fc.at(M_name_w)), std::out_of_range); } /* ---------------------------------------------------------------------- */ //! Check whether fields can be initialized using mult_collections_t = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>>; using mult_collections_f = boost::mpl::list, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_glob, F, mult_collections_t, F) { Ccoord_t size; for (auto && s: size) { s = 3; } BOOST_CHECK_NO_THROW(F::fc.initialise(size)); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca, F, mult_collections_f, F) { testGoodies::RandRange rng; for (int i = 0; i < 7; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(init_test_loca_with_push_back, F, mult_collections_f, F) { constexpr auto mdim = F::mdim(); constexpr int nb_pix = 7; testGoodies::RandRange rng; using ftype = internal::TypedSizedFieldBase< decltype(F::fc), Real, mdim*mdim*mdim*mdim>; using stype = Eigen::Array; auto & field = reinterpret_cast(F::fc["Tensorfield Real o4"]); field.push_back(stype()); for (int i = 0; i < nb_pix; ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } F::fc.add_pixel(pixel); } BOOST_CHECK_THROW(F::fc.initialise(), FieldCollectionError); for (int i = 0; i < nb_pix-1; ++i) { field.push_back(stype()); } BOOST_CHECK_NO_THROW(F::fc.initialise()); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_field_collections_header.hh b/tests/test_field_collections_header.hh index 695e8b0..4a1893e 100644 --- a/tests/test_field_collections_header.hh +++ b/tests/test_field_collections_header.hh @@ -1,156 +1,162 @@ /** * @file test_field_collections_header.hh * * @author Till Junge * * @date 23 Nov 2017 * * @brief declares fixtures for field_collection tests, so that they can be split * * 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 TEST_FIELD_COLLECTIONS_HEADER_H #define TEST_FIELD_COLLECTIONS_HEADER_H #include #include #include #include #include #include #include "common/common.hh" #include "common/ccoord_operations.hh" #include "tests/test_goodies.hh" #include "tests.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common/field_map.hh" namespace muSpectre { //! Test fixture for simple tests on single field in collection template struct FC_fixture: - public FieldCollection { + public std::conditional_t, + LocalFieldCollection> { FC_fixture() :fc() {} inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} - using FC_t = FieldCollection; + using FC_t = std::conditional_t, + LocalFieldCollection>; FC_t fc; }; using test_collections = boost::mpl::list, FC_fixture<2, 3, true>, FC_fixture<3, 3, true>, FC_fixture<2, 2, false>, FC_fixture<2, 3, false>, FC_fixture<3, 3, false>>; constexpr Dim_t order{4}, matrix_order{2}; //! Test fixture for multiple fields in the collection template struct FC_multi_fixture{ - using FC_t = FieldCollection; + using FC_t = std::conditional_t, + LocalFieldCollection>; using T4_t = TensorField; using T2_t = TensorField; using Sc_t = ScalarField; using M2_t = MatrixField; FC_multi_fixture() :fc(), t4_field{make_field("Tensorfield Real o4", fc)},//Real tensor field t2_field{make_field("Tensorfield Real o2", fc)},//Real tensor field sc_field{make_field("integer Scalar", fc)}, // integer scalar field m2_field{make_field("Matrixfield Complex sdim x mdim", fc)} //complex matrix field { } inline static constexpr Dim_t sdim(){return DimS;} inline static constexpr Dim_t mdim(){return DimM;} inline static constexpr bool global(){return Global;} FC_t fc; T4_t & t4_field; T2_t & t2_field; Sc_t & sc_field; M2_t & m2_field; }; using mult_collections = boost::mpl::list, FC_multi_fixture<2, 3, true>, FC_multi_fixture<3, 3, true>, FC_multi_fixture<2, 2, false>, FC_multi_fixture<2, 3, false>, FC_multi_fixture<3, 3, false>>; //! Test fixture for iterators over multiple fields template struct FC_iterator_fixture : public FC_multi_fixture { using Parent = FC_multi_fixture; FC_iterator_fixture() :Parent() { this-> fill(); } template std::enable_if_t fill() { static_assert(Global==isGlobal, "You're breaking my SFINAE plan"); Ccoord_t size; for (auto && s: size) { s = cube_size(); } this->fc.initialise(size); } template std::enable_if_t fill (int dummy = 0) { static_assert(notGlobal != Global, "You're breaking my SFINAE plan"); testGoodies::RandRange rng; this->fc.add_pixel({0,0}); for (int i = 0*dummy; i < sele_size(); ++i) { Ccoord_t pixel; for (auto && s: pixel) { s = rng.randval(0, 7); } this->fc.add_pixel(pixel); } this->fc.initialise(); } constexpr static Dim_t cube_size() {return 3;} constexpr static Dim_t sele_size() {return 7;} }; using iter_collections = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>, FC_iterator_fixture<2, 2, false>, FC_iterator_fixture<2, 3, false>, FC_iterator_fixture<3, 3, false>>; using glob_iter_colls = boost::mpl::list, FC_iterator_fixture<2, 3, true>, FC_iterator_fixture<3, 3, true>>; } // muSpectre #endif /* TEST_FIELD_COLLECTIONS_HEADER_H */ diff --git a/tests/test_fields.cc b/tests/test_fields.cc index 88f4ef9..e97f746 100644 --- a/tests/test_fields.cc +++ b/tests/test_fields.cc @@ -1,49 +1,49 @@ /** * @file test_fields.cc * * @author Till Junge * * @date 20 Sep 2017 * * @brief Test Fields that are used in FieldCollections * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "tests.hh" #include "common/field_collection.hh" #include "common/field.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(field_test); BOOST_AUTO_TEST_CASE(simple_creation) { const Dim_t sdim = 2; const Dim_t mdim = 2; const Dim_t order = 4; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t fc; using TF_t = TensorField; make_field("TensorField 1", fc); } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_material_linear_elastic1.cc b/tests/test_material_linear_elastic1.cc index 8ce1b0b..339c30c 100644 --- a/tests/test_material_linear_elastic1.cc +++ b/tests/test_material_linear_elastic1.cc @@ -1,230 +1,230 @@ /** * @file test_material_linear_elastic1.cc * * @author Till Junge * * @date 28 Nov 2017 * * @brief Tests for the large-strain, objective Hooke's law, implemented in * the convenient strategy (i.e., using MaterialMuSpectre), also used * to test parts of MaterialLinearElastic2 * * 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 #include #include #include "materials/material_linear_elastic1.hh" #include "materials/material_linear_elastic2.hh" #include "tests.hh" #include "tests/test_goodies.hh" #include "common/field_collection.hh" #include "common/iterators.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_1); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real lambda{2}, mu{1.5}; constexpr static Real young{mu*(3*lambda + 2*mu)/(lambda + mu)}; constexpr static Real poisson{lambda/(2*(lambda + mu))}; MaterialFixture():mat("Name", young, poisson){}; constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; template struct has_internals {constexpr static bool value{false};}; template struct has_internals> { constexpr static bool value{true}; }; using mat_list = boost::mpl::list >, MaterialFixture>, MaterialFixture>, MaterialFixture>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_add_pixel, Fix, mat_list, Fix) { auto & mat{Fix::mat}; constexpr Dim_t sdim{Fix::sdim}; testGoodies::RandRange rng;; const Dim_t nb_pixel{7}, box_size{17}; using Ccoord = Ccoord_t; for (Dim_t i = 0; i < nb_pixel; ++i) { Ccoord c; for (Dim_t j = 0; j < sdim; ++j) { c[j] = rng.randval(0, box_size); } if (!has_internals::value) { BOOST_CHECK_NO_THROW(mat.add_pixel(c)); } } BOOST_CHECK_NO_THROW(mat.initialise()); } template struct MaterialFixtureFilled: public MaterialFixture { using Mat = Mat_t; constexpr static Dim_t box_size{3}; MaterialFixtureFilled():MaterialFixture(){ using Ccoord = Ccoord_t; Ccoord cube{CcoordOps::get_cube(box_size)}; CcoordOps::Pixels pixels(cube); for (auto pixel: pixels) { this->mat.add_pixel(pixel); } this->mat.initialise(); }; }; using mat_fill = boost::mpl::list >, MaterialFixtureFilled>, MaterialFixtureFilled>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_law, Fix, mat_fill, Fix) { constexpr auto cube{CcoordOps::get_cube(Fix::box_size)}; auto & mat{Fix::mat}; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t globalfields; auto & F{make_field ("Transformation Gradient", globalfields)}; auto & P1 = make_field ("Nominal Stress1", globalfields); // to be computed alone auto & P2 = make_field ("Nominal Stress2", globalfields); // to be computed with tangent auto & K = make_field ("Tangent Moduli", globalfields); // to be computed with tangent auto & Pr = make_field ("Nominal Stress reference", globalfields); auto & Kr = make_field ("Tangent Moduli reference", globalfields); // to be computed with tangent globalfields.initialise(cube); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); using traits = MaterialMuSpectre_traits; { // block to contain not-constant gradient map typename traits::StressMap_t grad_map (globalfields["Transformation Gradient"]); for (auto F_: grad_map) { F_.setRandom(); } grad_map[0] = grad_map[0].Identity(); // identifiable gradients for debug grad_map[1] = 1.2*grad_map[1].Identity(); // ditto } //compute stresses using material mat.compute_stresses(globalfields["Transformation Gradient"], globalfields["Nominal Stress1"], Formulation::finite_strain); //compute stresses and tangent moduli using material BOOST_CHECK_THROW (mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Nominal Stress2"], Formulation::finite_strain), std::runtime_error); mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Tangent Moduli"], Formulation::finite_strain); typename traits::StrainMap_t Fmap(globalfields["Transformation Gradient"]); typename traits::StressMap_t Pmap_ref(globalfields["Nominal Stress reference"]); typename traits::TangentMap_t Kmap_ref(globalfields["Tangent Moduli reference"]); for (auto tup: akantu::zip(Fmap, Pmap_ref, Kmap_ref)) { auto F_ = std::get<0>(tup); auto P_ = std::get<1>(tup); auto K_ = std::get<2>(tup); std::tie(P_,K_) = testGoodies::objective_hooke_explicit (Fix::lambda, Fix::mu, F_); } typename traits::StressMap_t Pmap_1(globalfields["Nominal Stress1"]); for (auto tup: akantu::zip(Pmap_ref, Pmap_1)) { auto P_r = std::get<0>(tup); auto P_1 = std::get<1>(tup); Real error = (P_r - P_1).norm(); BOOST_CHECK_LT(error, tol); } typename traits::StressMap_t Pmap_2(globalfields["Nominal Stress2"]); typename traits::TangentMap_t Kmap(globalfields["Tangent Moduli"]); for (auto tup: akantu::zip(Pmap_ref, Pmap_2, Kmap_ref, Kmap)) { auto P_r = std::get<0>(tup); auto P = std::get<1>(tup); Real error = (P_r - P).norm(); BOOST_CHECK_LT(error, tol); auto K_r = std::get<2>(tup); auto K = std::get<3>(tup); error = (K_r - K).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_material_linear_elastic2.cc b/tests/test_material_linear_elastic2.cc index 0d3f4af..204df83 100644 --- a/tests/test_material_linear_elastic2.cc +++ b/tests/test_material_linear_elastic2.cc @@ -1,273 +1,273 @@ /** * @file test_material_linear_elastic2.cc * * @author Till Junge * * @date 04 Feb 2018 * * @brief Tests for the objective Hooke's law with eigenstrains, * (tests that do not require add_pixel are integrated into * `test_material_linear_elastic1.cc` * * @section LICENSE * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include "materials/material_linear_elastic2.hh" #include "tests.hh" #include "tests/test_goodies.hh" #include "common/field_collection.hh" #include "common/iterators.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(material_linear_elastic_2); template struct MaterialFixture { using Mat = Mat_t; constexpr static Real lambda{2}, mu{1.5}; constexpr static Real young{mu*(3*lambda + 2*mu)/(lambda + mu)}; constexpr static Real poisson{lambda/(2*(lambda + mu))}; MaterialFixture():mat("Name", young, poisson){}; constexpr static Dim_t sdim{Mat_t::sdim()}; constexpr static Dim_t mdim{Mat_t::mdim()}; Mat_t mat; }; using mat_list = boost::mpl::list< MaterialFixture>, MaterialFixture>, MaterialFixture>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_constructor, Fix, mat_list, Fix) { BOOST_CHECK_EQUAL("Name", Fix::mat.get_name()); auto & mat{Fix::mat}; auto sdim{Fix::sdim}; auto mdim{Fix::mdim}; BOOST_CHECK_EQUAL(sdim, mat.sdim()); BOOST_CHECK_EQUAL(mdim, mat.mdim()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_add_pixel, Fix, mat_list, Fix) { auto & mat{Fix::mat}; constexpr Dim_t sdim{Fix::sdim}; testGoodies::RandRange rng;; const Dim_t nb_pixel{7}, box_size{17}; using Ccoord = Ccoord_t; for (Dim_t i = 0; i < nb_pixel; ++i) { Ccoord c; for (Dim_t j = 0; j < sdim; ++j) { c[j] = rng.randval(0, box_size); } Eigen::Matrix Zero = Eigen::Matrix::Zero(); BOOST_CHECK_NO_THROW (mat.add_pixel(c, Zero)); } BOOST_CHECK_NO_THROW(mat.initialise()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_eigenstrain_equivalence, Fix, mat_list, Fix) { auto & mat{Fix::mat}; const Dim_t nb_pixel{2}; constexpr auto cube{CcoordOps::get_cube(nb_pixel)}; using Mat_t = Eigen::Matrix; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t globalfields; auto & F_f{make_field ("Transformation Gradient", globalfields)}; auto & P1_f = make_field ("Nominal Stress1", globalfields); // to be computed alone auto & K_f = make_field ("Tangent Moduli", globalfields); // to be computed with tangent globalfields.initialise(cube); Mat_t zero{Mat_t::Zero()}; Mat_t F{Mat_t::Random()/100 + Mat_t::Identity()}; Mat_t strain{-.5*(F+F.transpose())-Mat_t::Identity()}; using Ccoord = Ccoord_t; Ccoord pix0{0}; Ccoord pix1{1}; mat.add_pixel(pix0, zero); mat.add_pixel(pix1, strain); mat.initialise(); F_f.get_map()[pix0] = -strain; F_f.get_map()[pix1] = zero; mat.compute_stresses_tangent(F_f, P1_f, K_f, Formulation::small_strain); Real error{(P1_f.get_map()[pix0]-P1_f.get_map()[pix1]).norm()}; Real tol{1e-12}; if (error >= tol) { std::cout << "error = " << error << " >= " << tol << " = tol" << std::endl; std::cout << "P(0) =" << std::endl << P1_f.get_map()[pix0] << std::endl; std::cout << "P(1) =" << std::endl << P1_f.get_map()[pix1] << std::endl; } BOOST_CHECK_LT(error, tol); } template struct MaterialFixtureFilled: public MaterialFixture { using Par = MaterialFixture; using Mat = Mat_t; constexpr static Dim_t box_size{3}; MaterialFixtureFilled():MaterialFixture(){ using Ccoord = Ccoord_t; Ccoord cube{CcoordOps::get_cube(box_size)}; CcoordOps::Pixels pixels(cube); for (auto pixel: pixels) { Eigen::Matrix Zero = Eigen::Matrix::Zero(); this->mat.add_pixel(pixel, Zero); } this->mat.initialise(); }; }; using mat_fill = boost::mpl::list >, MaterialFixtureFilled>, MaterialFixtureFilled>>; BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_evaluate_law, Fix, mat_fill, Fix) { constexpr auto cube{CcoordOps::get_cube(Fix::box_size)}; auto & mat{Fix::mat}; - using FC_t = FieldCollection; + using FC_t = GlobalFieldCollection; FC_t globalfields; auto & F{make_field ("Transformation Gradient", globalfields)}; auto & P1 = make_field ("Nominal Stress1", globalfields); // to be computed alone auto & P2 = make_field ("Nominal Stress2", globalfields); // to be computed with tangent auto & K = make_field ("Tangent Moduli", globalfields); // to be computed with tangent auto & Pr = make_field ("Nominal Stress reference", globalfields); auto & Kr = make_field ("Tangent Moduli reference", globalfields); // to be computed with tangent globalfields.initialise(cube); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); static_assert(std::is_same::value, "oh oh"); using traits = MaterialMuSpectre_traits; { // block to contain not-constant gradient map typename traits::StressMap_t grad_map (globalfields["Transformation Gradient"]); for (auto F_: grad_map) { F_.setRandom(); } grad_map[0] = grad_map[0].Identity(); // identifiable gradients for debug grad_map[1] = 1.2*grad_map[1].Identity(); // ditto } //compute stresses using material mat.compute_stresses(globalfields["Transformation Gradient"], globalfields["Nominal Stress1"], Formulation::finite_strain); //compute stresses and tangent moduli using material BOOST_CHECK_THROW (mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Nominal Stress2"], Formulation::finite_strain), std::runtime_error); mat.compute_stresses_tangent(globalfields["Transformation Gradient"], globalfields["Nominal Stress2"], globalfields["Tangent Moduli"], Formulation::finite_strain); typename traits::StrainMap_t Fmap(globalfields["Transformation Gradient"]); typename traits::StressMap_t Pmap_ref(globalfields["Nominal Stress reference"]); typename traits::TangentMap_t Kmap_ref(globalfields["Tangent Moduli reference"]); for (auto tup: akantu::zip(Fmap, Pmap_ref, Kmap_ref)) { auto F_ = std::get<0>(tup); auto P_ = std::get<1>(tup); auto K_ = std::get<2>(tup); std::tie(P_,K_) = testGoodies::objective_hooke_explicit (Fix::lambda, Fix::mu, F_); } typename traits::StressMap_t Pmap_1(globalfields["Nominal Stress1"]); for (auto tup: akantu::zip(Pmap_ref, Pmap_1)) { auto P_r = std::get<0>(tup); auto P_1 = std::get<1>(tup); Real error = (P_r - P_1).norm(); BOOST_CHECK_LT(error, tol); } typename traits::StressMap_t Pmap_2(globalfields["Nominal Stress2"]); typename traits::TangentMap_t Kmap(globalfields["Tangent Moduli"]); for (auto tup: akantu::zip(Pmap_ref, Pmap_2, Kmap_ref, Kmap)) { auto P_r = std::get<0>(tup); auto P = std::get<1>(tup); Real error = (P_r - P).norm(); if (error >= tol) { std::cout << "error = " << error << " >= " << tol << " = tol" << std::endl; std::cout << "P(0) =" << std::endl << P_r << std::endl; std::cout << "P(1) =" << std::endl << P << std::endl; } BOOST_CHECK_LT(error, tol); auto K_r = std::get<2>(tup); auto K = std::get<3>(tup); error = (K_r - K).norm(); BOOST_CHECK_LT(error, tol); } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 85d4ada..7c6e83c 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,127 +1,127 @@ /** * @file test_projection_finite.cc * * @author Till Junge * * @date 07 Dec 2017 * * @brief tests for standard finite strain 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. */ #include "fft/projection_finite_strain.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/fft_utils.hh" #include "test_projection.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_finite_strain); /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list< ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrain>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>, ProjectionFixture, ProjectionFiniteStrainFast>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { // create a gradient field with a zero mean gradient and verify // that the projection preserves it constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim}; static_assert(dim == fix::mdim, "These tests assume that the material and spatial dimension are " "identical"); - using Fields = FieldCollection; + using Fields = GlobalFieldCollection; using FieldT = TensorField; using FieldMap = MatrixFieldMap; using Vector = Eigen::Matrix; Fields fields{}; FieldT & f_grad{make_field("gradient", fields)}; FieldT & f_var{make_field("working field", fields)}; FieldMap grad(f_grad); FieldMap var(f_var); fields.initialise(fix::projector.get_resolutions()); FFT_freqs freqs{fix::projector.get_resolutions(), fix::projector.get_lengths()}; Vector k; for (Dim_t i = 0; i < dim; ++i) { // the wave vector has to be such that it leads to an integer // number of periods in each length of the domain k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i]; ; } for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_resolutions()); g.row(0) = k.transpose() * cos(k.dot(vec)); v.row(0) = g.row(0); } fix::projector.initialise(FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_resolutions()); Real error = (g-v).norm(); BOOST_CHECK_LT(error, tol); if (error >=tol) { std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl; std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl; std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl; std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl; } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_projection_small.cc b/tests/test_projection_small.cc index 8b0d3ba..ae568b8 100644 --- a/tests/test_projection_small.cc +++ b/tests/test_projection_small.cc @@ -1,129 +1,129 @@ /** * @file test_projection_small.cc * * @author Till Junge * * @date 16 Jan 2018 * * @brief tests for standard small strain projection operator * * 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 "fft/projection_small_strain.hh" #include "test_projection.hh" #include "fft/fft_utils.hh" #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_small_strain); using fixlist = boost::mpl::list< ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>, ProjectionFixture, ProjectionSmallStrain>>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { // create a gradient field with a zero mean gradient and verify // that the projection preserves it constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim}; static_assert(dim == fix::mdim, "These tests assume that the material and spatial dimension are " "identical"); - using Fields = FieldCollection; + using Fields = GlobalFieldCollection; using FieldT = TensorField; using FieldMap = MatrixFieldMap; using Vector = Eigen::Matrix; Fields fields{}; FieldT & f_grad{make_field("strain", fields)}; FieldT & f_var{make_field("working field", fields)}; FieldMap grad(f_grad); FieldMap var(f_var); fields.initialise(fix::projector.get_resolutions()); FFT_freqs freqs{fix::projector.get_resolutions(), fix::projector.get_lengths()}; Vector k; for (Dim_t i = 0; i < dim; ++i) { // the wave vector has to be such that it leads to an integer // number of periods in each length of the domain k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i]; } for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_resolutions()); g.row(0) << k.transpose() * cos(k.dot(vec)); // We need to add I to the term, because this field has a net // zero gradient, which leads to a net -I strain g = 0.5*((g-g.Identity()).transpose() + (g-g.Identity())).eval()+g.Identity(); v = g; } fix::projector.initialise(FFT_PlanFlags::estimate); fix::projector.apply_projection(f_var); constexpr bool verbose{false}; for (auto && tup: akantu::zip(fields, grad, var)) { auto & ccoord = std::get<0>(tup); auto & g = std::get<1>(tup); auto & v = std::get<2>(tup); Vector vec = CcoordOps::get_vector(ccoord, fix::projector.get_lengths()/ fix::projector.get_resolutions()); Real error = (g-v).norm(); BOOST_CHECK_LT(error, tol); if ((error >=tol) || verbose) { std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl; std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl; std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl; std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl; std::cout << "means:" << std::endl << ":" << std::endl << grad.mean() << std::endl << ":" << std::endl << var.mean(); } } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre