diff --git a/src/fft/projection_finite_strain.cc b/src/fft/projection_finite_strain.cc index bc020c5..292ab27 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) :Parent{std::move(engine), Formulation::finite_strain} { - for (auto res: this->fft_engine->get_resolutions()) { + 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->fft_engine->get_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_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 89b2cf2..08337cf 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) :Parent{std::move(engine), Formulation::finite_strain}, xiField{make_field("Projection Operator", this->projection_container)}, xis(xiField) { - for (auto res: this->fft_engine->get_resolutions()) { + 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->fft_engine->get_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_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 f045ce6..7279fc6 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) : Parent{std::move(engine), Formulation::small_strain} { - for (auto res: this->fft_engine->get_resolutions()) { + 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->fft_engine->get_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_locations() == Ccoord{}) { this->Ghat[0].setZero(); } } template class ProjectionSmallStrain; template class ProjectionSmallStrain; } // muSpectre