diff --git a/src/common/ccoord_operations.hh b/src/common/ccoord_operations.hh index 11af0e2..b489c2b 100644 --- a/src/common/ccoord_operations.hh +++ b/src/common/ccoord_operations.hh @@ -1,170 +1,208 @@ /** * file ccoord_operations.hh * * @author Till Junge * * @date 29 Sep 2017 * * @brief common operations on pixel addressing * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include +#include + #include "common/common.hh" #ifndef CCOORD_OPERATIONS_H #define CCOORD_OPERATIONS_H namespace muSpectre { namespace CcoordOps { namespace internal { constexpr Dim_t ret(Dim_t val, size_t /*dummy*/) {return val;} template constexpr Ccoord_t funt(Dim_t val, std::index_sequence) { return Ccoord_t{ret(val, I)...}; } template constexpr Ccoord_t herm(const Ccoord_t & full_sizes, std::index_sequence) { return Ccoord_t{full_sizes[I]..., (full_sizes.back()+1)/2}; } } // internal //----------------------------------------------------------------------------// template constexpr Ccoord_t get_cube(Dim_t size) { return internal::funt(size, std::make_index_sequence{}); } /* ---------------------------------------------------------------------- */ template constexpr Ccoord_t get_hermitian_sizes(Ccoord_t full_sizes) { return internal::herm(full_sizes, std::make_index_sequence{}); } + /* ---------------------------------------------------------------------- */ + template + Eigen::Matrix get_vector(const Ccoord_t & ccoord, T pix_size = T{1.}) { + Eigen::Matrix retval; + for (size_t i = 0; i < dim; ++i) { + retval[i] = pix_size * ccoord[i]; + } + return retval; + } + + /* ---------------------------------------------------------------------- */ + template + Eigen::Matrix get_vector(const Ccoord_t & ccoord, + Eigen::Matrix pix_size) { + Eigen::Matrix retval = pix_size; + for (size_t i = 0; i < dim; ++i) { + retval[i] *= ccoord[i]; + } + return retval; + } + + //----------------------------------------------------------------------------// template constexpr Ccoord_t get_ccoord(const Ccoord_t & sizes, Dim_t index) { Ccoord_t retval{0}; Dim_t factor{1}; for (Dim_t i = dim-1; i >=0; --i) { retval[i] = index/factor%sizes[i]; if (i != 0 ) { factor *= sizes[i]; } } return retval; } //----------------------------------------------------------------------------// template constexpr Dim_t get_index(const Ccoord_t & sizes, const Ccoord_t & ccoord) { Dim_t retval{0}; Dim_t factor{1}; for (Dim_t i = dim-1; i >=0; --i) { retval += ccoord[i]*factor; if (i != 0) { factor *= sizes[i]; } } return retval; } //----------------------------------------------------------------------------// template constexpr size_t get_size(const Ccoord_t& sizes) { Dim_t retval{1}; for (size_t i = 0; i < dim; ++i) { retval *= sizes[i]; } return retval; } /* ---------------------------------------------------------------------- */ template class Pixels { public: Pixels(const Ccoord_t & sizes):sizes(sizes){}; + Pixels(const Pixels & other) = default; + Pixels & operator=(const Pixels & other) = default; virtual ~Pixels() = default; class iterator { public: using value_type = Ccoord_t; + using const_value_type = const value_type; + using pointer = value_type*; + using difference_type = std::ptrdiff_t; using iterator_category = std::forward_iterator_tag; + using reference = value_type; iterator(const Pixels & pixels, bool begin=true); virtual ~iterator() = default; inline value_type operator*() const; inline iterator & operator++(); inline bool operator!=(const iterator & other) const; + inline bool operator==(const iterator & other) const; protected: const Pixels& pixels; size_t index; }; inline iterator begin() const {return iterator(*this);} inline iterator end() const {return iterator(*this, false);} protected: - const Ccoord_t & sizes; + Ccoord_t sizes; }; /* ---------------------------------------------------------------------- */ template Pixels::iterator::iterator(const Pixels & pixels, bool begin) :pixels{pixels}, index{begin? 0: get_size(pixels.sizes)} {} /* ---------------------------------------------------------------------- */ template typename Pixels::iterator::value_type Pixels::iterator::operator*() const { return get_ccoord(pixels.sizes, this->index); } /* ---------------------------------------------------------------------- */ template bool Pixels::iterator::operator!=(const iterator &other) const { return (this->index != other.index) || (&this->pixels != &other.pixels); } + /* ---------------------------------------------------------------------- */ + template + bool + Pixels::iterator::operator==(const iterator &other) const { + return !(*this!= other); + } + /* ---------------------------------------------------------------------- */ template typename Pixels::iterator& Pixels::iterator::operator++() { ++this->index; return *this; } } // CcoordOps } // muSpectre #endif /* CCOORD_OPERATIONS_H */ diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh index aa158d2..d04e0c5 100644 --- a/src/common/field_collection_global.hh +++ b/src/common/field_collection_global.hh @@ -1,165 +1,186 @@ /** * file field_collection_global.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for global fields * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef FIELD_COLLECTION_GLOBAL_H #define FIELD_COLLECTION_GLOBAL_H #include "common/field_collection_base.hh" #include "common/ccoord_operations.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ //! DimS spatial dimension (dimension of problem) //! DimM material_dimension (dimension of constitutive law) template class GlobalFieldCollection: public FieldCollectionBase> { public: using Parent = FieldCollectionBase >; using Ccoord = typename Parent::Ccoord; using Field_p = typename Parent::Field_p; + using iterator = typename CcoordOps::Pixels::iterator; //! Default constructor GlobalFieldCollection(); //! Copy constructor GlobalFieldCollection(const GlobalFieldCollection &other) = delete; //! Move constructor GlobalFieldCollection (GlobalFieldCollection &&other) noexcept = delete; //! Destructor virtual ~GlobalFieldCollection() noexcept = default; //! Copy assignment operator GlobalFieldCollection& operator=(const GlobalFieldCollection &other) = delete; //! Move assignment operator GlobalFieldCollection& operator=(GlobalFieldCollection &&other) noexcept = delete; /** allocate memory, etc. At this point, the collection is informed aboud the size and shape of the domain (through the sizes parameter). The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' any field of a different size is wrong TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(Ccoord sizes); //! return the pixel sizes inline const Ccoord & get_sizes() const; //! returns the linear index corresponding to cell coordinates template inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; + inline iterator begin(); + inline iterator end(); + + static constexpr inline Dim_t spatial_dim() {return DimS;} static constexpr inline Dim_t material_dim() {return DimM;} protected: //! number of discretisation cells in each of the DimS spatial directions Ccoord sizes{0}; + CcoordOps::Pixels pixels; private: }; /* ---------------------------------------------------------------------- */ template GlobalFieldCollection::GlobalFieldCollection() - :Parent(){} + :Parent(), pixels{{0}}{} /* ---------------------------------------------------------------------- */ template void GlobalFieldCollection:: initialise(Ccoord sizes) { + this->pixels = CcoordOps::Pixels(sizes); this->size_ = std::accumulate(sizes.begin(), sizes.end(), 1, std::multiplies()); this->sizes = sizes; std::for_each(std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! return the pixel sizes template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_sizes() const { return this->sizes; } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename GlobalFieldCollection::Ccoord GlobalFieldCollection::get_ccoord(size_t index) const { return CcoordOps::get_ccoord(this->get_sizes(), std::move(index)); } + + /* ---------------------------------------------------------------------- */ + template + typename GlobalFieldCollection::iterator + GlobalFieldCollection::begin() { + return this->pixels.begin(); + } + + /* ---------------------------------------------------------------------- */ + template + typename GlobalFieldCollection::iterator + GlobalFieldCollection::end() { + return this->pixels.end(); + } //----------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template template size_t GlobalFieldCollection::get_index(CcoordRef && ccoord) const { static_assert(std::is_same< Ccoord, std::remove_const_t< std::remove_reference_t>>::value, "can only be called with values or references of Ccoord"); return CcoordOps::get_index(this->get_sizes(), std::forward(ccoord)); } } // muSpectre #endif /* FIELD_COLLECTION_GLOBAL_H */ diff --git a/src/fft/projection_finite_strain.cc b/src/fft/projection_finite_strain.cc index b6ba58d..e8212bb 100644 --- a/src/fft/projection_finite_strain.cc +++ b/src/fft/projection_finite_strain.cc @@ -1,87 +1,88 @@ /** * file projection_finite_strain.cc * * @author Till Junge * * @date 05 Dec 2017 * * @brief implementation of standard finite strain projection operator * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include "Eigen/Dense" #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" namespace muSpectre { /* ---------------------------------------------------------------------- */ template ProjectionFiniteStrain:: ProjectionFiniteStrain(FFT_Engine & engine) :Parent{engine}, Ghat{make_field("Projection Operator", this->projection_container)} {} /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrain:: initialise(FFT_PlanFlags flags) { Parent::initialise(flags); FFT_freqs fft_freqs(this->fft_engine.get_sizes()); for (auto && tup: boost::combine(this->fft_engine, this->Ghat)) { const auto & ccoord = boost::get<0> (tup); auto & G = boost::get<1>(tup); auto xi = fft_freqs.get_unit_xi(ccoord); //! this is simplyfiable usinc Curnier's Méthodes numériques, 6.69(c) - G = Matrices::outer_over(Matrices::I2(), xi*xi.transpose()); + G = Matrices::outer_under(Matrices::I2(), xi*xi.transpose()); // 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); // } // } // } } this->Ghat[0].setZero(); } /* ---------------------------------------------------------------------- */ template void ProjectionFiniteStrain::apply_projection(Field_t & field) { Vector_map field_map{this->fft_engine.fft(field)}; Real factor = this->fft_engine.normalisation(); for (auto && tup: boost::combine(this->Ghat, field_map)) { auto & G{boost::get<0>(tup)}; auto & f{boost::get<1>(tup)}; f = factor * (G*f).eval(); } + this->fft_engine.ifft(field); } template class ProjectionFiniteStrain; template class ProjectionFiniteStrain; } // muSpectre diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 20b2f75..412ca91 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,82 +1,141 @@ /** * file test_projection_finite.cc * * @author Till Junge * * @date 07 Dec 2017 * * @brief tests for standard finite strain projection operator * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include +#include "boost/range/combine.hpp" #include "tests.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain.hh" #include "common/common.hh" #include "common/field_collection.hh" +#include "common/field_map.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(projection_finite_strain); /* ---------------------------------------------------------------------- */ template struct Sizes { }; template<> struct Sizes { constexpr static Ccoord_t get_value() { return Ccoord_t{3, 5};} }; template<> struct Sizes { constexpr static Ccoord_t get_value() { return Ccoord_t{3, 5, 7};} }; /* ---------------------------------------------------------------------- */ template struct ProjectionFixture { using Engine = FFTW_Engine; using Parent = ProjectionFiniteStrain; + constexpr static Dim_t sdim{DimS}; + constexpr static Dim_t mdim{DimM}; ProjectionFixture(): engine{Sizes::get_value()}, - parent(engine){} + projector(engine){} Engine engine; - Parent parent; + Parent projector; }; /* ---------------------------------------------------------------------- */ using fixlist = boost::mpl::list, ProjectionFixture>; /* ---------------------------------------------------------------------- */ BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) { - BOOST_CHECK_NO_THROW(fix::parent.initialise(FFT_PlanFlags::estimate)); + BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate)); } /* ---------------------------------------------------------------------- */ - BOOST_FIXTURE_TEST_CASE_TEMPLATE(Ctest_name, type_name, TL, F) + using F3 = ProjectionFixture; + BOOST_FIXTURE_TEST_CASE(Helmholtz_decomposition, F3) { + // create a field that is explicitely the sum of a gradient and a + // curl, then verify that the projected field corresponds to the + // gradient (without the zero'th component) + using Fields = FieldCollection; + Fields fields{}; + using FieldT = TensorField; + using FieldMap = MatrixFieldMap; + FieldT & f_grad{make_field("gradient", fields)}; + FieldT & f_curl{make_field("curl", fields)}; + FieldT & f_sum{make_field("sum", fields)}; + FieldT & f_var{make_field("working field", fields)}; + + FieldMap grad(f_grad); + FieldMap curl(f_curl); + FieldMap sum(f_sum); + FieldMap var(f_var); + + fields.initialise(Sizes::get_value()); + + for (auto && tup: boost::combine(fields, grad, curl, sum, var)) { + auto & ccoord = boost::get<0>(tup); + auto & g = boost::get<1>(tup); + auto & c = boost::get<2>(tup); + auto & s = boost::get<3>(tup); + auto & v = boost::get<4>(tup); + auto vec = CcoordOps::get_vector(ccoord); + g.row(0) << vec(1)*vec(2), vec(0)*vec(2), vec(1)*vec(2); + c.row(0) << 0, vec(0)*vec(1), -vec(0)*vec(2); + v.row(0) = s.row(0) = g.row(0) + c.row(0); + } + + projector.initialise(FFT_PlanFlags::estimate); + projector.apply_projection(f_var); + + for (auto && tup: boost::combine(fields, grad, curl, sum, var)) { + auto & ccoord = boost::get<0>(tup); + auto & g = boost::get<1>(tup); + auto & c = boost::get<2>(tup); + auto & s = boost::get<3>(tup); + auto & v = boost::get<4>(tup); + auto vec = CcoordOps::get_vector(ccoord); + + Real error = (s-g).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 << "curl :" << std::endl << c << std::endl; + std::cout << std::endl << "sum :" << std::endl << s << 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