diff --git a/src/common/common.hh b/src/common/common.hh index 61da841..2af06e4 100644 --- a/src/common/common.hh +++ b/src/common/common.hh @@ -1,169 +1,192 @@ /** * file common.hh * * @author Till Junge * * @date 01 May 2017 * * @brief Small definitions of commonly used types througout µSpectre * * @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 #ifndef COMMON_H #define COMMON_H namespace muSpectre { //! Eigen uses signed integers for dimensions. For consistency, µSpectre uses //! them througout the code using Dim_t = int;// needs to represent -1 for eigen const Dim_t oneD{1}; const Dim_t twoD{2}; const Dim_t threeD{3}; const Dim_t secondOrder{2}; const Dim_t fourthOrder{4}; + //! Scalar types used for mathematical calculations + using Uint = unsigned int; + using Int = int; + using Real = double; + using Complex = std::complex; + //! Ccoord_t are cell coordinates, i.e. integer coordinates template using Ccoord_t = std::array; - template - std::ostream & operator << (std::ostream & os, const Ccoord_t & index) { + //! Real space coordinates + template + using Rcoord_t = std::array; + + template + std::ostream & operator << (std::ostream & os, + const std::array & index) { os << "("; for (size_t i = 0; i < dim-1; ++i) { os << index[i] << ", "; } os << index.back() << ")"; return os; } - //! Scalar types used for mathematical calculations - using Uint = unsigned int; - using Int = int; - using Real = double; - using Complex = std::complex; + template + Rcoord_t operator/(const Rcoord_t & a, const Rcoord_t & b) { + Rcoord_t retval{a}; + for (size_t i = 0; i < dim; ++i) { + retval[i]/=b[i]; + } + return retval; + } + + template + Rcoord_t operator/(const Rcoord_t & a, const Ccoord_t & b) { + Rcoord_t retval{a}; + for (size_t i = 0; i < dim; ++i) { + retval[i]/=b[i]; + } + return retval; + } //! convenience definitions constexpr Real pi{atan(1.)*4}; //! compile-time potentiation required for field-size computations template constexpr I ipow(I base, I exponent) { static_assert(std::is_integral::value, "Type must be integer"); I retval{1}; for (I i = 0; i < exponent; ++i) { retval *= base; } return retval; } //! continuum mechanics flags enum class Formulation{finite_strain, small_strain}; std::ostream & operator<<(std::ostream & os, Formulation f); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of stress measure they provide, //! and µSpectre will handle conversions enum class StressMeasure { Cauchy, PK1, PK2, Kirchhoff, Biot, Mandel, __nostress__}; std::ostream & operator<<(std::ostream & os, StressMeasure s); /* ---------------------------------------------------------------------- */ //! Material laws can declare which type of strain measure they require and //! µSpectre will provide it enum class StrainMeasure { Gradient, Infinitesimal, GreenLagrange, Biot, Log, Almansi, RCauchyGreen, LCauchyGreen, __nostrain__}; std::ostream & operator<<(std::ostream & os, StrainMeasure s); /* ---------------------------------------------------------------------- */ /** Compile-time functions to set the stress and strain measures stored by mu_spectre depending on the formulation **/ constexpr StrainMeasure get_stored_strain_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StrainMeasure::Gradient; break; } case Formulation::small_strain: { return StrainMeasure::Infinitesimal; break; } default: return StrainMeasure::__nostrain__; break; } } constexpr StressMeasure get_stored_stress_type(Formulation form) { switch (form) { case Formulation::finite_strain: { return StressMeasure::PK1; break; } case Formulation::small_strain: { return StressMeasure::Cauchy; break; } default: return StressMeasure::__nostress__; break; } } /* ---------------------------------------------------------------------- */ /** Compile-time functions to get the stress and strain measures after they may have been modified by choosing a formulation. For instance, a law that expecs a Green-Lagrange strain as input will get the infinitesimal strain tensor instead in a small strain computation **/ constexpr StrainMeasure get_formulation_strain_type(Formulation form, StrainMeasure expected) { switch (form) { case Formulation::finite_strain: { return expected; break; } case Formulation::small_strain: { return get_stored_strain_type(form); break; } default: return StrainMeasure::__nostrain__; break; } } } // muSpectre #ifndef EXPLICITLY_TURNED_ON_CXX17 #include "common/utilities.hh" #endif #endif /* COMMON_H */ diff --git a/src/fft/fft_engine_base.cc b/src/fft/fft_engine_base.cc index e2b258a..b869c65 100644 --- a/src/fft/fft_engine_base.cc +++ b/src/fft/fft_engine_base.cc @@ -1,64 +1,64 @@ /** * file fft_engine_base.cc * * @author Till Junge * * @date 03 Dec 2017 * * @brief implementation for FFT engine base class * * @section LICENCE * * Copyright (C) 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "fft/fft_engine_base.hh" namespace muSpectre { /* ---------------------------------------------------------------------- */ template - FFT_Engine_base::FFT_Engine_base(Ccoord resolutions, Lengths lengths) + FFT_Engine_base::FFT_Engine_base(Ccoord resolutions, Rcoord lengths) :resolutions{resolutions}, lengths{lengths}, work{make_field("work space", work_space_container)}, norm_factor{1./CcoordOps::get_size(resolutions)} {} /* ---------------------------------------------------------------------- */ template void FFT_Engine_base::initialise(FFT_PlanFlags /*plan_flags*/) { this->work_space_container.initialise(); } /* ---------------------------------------------------------------------- */ template size_t FFT_Engine_base::size() const { return CcoordOps::get_size(this->resolutions); } /* ---------------------------------------------------------------------- */ template size_t FFT_Engine_base::workspace_size() const { return this->work_space_container.size(); } template class FFT_Engine_base; template class FFT_Engine_base; template class FFT_Engine_base; } // muSpectre diff --git a/src/fft/fft_engine_base.hh b/src/fft/fft_engine_base.hh index 1b5b6d4..093d805 100644 --- a/src/fft/fft_engine_base.hh +++ b/src/fft/fft_engine_base.hh @@ -1,122 +1,122 @@ /** * file fft_engine_base.hh * * @author Till Junge * * @date 01 Dec 2017 * * @brief Interface for FFT engines * * @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 "common/common.hh" #include "common/field_collection.hh" #ifndef FFT_ENGINE_BASE_H #define FFT_ENGINE_BASE_H namespace muSpectre { enum class FFT_PlanFlags {estimate, measure, patient}; template class FFT_Engine_base { public: constexpr static Dim_t sdim{DimS}; constexpr static Dim_t mdim{DimM}; using Ccoord = Ccoord_t; - using Lengths = std::array; + using Rcoord = std::array; using GFieldCollection_t = FieldCollection; using LFieldCollection_t = FieldCollection; using Field_t = TensorField; using Workspace_t = TensorField; using iterator = typename LFieldCollection_t::iterator; //! Default constructor FFT_Engine_base() = delete; //! Constructor with system resolutions - FFT_Engine_base(Ccoord resolutions, Lengths lengths); + 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) noexcept = default; //! Destructor virtual ~FFT_Engine_base() noexcept = 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) noexcept = default; // compute the plan, etc virtual void initialise(FFT_PlanFlags /*plan_flags*/); //! forward transform (dummy for interface) virtual Workspace_t & fft(const Field_t & /*field*/) = 0; //! inverse transform (dummy for interface) virtual void ifft(Field_t & /*field*/) const = 0; /** * iterators over only thos pixels that exist in frequency space * (i.e. about half of all pixels, see rfft) */ inline iterator begin() {return this->work_space_container.begin();} inline iterator end() {return this->work_space_container.end();} //! nb of pixels (mostly for debugging) size_t size() const; size_t workspace_size() const; //! const Ccoord & get_resolutions() const {return this->resolutions;} - const Lengths & get_lengths() const {return this->lengths;} + const Rcoord & get_lengths() const {return this->lengths;} LFieldCollection_t & get_field_collection() { return this->work_space_container;} //! 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: LFieldCollection_t work_space_container{}; const Ccoord resolutions; - const Lengths lengths; + const Rcoord lengths; Workspace_t & work; const Real norm_factor; private: }; } // muSpectre #endif /* FFT_ENGINE_BASE_H */ diff --git a/src/fft/fftw_engine.cc b/src/fft/fftw_engine.cc index ccc20a9..e937117 100644 --- a/src/fft/fftw_engine.cc +++ b/src/fft/fftw_engine.cc @@ -1,128 +1,131 @@ /** * file fftw_engine.cc * * @author Till Junge * * @date 03 Dec 2017 * * @brief implements the fftw engine * * @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 "common/ccoord_operations.hh" #include "fft/fftw_engine.hh" namespace muSpectre { template - FFTW_Engine::FFTW_Engine(Ccoord resolutions, Lengths lengths) + FFTW_Engine::FFTW_Engine(Ccoord resolutions, Rcoord lengths) :Parent{resolutions, lengths}, hermitian_resolutions{CcoordOps::get_hermitian_sizes(resolutions)} { for (auto && pixel: CcoordOps::Pixels(this->hermitian_resolutions)) { this->work_space_container.add_pixel(pixel); } } /* ---------------------------------------------------------------------- */ template void FFTW_Engine::initialise(FFT_PlanFlags plan_flags) { Parent::initialise(plan_flags); const int & rank = DimS; std::array narr; const int * const n = &narr[0]; std::copy(this->resolutions.begin(), this->resolutions.end(), narr.begin()); int howmany = Field_t::nb_components; //temporary buffer for plan Real * r_work_space = fftw_alloc_real(CcoordOps::get_size(this->resolutions) *howmany); Real * in = r_work_space; const int * const inembed = nullptr;//nembed are tricky: they refer to physical layout int istride = howmany; int idist = 1; fftw_complex * out = reinterpret_cast(this->work.data()); const int * const onembed = nullptr; int ostride = howmany; int odist = idist; unsigned int flags; if (plan_flags == FFT_PlanFlags::estimate) { flags = FFTW_ESTIMATE; } if (plan_flags == FFT_PlanFlags::measure) { flags = FFTW_MEASURE; } if (plan_flags == FFT_PlanFlags::patient) { flags = FFTW_PATIENT; } this->plan_fft = fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist, out, onembed, ostride, odist, FFTW_PRESERVE_INPUT | flags); if (this->plan_fft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_complex * i_in = reinterpret_cast(this->work.data()); Real * i_out = r_work_space; this->plan_ifft = fftw_plan_many_dft_c2r(rank, n, howmany, i_in, inembed, istride, idist, i_out, onembed, ostride, odist, flags); if (this->plan_ifft == nullptr) { throw std::runtime_error("Plan failed"); } fftw_free(r_work_space); } /* ---------------------------------------------------------------------- */ template FFTW_Engine::~FFTW_Engine() noexcept { fftw_destroy_plan(this->plan_fft); fftw_destroy_plan(this->plan_ifft); fftw_cleanup(); } /* ---------------------------------------------------------------------- */ template typename FFTW_Engine::Workspace_t & FFTW_Engine::fft (const Field_t & field) { fftw_execute_dft_r2c(this->plan_fft, const_cast(field.data()), reinterpret_cast(this->work.data())); return this->work; } /* ---------------------------------------------------------------------- */ template void FFTW_Engine::ifft (Field_t & field) const { + if (field.size() != CcoordOps::get_size(this->resolutions)) { + throw std::runtime_error("size mismatch"); + } fftw_execute_dft_c2r(this->plan_ifft, reinterpret_cast(this->work.data()), field.data()); } template class FFTW_Engine; template class FFTW_Engine; template class FFTW_Engine; } // muSpectre diff --git a/src/fft/fftw_engine.hh b/src/fft/fftw_engine.hh index 46c85e0..e0066f8 100644 --- a/src/fft/fftw_engine.hh +++ b/src/fft/fftw_engine.hh @@ -1,87 +1,87 @@ /** * file fftw_engine.hh * * @author Till Junge * * @date 03 Dec 2017 * * @brief FFT engine using FFTW * * @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 "fft/fft_engine_base.hh" #ifndef FFTW_ENGINE_H #define FFTW_ENGINE_H namespace muSpectre { template class FFTW_Engine: public FFT_Engine_base { public: using Parent = FFT_Engine_base; using Ccoord = typename Parent::Ccoord; - using Lengths = typename Parent::Lengths; + using Rcoord = typename Parent::Rcoord; using Workspace_t = typename Parent::Workspace_t; using Field_t = typename Parent::Field_t; //! Default constructor FFTW_Engine() = delete; //! Constructor with system resolutions - FFTW_Engine(Ccoord resolutions, Lengths lengths); + FFTW_Engine(Ccoord resolutions, Rcoord lengths); //! Copy constructor FFTW_Engine(const FFTW_Engine &other) = delete; //! Move constructor FFTW_Engine(FFTW_Engine &&other) = default; //! Destructor virtual ~FFTW_Engine() noexcept; //! Copy assignment operator FFTW_Engine& operator=(const FFTW_Engine &other) = delete; //! Move assignment operator FFTW_Engine& operator=(FFTW_Engine &&other) = default; // compute the plan, etc virtual void initialise(FFT_PlanFlags plan_flags) override; //! forward transform virtual Workspace_t & fft(const Field_t & field) override; //! inverse transform virtual void ifft(Field_t & field) const override; protected: Ccoord hermitian_resolutions; fftw_plan plan_fft{}; fftw_plan plan_ifft{}; private: }; } // muSpectre #endif /* FFTW_ENGINE_H */ diff --git a/tests/test_ccoord_operations.cc b/tests/test_ccoord_operations.cc index 50cc5a4..002ed00 100644 --- a/tests/test_ccoord_operations.cc +++ b/tests/test_ccoord_operations.cc @@ -1,120 +1,120 @@ /** * file test_ccoord_operations.cc * * @author Till Junge * * @date 03 Dec 2017 * * @brief tests for cell coordinate operations * * @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 "common/common.hh" #include "common/ccoord_operations.hh" #include "common/test_goodies.hh" #include "tests.hh" namespace muSpectre { BOOST_AUTO_TEST_SUITE(ccoords_operations); BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_cube, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); Ccoord ref_cube; for (Dim_t i = 0; i < dim; ++i) { ref_cube[i] = size; } BOOST_CHECK_EQUAL_COLLECTIONS(ref_cube.begin(), ref_cube.end(), cube.begin(), cube.end()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_hermitian, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); constexpr Ccoord herm = CcoordOps::get_hermitian_sizes(cube); Ccoord ref_cube = cube; ref_cube.back() = (cube.back() + 1) / 2; BOOST_CHECK_EQUAL_COLLECTIONS(ref_cube.begin(), ref_cube.end(), herm.begin(), herm.end()); } BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_size, Fix, testGoodies::dimlist, Fix) { constexpr auto dim{Fix::dim}; using Ccoord = Ccoord_t; constexpr Dim_t size{5}; constexpr Ccoord cube = CcoordOps::get_cube(size); BOOST_CHECK_EQUAL(CcoordOps::get_size(cube), ipow(size, dim)); } BOOST_AUTO_TEST_CASE(vector_test) { constexpr Ccoord_t c3{1, 2, 3}; constexpr Ccoord_t c2{c3[0], c3[1]}; - constexpr std::array s3{1.3, 2.8, 5.7}; - constexpr std::array s2{s3[0], s3[1]}; + constexpr Rcoord_t s3{1.3, 2.8, 5.7}; + constexpr Rcoord_t s2{s3[0], s3[1]}; Eigen::Matrix v2; v2 << s3[0], s3[1]; Eigen::Matrix v3; v3 << s3[0], s3[1], s3[2]; auto vec2{CcoordOps::get_vector(c2, v2(1))}; auto vec3{CcoordOps::get_vector(c3, v3(1))}; for (Dim_t i = 0; i < twoD; ++i) { BOOST_CHECK_EQUAL(c2[i]*v2(1), vec2[i]); } for (Dim_t i = 0; i < threeD; ++i) { BOOST_CHECK_EQUAL(c3[i]*v3(1), vec3[i]); } vec2 = CcoordOps::get_vector(c2, v2); vec3 = CcoordOps::get_vector(c3, v3); for (Dim_t i = 0; i < twoD; ++i) { BOOST_CHECK_EQUAL(c2[i]*v2(i), vec2[i]); BOOST_CHECK_EQUAL(vec2[i], CcoordOps::get_vector(c2, s2)[i]); } for (Dim_t i = 0; i < threeD; ++i) { BOOST_CHECK_EQUAL(c3[i]*v3(i), vec3[i]); BOOST_CHECK_EQUAL(vec3[i], CcoordOps::get_vector(c3, s3)[i]); } } BOOST_AUTO_TEST_SUITE_END(); } // muSpectre diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc index 04d5bc6..53f3416 100644 --- a/tests/test_projection_finite.cc +++ b/tests/test_projection_finite.cc @@ -1,213 +1,237 @@ /** * 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 #include #include #include "tests.hh" #include "fft/fftw_engine.hh" #include "fft/projection_finite_strain.hh" #include "fft/fft_utils.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_resolution() { - return Ccoord_t{13, 15};} - constexpr static std::array get_lengths() { - return std::array{3.4, 5.8};} + return Ccoord_t{3, 5};} + constexpr static Rcoord_t get_lengths() { + return Rcoord_t{3.4, 5.8};} }; template<> struct Sizes { constexpr static Ccoord_t get_resolution() { - return Ccoord_t{13, 15, 17};} - constexpr static std::array get_lengths() { - return std::array{3.4, 5.8, 6.7};} + return Ccoord_t{3, 5, 7};} + constexpr static Rcoord_t get_lengths() { + return Rcoord_t{3.4, 5.8, 6.7};} }; + template + struct Squares { + }; + template<> + struct Squares { + constexpr static Ccoord_t get_resolution() { + return Ccoord_t{5, 5};} + constexpr static Rcoord_t get_lengths() { + return Rcoord_t{5, 5};} + }; + template<> + struct Squares { + constexpr static Ccoord_t get_resolution() { + return Ccoord_t{7, 7, 7};} + constexpr static Rcoord_t get_lengths() { + return Rcoord_t{7, 7, 7};} + }; + /* ---------------------------------------------------------------------- */ - template + 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_resolution(), - Sizes::get_lengths()}, + ProjectionFixture(): engine{SizeGiver::get_resolution(), + SizeGiver::get_lengths()}, projector(engine){} Engine engine; Parent projector; }; /* ---------------------------------------------------------------------- */ - using fixlist = boost::mpl::list, - ProjectionFixture>; + using fixlist = + boost::mpl::list>, + ProjectionFixture>, + ProjectionFixture>, + ProjectionFixture>>; /* ---------------------------------------------------------------------- */ 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 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(Sizes::get_resolution()); + fields.initialise(fix::engine.get_resolutions()); FFT_freqs freqs{fix::engine.get_resolutions(), fix::engine.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::engine.get_lengths()[i]; ; } for (auto && tup: boost::combine(fields, grad, var)) { auto & ccoord = boost::get<0>(tup); auto & g = boost::get<1>(tup); auto & v = boost::get<2>(tup); - Vector vec = CcoordOps::get_vector(ccoord, fix::engine.get_lengths()); - g.col(0) << 2*pi * k * cos(k.dot(vec)); - v.col(0) = g.col(0); + Vector vec = CcoordOps::get_vector(ccoord, + fix::engine.get_lengths()/ + fix::engine.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: boost::combine(fields, grad, var)) { auto & ccoord = boost::get<0>(tup); auto & g = boost::get<1>(tup); auto & v = boost::get<2>(tup); - auto vec = CcoordOps::get_vector(ccoord); - + Vector vec = CcoordOps::get_vector(ccoord, + fix::engine.get_lengths()/ + fix::engine.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; } } } - /* ---------------------------------------------------------------------- */ - 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) - constexpr Dim_t sdim{F3::sdim}, mdim{F3::mdim}; - 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_resolution()); - - 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 & x{vec(0)}; - Real & y{vec(1)}; - Real & z{vec(2)}; - g.col(0) << y*z, x*z, y*z; - c.col(0) << 0, x*y, -x*z; - v.col(0) = s.col(0) = g.col(0) + c.col(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; - } - } - } + // /* ---------------------------------------------------------------------- */ + // 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) + // constexpr Dim_t sdim{F3::sdim}, mdim{F3::mdim}; + // 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(F3::engine.get_resolutions()); + + // 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 & x{vec(0)}; + // Real & y{vec(1)}; + // Real & z{vec(2)}; + // g.row(0) << y*z, x*z, y*z; + // c.row(0) << 0, x*y, -x*z; + // 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