diff --git a/language_bindings/python/bind_py_projections.cc b/language_bindings/python/bind_py_projections.cc index acb5d96..f9560e6 100644 --- a/language_bindings/python/bind_py_projections.cc +++ b/language_bindings/python/bind_py_projections.cc @@ -1,198 +1,198 @@ /** * @file bind_py_projections.cc * * @author Till Junge * * @date 18 Jan 2018 * * @brief Python bindings for the Projection operators * * 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/projection_finite_strain.hh" #include "fft/projection_finite_strain_fast.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 #include #include #include #include using namespace muSpectre; namespace py=pybind11; using namespace pybind11::literals; /** * "Trampoline" class for handling the pure virtual methods, see * [http://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python] * for details */ template class PyProjectionBase: public ProjectionBase { public: //! base class using Parent = ProjectionBase; //! field type on which projection is applied using Field_t = typename Parent::Field_t; void apply_projection(Field_t & field) override { PYBIND11_OVERLOAD_PURE (void, Parent, apply_projection, field ); } Eigen::Map get_operator() override { PYBIND11_OVERLOAD_PURE (Eigen::Map, Parent, get_operator ); } }; template void add_proj_helper(py::module & mod, std::string name_start) { using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; using Field_t = typename Proj::Field_t; static_assert(DimS == DimM, "currently only for DimS==DimM"); std::stringstream name{}; name << name_start << '_' << DimS << 'd'; py::class_(mod, name.str().c_str()) #ifdef WITH_MPI .def(py::init([](Ccoord res, Rcoord lengths, const std::string & fft, size_t comm) { if (fft == "fftw") { auto engine = std::make_unique> - (res, Proj::NbComponents, std::move(Communicator(MPI_Comm(comm)))); + (res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); } #else .def(py::init([](Ccoord res, Rcoord lengths, const std::string & fft) { if (fft == "fftw") { auto engine = std::make_unique> - (res, Proj::NbComponents); + (res, Proj::NbComponents()); return Proj(std::move(engine), lengths); } #endif #ifdef WITH_FFTWMPI else if (fft == "fftwmpi") { auto engine = std::make_unique> - (res, Proj::NbComponents, + (res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); } #endif #ifdef WITH_PFFT else if (fft == "pfft") { auto engine = std::make_unique> - (res, Proj::NbComponents, + (res, Proj::NbComponents(), std::move(Communicator(MPI_Comm(comm)))); return Proj(std::move(engine), lengths); } #endif else { throw std::runtime_error("Unknown FFT engine '"+fft+"' specified."); } }), "resolutions"_a, "lengths"_a, #ifdef WITH_MPI "fft"_a="fftw", "communicator"_a=size_t(MPI_COMM_SELF)) #else "fft"_a="fftw") #endif .def("initialise", &Proj::initialise, "flags"_a=FFT_PlanFlags::estimate, "initialises the fft engine (plan the transform)") .def("apply_projection", [](Proj & proj, py::EigenDRef v){ typename FFTEngineBase::GFieldCollection_t coll{}; Eigen::Index subdomain_size = CcoordOps::get_size(proj.get_subdomain_resolutions()); if (v.rows() != DimS*DimM || v.cols() != subdomain_size) { throw std::runtime_error("Expected input array of shape ("+ std::to_string(DimS*DimM)+", "+ std::to_string(subdomain_size)+ "), but input array has shape ("+ std::to_string(v.rows())+", "+ std::to_string(v.cols())+")."); } coll.initialise(proj.get_subdomain_resolutions(), proj.get_subdomain_locations()); Field_t & temp{make_field("temp_field", coll, proj.get_nb_components())}; temp.eigen() = v; proj.apply_projection(temp); return Eigen::ArrayXXd{temp.eigen()}; }) .def("get_operator", &Proj::get_operator) .def("get_formulation", &Proj::get_formulation, "return a Formulation enum indicating whether the projection is small" " or finite strain") .def("get_subdomain_resolutions", &Proj::get_subdomain_resolutions) .def("get_subdomain_locations", &Proj::get_subdomain_locations) .def("get_domain_resolutions", &Proj::get_domain_resolutions) .def("get_domain_lengths", &Proj::get_domain_resolutions); } void add_proj_dispatcher(py::module & mod) { add_proj_helper< ProjectionSmallStrain< twoD, twoD>, twoD>(mod, "ProjectionSmallStrain"); add_proj_helper< ProjectionSmallStrain, threeD>(mod, "ProjectionSmallStrain"); add_proj_helper< ProjectionFiniteStrain< twoD, twoD>, twoD>(mod, "ProjectionFiniteStrain"); add_proj_helper< ProjectionFiniteStrain, threeD>(mod, "ProjectionFiniteStrain"); add_proj_helper< ProjectionFiniteStrainFast< twoD, twoD>, twoD>(mod, "ProjectionFiniteStrainFast"); add_proj_helper< ProjectionFiniteStrainFast, threeD>(mod, "ProjectionFiniteStrainFast"); } void add_projections(py::module & mod) { add_proj_dispatcher(mod); } diff --git a/src/fft/projection_default.hh b/src/fft/projection_default.hh index 18e6029..4a5248e 100644 --- a/src/fft/projection_default.hh +++ b/src/fft/projection_default.hh @@ -1,108 +1,109 @@ /** * @file projection_default.hh * * @author Till Junge * * @date 14 Jan 2018 * * @brief virtual base class for default projection implementation, where the * projection operator is stored as a full fourth-order tensor per * k-space point (as opposed to 'smart' faster implementations, such as * ProjectionFiniteStrainFast * * Copyright (C) 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef PROJECTION_DEFAULT_H #define PROJECTION_DEFAULT_H #include "fft/projection_base.hh" namespace muSpectre { /** * base class to inherit from if one implements a projection * operator that is stored in form of a fourth-order tensor of real * values per k-grid point */ template class ProjectionDefault: public ProjectionBase { public: - constexpr static Dim_t NbComponents{ipow(DimM, 2)}; using Parent = ProjectionBase; //!< base class using Vector_t = typename Parent::Vector_t; //!< to represent fields //! polymorphic FFT pointer type using FFTEngine_ptr = typename Parent::FFTEngine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type //! global field collection using GFieldCollection_t = GlobalFieldCollection; //! local field collection for Fourier-space fields using LFieldCollection_t = LocalFieldCollection; //! Real space second order tensor fields (to be projected) using Field_t = TypedField; //! Fourier-space field containing the projection operator itself using Proj_t = TensorField; //! iterable form of the operator using Proj_map = T4MatrixFieldMap; //! vectorized version of the Fourier-space second-order tensor field using Vector_map = MatrixFieldMap; //! Default constructor ProjectionDefault() = delete; //! Constructor with cell sizes and formulation ProjectionDefault(FFTEngine_ptr engine, Rcoord lengths, Formulation form); //! Copy constructor ProjectionDefault(const ProjectionDefault &other) = delete; //! Move constructor ProjectionDefault(ProjectionDefault &&other) = default; //! Destructor virtual ~ProjectionDefault() = default; //! Copy assignment operator ProjectionDefault& operator=(const ProjectionDefault &other) = delete; //! Move assignment operator ProjectionDefault& operator=(ProjectionDefault &&other) = delete; //! apply the projection operator to a field void apply_projection(Field_t & field) override final; Eigen::Map get_operator() override final; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const override final; + constexpr static Dim_t NbComponents() {return ipow(DimM, 2);} + protected: Proj_t & Gfield; //!< field holding the operator Proj_map Ghat; //!< iterable version of operator private: }; } // muSpectre #endif /* PROJECTION_DEFAULT_H */ diff --git a/src/fft/projection_finite_strain_fast.hh b/src/fft/projection_finite_strain_fast.hh index 444d16f..861ef85 100644 --- a/src/fft/projection_finite_strain_fast.hh +++ b/src/fft/projection_finite_strain_fast.hh @@ -1,114 +1,115 @@ /** * @file projection_finite_strain_fast.hh * * @author Till Junge * * @date 12 Dec 2017 * * @brief Faster alternative to ProjectionFinitestrain * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with GNU Emacs; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #ifndef PROJECTION_FINITE_STRAIN_FAST_H #define PROJECTION_FINITE_STRAIN_FAST_H #include "fft/projection_base.hh" #include "common/common.hh" #include "common/field_collection.hh" #include "common/field_map.hh" namespace muSpectre { /** * replaces `muSpectre::ProjectionFiniteStrain` with a faster and * less memory-hungry alternative formulation. Use this if you don't * have a very good reason not to (and tell me (author) about it, * I'd be interested to hear it). */ template class ProjectionFiniteStrainFast: public ProjectionBase { public: - constexpr static Dim_t NbComponents{ipow(DimM, 2)}; using Parent = ProjectionBase; //!< base class //! polymorphic pointer to FFT engines using FFTEngine_ptr = typename Parent::FFTEngine_ptr; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type //! global field collection (for real-space representations) using GFieldCollection_t = GlobalFieldCollection; //! local field collection (for Fourier-space representations) using LFieldCollection_t = LocalFieldCollection; //! Real space second order tensor fields (to be projected) using Field_t = TypedField; //! Fourier-space field containing the projection operator itself using Proj_t = TensorField; //! iterable form of the operator using Proj_map = MatrixFieldMap; //! iterable Fourier-space second-order tensor field using Grad_map = MatrixFieldMap; //! Default constructor ProjectionFiniteStrainFast() = delete; //! Constructor with fft_engine ProjectionFiniteStrainFast(FFTEngine_ptr engine, Rcoord lengths); //! Copy constructor ProjectionFiniteStrainFast(const ProjectionFiniteStrainFast &other) = delete; //! Move constructor ProjectionFiniteStrainFast(ProjectionFiniteStrainFast &&other) = default; //! Destructor virtual ~ProjectionFiniteStrainFast() = default; //! Copy assignment operator ProjectionFiniteStrainFast& operator=(const ProjectionFiniteStrainFast &other) = delete; //! Move assignment operator ProjectionFiniteStrainFast& operator=(ProjectionFiniteStrainFast &&other) = default; //! initialises the fft engine (plan the transform) virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final; //! apply the projection operator to a field void apply_projection(Field_t & field) override final; Eigen::Map get_operator() override final; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const override final; + constexpr static Dim_t NbComponents(){return ipow(DimM, 2);} + protected: Proj_t & xiField; //!< field of normalised wave vectors Proj_map xis; //!< iterable normalised wave vectors private: }; } // muSpectre #endif /* PROJECTION_FINITE_STRAIN_FAST_H */