diff --git a/src/cell/cell_base.cc b/src/cell/cell_base.cc index 2af17f1..502ae4c 100644 --- a/src/cell/cell_base.cc +++ b/src/cell/cell_base.cc @@ -1,295 +1,295 @@ /** * @file cell_base.cc * * @author Till Junge * * @date 01 Nov 2017 * * @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 "cell/cell_base.hh" #include "common/ccoord_operations.hh" #include "common/iterators.hh" #include "common/tensor_algebra.hh" #include #include namespace muSpectre { /* ---------------------------------------------------------------------- */ template CellBase::CellBase(Projection_ptr projection_) :subdomain_resolutions{projection_->get_subdomain_resolutions()}, subdomain_locations{projection_->get_subdomain_locations()}, domain_resolutions{projection_->get_domain_resolutions()}, pixels(subdomain_resolutions, subdomain_locations), - lengths{projection_->get_lengths()}, + domain_lengths{projection_->get_domain_lengths()}, fields{std::make_unique()}, F{make_field("Gradient", *this->fields)}, P{make_field("Piola-Kirchhoff-1", *this->fields)}, projection{std::move(projection_)}, form{projection->get_formulation()} { } /* ---------------------------------------------------------------------- */ template typename CellBase::Material_t & CellBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template typename CellBase::FullResponse_t CellBase::evaluate_stress_tangent(StrainField_t & grad) { if (this->initialised == false) { this->initialise(); } //! High level compatibility checks if (grad.size() != this->F.size()) { throw std::runtime_error("Size mismatch"); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(grad, this->P, this->K.value(), this->form); } return std::tie(this->P, this->K.value()); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::directional_stiffness(const TangentField_t &K, const StrainField_t &delF, StressField_t &delP) { for (auto && tup: akantu::zip(K.get_map(), delF.get_map(), delP.get_map())){ auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return this->project(delP); } /* ---------------------------------------------------------------------- */ template typename CellBase::SolvVectorOut CellBase::directional_stiffness_vec(const SolvVectorIn &delF) { if (!this->K) { throw std::runtime_error ("corrently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); SolvVectorOut(in_tempref.data(), this->nb_dof()) = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return SolvVectorOut(out_tempref.data(), this->nb_dof()); } /* ---------------------------------------------------------------------- */ template Eigen::ArrayXXd CellBase:: directional_stiffness_with_copy (Eigen::Ref delF) { if (!this->K) { throw std::runtime_error ("corrently only implemented for cases where a stiffness matrix " "exists"); } const std::string out_name{"temp output for directional stiffness"}; const std::string in_name{"temp input for directional stiffness"}; auto & out_tempref = this->get_managed_field(out_name); auto & in_tempref = this->get_managed_field(in_name); in_tempref.eigen() = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return out_tempref.eigen(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StressField_t & CellBase::project(StressField_t &field) { this->projection->apply_projection(field); return field; } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_strain() { if (this->initialised == false) { this->initialise(); } return this->F; } /* ---------------------------------------------------------------------- */ template const typename CellBase::StressField_t & CellBase::get_stress() const { return this->P; } /* ---------------------------------------------------------------------- */ template const typename CellBase::TangentField_t & CellBase::get_tangent(bool create) { if (!this->K) { if (create) { this->K = make_field("Tangent Stiffness", *this->fields); } else { throw std::runtime_error ("K does not exist"); } } return this->K.value(); } /* ---------------------------------------------------------------------- */ template typename CellBase::StrainField_t & CellBase::get_managed_field(std::string unique_name) { if (!this->fields->check_field_exists(unique_name)) { return make_field(unique_name, *this->fields); } else { return static_cast(this->fields->at(unique_name)); } } /* ---------------------------------------------------------------------- */ template void CellBase::initialise(FFT_PlanFlags flags) { // check that all pixels have been assigned exactly one material this->check_material_coverage(); // resize all global fields (strain, stress, etc) this->fields->initialise(this->subdomain_resolutions, this->subdomain_locations); // initialise the projection and compute the fft plan this->projection->initialise(flags); this->initialised = true; } /* ---------------------------------------------------------------------- */ template void CellBase::initialise_materials(bool stiffness) { for (auto && mat: this->materials) { mat->initialise(stiffness); } } /* ---------------------------------------------------------------------- */ template void CellBase::save_history_variables() { for (auto && mat: this->materials) { mat->save_history_variables(); } } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::begin() { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename CellBase::iterator CellBase::end() { return this->pixels.end(); } /* ---------------------------------------------------------------------- */ template CellAdaptor> CellBase::get_adaptor() { return Adaptor(*this); } /* ---------------------------------------------------------------------- */ template void CellBase::check_material_coverage() { auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions); std::vector*> assignments(nb_pixels, nullptr); for (auto & mat: this->materials) { for (auto & pixel: *mat) { auto index = CcoordOps::get_index(this->subdomain_resolutions, this->subdomain_locations, pixel); auto& assignment{assignments.at(index)}; if (assignment != nullptr) { std::stringstream err{}; err << "Pixel " << pixel << "is already assigned to material '" << assignment->get_name() << "' and cannot be reassigned to material '" << mat->get_name(); throw std::runtime_error(err.str()); } else { assignments[index] = mat.get(); } } } // find and identify unassigned pixels std::vector unassigned_pixels; for (size_t i = 0; i < assignments.size(); ++i) { if (assignments[i] == nullptr) { unassigned_pixels.push_back( CcoordOps::get_ccoord(this->subdomain_resolutions, this->subdomain_locations, i)); } } if (unassigned_pixels.size() != 0) { std::stringstream err {}; err << "The following pixels have were not assigned a material: "; for (auto & pixel: unassigned_pixels) { err << pixel << ", "; } err << "and that cannot be handled"; throw std::runtime_error(err.str()); } } template class CellBase; template class CellBase; } // muSpectre diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh index 8cc9a2a..018ae27 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,313 +1,313 @@ /** * @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 namespace muSpectre { template class CellAdaptor; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase { public: 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; //! input vector for solvers using SolvVectorIn = Eigen::Ref; //! output vector for solvers using SolvVectorOut = Eigen::Map; //! sparse matrix emulation using Adaptor = CellAdaptor; //! Default constructor CellBase() = delete; //! constructor using sizes and resolution CellBase(Projection_ptr projection); //! Copy constructor CellBase(const CellBase &other) = delete; //! Move constructor CellBase(CellBase &&other) = default; //! 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); /** * 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 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); /** * initialise materials (including resetting any history variables) */ void initialise_materials(bool stiffness=false); /** * 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(); 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_lengths() const {return this->lengths;} + const Rcoord & get_domain_lengths() const {return this->domain_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);}; //! return the communicator object const Communicator & get_communicator() const { return this->projection->get_communicator(); } protected: //! make sure that every pixel is assigned to one and only one material 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 & lengths; //!< the cell's lengths + 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 const Formulation form; //!< formulation for solution private: }; /** * lightweight resource handle wrapping a `muSpectre::CellBase` 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.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.directional_stiffness_vec(rhs); } }; } } #endif /* CELL_BASE_H */ diff --git a/src/fft/projection_base.cc b/src/fft/projection_base.cc index 351e62c..c291fc6 100644 --- a/src/fft/projection_base.cc +++ b/src/fft/projection_base.cc @@ -1,57 +1,57 @@ /** * @file projection_base.cc * * @author Till Junge * * @date 06 Dec 2017 * * @brief implementation of base class for projections * * 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_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionBase::ProjectionBase(FFTEngine_ptr engine, - Rcoord lengths, + Rcoord domain_lengths, Formulation form) - : fft_engine{std::move(engine)}, lengths{lengths}, form{form}, + : fft_engine{std::move(engine)}, domain_lengths{domain_lengths}, form{form}, projection_container{this->fft_engine->get_field_collection()} { static_assert((DimS == FFTEngine::sdim), "spatial dimensions are incompatible"); static_assert((DimM == FFTEngine::mdim), "material dimensions are incompatible"); } /* ---------------------------------------------------------------------- */ template void ProjectionBase:: initialise(FFT_PlanFlags flags) { fft_engine->initialise(flags); } template class ProjectionBase; template class ProjectionBase; template class ProjectionBase; } // muSpectre diff --git a/src/fft/projection_base.hh b/src/fft/projection_base.hh index 0b27e7c..6ff89b5 100644 --- a/src/fft/projection_base.hh +++ b/src/fft/projection_base.hh @@ -1,169 +1,169 @@ /** * @file projection_base.hh * * @author Till Junge * * @date 03 Dec 2017 * * @brief Base class for Projection operators * * 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_BASE_H #define PROJECTION_BASE_H #include "common/common.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "fft/fft_engine_base.hh" #include namespace muSpectre { /** * base class for projection related exceptions */ class ProjectionError: public std::runtime_error { public: //! constructor explicit ProjectionError(const std::string& what) :std::runtime_error(what){} //! constructor explicit ProjectionError(const char * what) :std::runtime_error(what){} }; template struct Projection_traits { }; /** * defines the interface which must be implemented by projection operators */ template class ProjectionBase { public: //! type of fft_engine used using FFTEngine = FFTEngineBase; //! reference to fft engine is safely managed through a `std::unique_ptr` using FFTEngine_ptr = std::unique_ptr; //! cell coordinates type using Ccoord = typename FFTEngine::Ccoord; //! spatial coordinates type using Rcoord = std::array; //! global FieldCollection using GFieldCollection_t = typename FFTEngine::GFieldCollection_t; //! local FieldCollection (for Fourier-space pixels) using LFieldCollection_t = typename FFTEngine::LFieldCollection_t; //! Field type on which to apply the projection using Field_t = typename FFTEngine::Field_t; /** * iterator over all pixels. This is taken from the FFT engine, * because depending on the real-to-complex FFT employed, only * roughly half of the pixels are present in Fourier space * (because of the hermitian nature of the transform) */ using iterator = typename FFTEngine::iterator; //! Default constructor ProjectionBase() = delete; //! Constructor with cell sizes - ProjectionBase(FFTEngine_ptr engine, Rcoord lengths, Formulation form); + ProjectionBase(FFTEngine_ptr engine, Rcoord domain_lengths, Formulation form); //! Copy constructor ProjectionBase(const ProjectionBase &other) = delete; //! Move constructor ProjectionBase(ProjectionBase &&other) = default; //! Destructor virtual ~ProjectionBase() = default; //! Copy assignment operator ProjectionBase& operator=(const ProjectionBase &other) = delete; //! Move assignment operator ProjectionBase& operator=(ProjectionBase &&other) = default; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); //! apply the projection operator to a field virtual void apply_projection(Field_t & field) = 0; //! returns the process-local resolutions of the cell const Ccoord & get_subdomain_resolutions() const { return this->fft_engine->get_subdomain_resolutions();} //! returns the process-local locations of the cell const Ccoord & get_subdomain_locations() const { return this->fft_engine->get_subdomain_locations();} //! returns the resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->fft_engine->get_domain_resolutions();} //! returns the physical sizes of the cell - const Rcoord & get_lengths() const {return this->lengths;} + const Rcoord & get_domain_lengths() const {return this->domain_lengths;} /** * return the `muSpectre::Formulation` that is used in solving * this cell. This allows tho check whether a projection is * compatible with the chosen formulation */ const Formulation & get_formulation() const {return this->form;} //! return the raw projection operator. This is mainly intended //! for maintenance and debugging and should never be required in //! regular use virtual Eigen::Map get_operator() = 0; //! return the communicator object const Communicator & get_communicator() const { return this->fft_engine->get_communicator(); } protected: //! handle on the fft_engine used FFTEngine_ptr fft_engine; - const Rcoord lengths; //!< physical sizes of the cell + const Rcoord domain_lengths; //!< physical sizes of the cell /** * formulation this projection can be applied to (determines * whether the projection enforces gradients, small strain tensor * or symmetric smal strain tensor */ const Formulation form; /** * A local `muSpectre::FieldCollection` to store the projection * operator per k-space point. This is a local rather than a * global collection, since the pixels considered depend on the * FFT implementation. See * http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data * for an example */ LFieldCollection_t & projection_container{}; private: }; } // muSpectre #endif /* PROJECTION_BASE_H */ diff --git a/src/fft/projection_finite_strain.cc b/src/fft/projection_finite_strain.cc index 36f66d3..55695ef 100644 --- a/src/fft/projection_finite_strain.cc +++ b/src/fft/projection_finite_strain.cc @@ -1,92 +1,92 @@ /** * @file projection_finite_strain.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief implementation of 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/fftw_engine.hh" #include "fft/fft_utils.hh" #include "common/field_map.hh" #include "common/tensor_algebra.hh" #include "common/iterators.hh" #include "Eigen/Dense" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrain:: ProjectionFiniteStrain(FFTEngine_ptr engine, Rcoord lengths) :Parent{std::move(engine), lengths, Formulation::finite_strain} { for (auto res: this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError ("Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrain:: initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->lengths); + this->domain_lengths); for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0> (tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); //! this is simplifiable using Curnier's Méthodes numériques, 6.69(c) G = Matrices::outer_under(Matrices::I2(), xi*xi.transpose()); // The commented block below corresponds to the original // definition of the operator in de Geus et // al. (https://doi.org/10.1016/j.cma.2016.12.032). However, // they use a bizarre definition of the double contraction // between fourth-order and second-order tensors that has a // built-in transpose operation (i.e., C = A:B <-> AᵢⱼₖₗBₗₖ = // Cᵢⱼ , note the inverted ₗₖ instead of ₖₗ), here, we define // the double contraction without the transposition. As a // result, the Projection operator produces the transpose of de // Geus's // for (Dim_t im = 0; im < DimS; ++im) { // for (Dim_t j = 0; j < DimS; ++j) { // for (Dim_t l = 0; l < DimS; ++l) { // get(G, im, j, l, im) = xi(j)*xi(l); // } // } // } } if (this->get_subdomain_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionFiniteStrain; template class ProjectionFiniteStrain; } // muSpectre diff --git a/src/fft/projection_finite_strain_fast.cc b/src/fft/projection_finite_strain_fast.cc index 6b20e08..dd3300c 100644 --- a/src/fft/projection_finite_strain_fast.cc +++ b/src/fft/projection_finite_strain_fast.cc @@ -1,92 +1,92 @@ /** * @file projection_finite_strain_fast.cc * * @author Till Junge * * @date 12 Dec 2017 * * @brief implementation for fast projection in finite strain * * 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_fast.hh" #include "fft/fft_utils.hh" #include "common/tensor_algebra.hh" #include "common/iterators.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrainFast:: ProjectionFiniteStrainFast(FFTEngine_ptr engine, Rcoord lengths) :Parent{std::move(engine), lengths, Formulation::finite_strain}, xiField{make_field("Projection Operator", this->projection_container)}, xis(xiField) { for (auto res: this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError ("Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrainFast:: initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->lengths); + this->domain_lengths); for (auto && tup: akantu::zip(*this->fft_engine, this->xis)) { const auto & ccoord = std::get<0> (tup); auto & xi = std::get<1>(tup); xi = fft_freqs.get_unit_xi(ccoord); } if (this->get_subdomain_locations() == Ccoord{}) { this->xis[0].setZero(); } } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrainFast::apply_projection(Field_t & field) { Grad_map field_map{this->fft_engine->fft(field)}; Real factor = this->fft_engine->normalisation(); for (auto && tup: akantu::zip(this->xis, field_map)) { auto & xi{std::get<0>(tup)}; auto & f{std::get<1>(tup)}; f = factor * ((f*xi).eval()*xi.transpose()); } this->fft_engine->ifft(field); } /* ---------------------------------------------------------------------- */ template Eigen::Map ProjectionFiniteStrainFast:: get_operator() { return this->xiField.dyn_eigen(); } template class ProjectionFiniteStrainFast; template class ProjectionFiniteStrainFast; } // muSpectre diff --git a/src/fft/projection_small_strain.cc b/src/fft/projection_small_strain.cc index f3581f7..9d86625 100644 --- a/src/fft/projection_small_strain.cc +++ b/src/fft/projection_small_strain.cc @@ -1,83 +1,83 @@ /** * @file projection_small_strain.cc * * @author Till Junge * * @date 14 Jan 2018 * * @brief Implementation for ProjectionSmallStrain * * 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 "fft/fft_utils.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionSmallStrain:: ProjectionSmallStrain(FFTEngine_ptr engine, Rcoord lengths) : Parent{std::move(engine), lengths, Formulation::small_strain} { for (auto res: this->fft_engine->get_domain_resolutions()) { if (res % 2 == 0) { throw ProjectionError ("Only an odd number of gridpoints in each direction is supported"); } } } /* ---------------------------------------------------------------------- */ template void ProjectionSmallStrain::initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine->get_domain_resolutions(), - this->lengths); + this->domain_lengths); for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) { const auto & ccoord = std::get<0> (tup); auto & G = std::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); auto kron = [](const Dim_t i, const Dim_t j) -> Real{ return (i==j) ? 1. : 0.; }; for (Dim_t i{0}; i < DimS; ++i) { for (Dim_t j{0}; j < DimS; ++j) { for (Dim_t l{0}; l < DimS; ++l) { for (Dim_t m{0}; m < DimS; ++m ) { Real & g = get(G, i, j, l, m); g = 0.5* (xi(i) * kron(j, l) * xi(m) + xi(i) * kron(j, m) * xi(l) + xi(j) * kron(i, l) * xi(m) + xi(j) * kron(i, m) * xi(l)) - xi(i)*xi(j)*xi(l)*xi(m); } } } } } if (this->get_subdomain_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionSmallStrain; template class ProjectionSmallStrain; } // muSpectre diff --git a/tests/mpi_test_projection_finite.cc b/tests/mpi_test_projection_finite.cc index ac4e832..c9ad967 100644 --- a/tests/mpi_test_projection_finite.cc +++ b/tests/mpi_test_projection_finite.cc @@ -1,189 +1,189 @@ /** * @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. */ #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 50 #include "fft/projection_finite_strain.hh" #include "fft/projection_finite_strain_fast.hh" #include "fft/fft_utils.hh" #include "mpi_test_projection.hh" #include "fft/fftw_engine.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(mpi_projection_finite_strain); /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list< #ifdef WITH_FFTWMPI ProjectionFixture, ProjectionFiniteStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrainFast, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrainFast, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrainFast, FFTWMPIEngine>, ProjectionFixture, ProjectionFiniteStrainFast, FFTWMPIEngine>, #endif #ifdef WITH_PFFT ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrain, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, ProjectionFixture, ProjectionFiniteStrainFast, PFFTEngine>, #endif ProjectionFixture, ProjectionFiniteStrain, FFTWEngine, false> >; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { if (fix::is_parallel || fix::projector.get_communicator().size() == 1) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } else { BOOST_CHECK_THROW(fix::projector.initialise(FFT_PlanFlags::estimate), std::runtime_error); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) { return; } // 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 = 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_subdomain_resolutions(), fix::projector.get_subdomain_locations()); 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]; + k(i) = (i+1)*2*pi/fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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); Real error = (g-v).norm(); BOOST_CHECK_LT(error, tol); if (error >=tol) { Vector vec = CcoordOps::get_vector(ccoord, - fix::projector.get_lengths()/ + fix::projector.get_domain_lengths()/ fix::projector.get_domain_resolutions()); 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/mpi_test_projection_small.cc b/tests/mpi_test_projection_small.cc index f460802..526d99c 100644 --- a/tests/mpi_test_projection_small.cc +++ b/tests/mpi_test_projection_small.cc @@ -1,173 +1,173 @@ /** * @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. */ #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS #define BOOST_MPL_LIMIT_LIST_SIZE 50 #include "fft/projection_small_strain.hh" #include "mpi_test_projection.hh" #include "fft/fft_utils.hh" #include "fft/fftw_engine.hh" #ifdef WITH_FFTWMPI #include "fft/fftwmpi_engine.hh" #endif #ifdef WITH_PFFT #include "fft/pfft_engine.hh" #endif #include namespace muSpectre { BOOST_AUTO_TEST_SUITE(mpi_projection_small_strain); using fixlist = boost::mpl::list< #ifdef WITH_FFTWMPI ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, ProjectionFixture, ProjectionSmallStrain, FFTWMPIEngine>, #endif #ifdef WITH_PFFT ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, ProjectionFixture, ProjectionSmallStrain, PFFTEngine>, #endif ProjectionFixture, ProjectionSmallStrain, FFTWEngine, false> >; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { if (fix::is_parallel || fix::projector.get_communicator().size() == 1) { BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } else { BOOST_CHECK_THROW(fix::projector.initialise(FFT_PlanFlags::estimate), std::runtime_error); } } /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test, fix, fixlist, fix) { if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) { return; } // 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 = 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_subdomain_resolutions(), fix::projector.get_subdomain_locations()); 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]; + k(i) = (i+1)*2*pi/fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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 diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 1c21380..916f0b3 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,137 +1,137 @@ /** * @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_AUTO_TEST_CASE(even_grid_test) { using Engine = FFTWEngine; using proj = ProjectionFiniteStrainFast; auto engine = std::make_unique(Ccoord_t{2, 2}); BOOST_CHECK_THROW(proj(std::move(engine), Rcoord_t{4.3, 4.3}), std::runtime_error); } /* ---------------------------------------------------------------------- */ 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 = 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_subdomain_resolutions(), fix::projector.get_subdomain_locations()); FFT_freqs freqs{fix::projector.get_domain_resolutions(), - fix::projector.get_lengths()}; + fix::projector.get_domain_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]; ; + k(i) = (i+1)*2*pi/fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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 3ae0f1d..f8ae988 100644 --- a/tests/test_projection_small.cc +++ b/tests/test_projection_small.cc @@ -1,130 +1,130 @@ /** * @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 = 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_subdomain_resolutions(), fix::projector.get_subdomain_locations()); FFT_freqs freqs{fix::projector.get_domain_resolutions(), - fix::projector.get_lengths()}; + fix::projector.get_domain_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]; + k(i) = (i+1)*2*pi/fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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_domain_lengths()/ fix::projector.get_domain_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