diff --git a/src/cell/cell_base.cc b/src/cell/cell_base.cc index 9801afa..3e69ebd 100644 --- a/src/cell/cell_base.cc +++ b/src/cell/cell_base.cc @@ -1,405 +1,426 @@ /** * @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), 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()} + projection{std::move(projection_)} + { } + + /** + * turns out that the default move container in combination with + * clang segfaults under certain (unclear) cicumstances, because the + * move constructor of the optional appears to be busted in gcc + * 7.2. Copying it (K) instead of moving it fixes the issue, and + * since it is a reference, the cost is practically nil + */ + template + CellBase::CellBase(CellBase && other): + subdomain_resolutions{std::move(other.subdomain_resolutions)}, + subdomain_locations{std::move(other.subdomain_locations)}, + domain_resolutions{std::move(other.domain_resolutions)}, + pixels{std::move(other.pixels)}, + domain_lengths{std::move(other.domain_lengths)}, + fields{std::move(other.fields)}, + F{other.F}, + P{other.P}, + K{other.K}, // this seems to segfault under clang if it's not a move + materials{std::move(other.materials)}, + projection{std::move(other.projection)} { } /* ---------------------------------------------------------------------- */ template typename CellBase::Material_t & CellBase::add_material(Material_ptr mat) { this->materials.push_back(std::move(mat)); return *this->materials.back(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_strain_vector() -> Vector_ref { return this->get_strain().eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase::get_stress_vector() const -> ConstVector_ref { return this->get_stress().eigenvec(); } /* ---------------------------------------------------------------------- */ template void CellBase:: set_uniform_strain(const Eigen::Ref & strain) { this->F.get_map() = strain; } /* ---------------------------------------------------------------------- */ template auto CellBase::evaluate_stress() -> ConstVector_ref { if (not this->initialised) { this->initialise(); } for (auto & mat: this->materials) { - mat->compute_stresses(this->F, this->P, this->form); + mat->compute_stresses(this->F, this->P, this->get_formulation()); } return this->P.const_eigenvec(); } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_stress_tangent() -> std::array { if (not this->initialised) { this->initialise(); } constexpr bool create_tangent{true}; this->get_tangent(create_tangent); for (auto & mat: this->materials) { mat->compute_stresses_tangent(this->F, this->P, this->K.value(), - this->form); + this->get_formulation()); } const TangentField_t & k = this->K.value(); return std::array{ this->P.const_eigenvec(), k.const_eigenvec()}; } /* ---------------------------------------------------------------------- */ template auto CellBase:: evaluate_projected_directional_stiffness (Eigen::Ref delF) -> Vector_ref { // the following const_cast should be safe, as long as the // constructed delF_field is const itself const TypedField delF_field ("Proxied raw memory for strain increment", *this->fields, Eigen::Map(const_cast(delF.data()), delF.size()), this->F.get_nb_components()); if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_nb_dof() << " but I got " << delF.size(); throw std::runtime_error(err.str()); } const std::string out_name{"δP; temp output for directional stiffness"}; auto & delP = this->get_managed_field(out_name); auto Kmap{this->K.value().get().get_map()}; auto delPmap{delP.get_map()}; MatrixFieldMap delFmap(delF_field); for (auto && tup: akantu::zip(Kmap, delFmap, delPmap)) { auto & k = std::get<0>(tup); auto & df = std::get<1>(tup); auto & dp = std::get<2>(tup); dp = Matrices::tensmult(k, df); } return Vector_ref(this->project(delP).data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template std::array CellBase::get_strain_shape() const { return this->projection->get_strain_shape(); } /* ---------------------------------------------------------------------- */ template void CellBase::apply_projection(Eigen::Ref vec) { TypedField field("Proxy for projection", *this->fields, vec, this->F.get_nb_components()); this->projection->apply_projection(field); } /* ---------------------------------------------------------------------- */ 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); + this->get_formulation()); } 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::Vector_ref CellBase::directional_stiffness_vec(const Eigen::Ref &delF) { if (!this->K) { throw std::runtime_error ("currently only implemented for cases where a stiffness matrix " "exists"); } if (delF.size() != this->get_nb_dof()) { std::stringstream err{}; err << "input should be of size ndof = ¶(" << this->subdomain_resolutions << ") × " << DimS << "² = "<< this->get_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); Vector_ref(in_tempref.data(), this->get_nb_dof()) = delF; this->directional_stiffness(this->K.value(), in_tempref, out_tempref); return Vector_ref(out_tempref.data(), this->get_nb_dof()); } /* ---------------------------------------------------------------------- */ template Eigen::ArrayXXd CellBase:: directional_stiffness_with_copy (Eigen::Ref delF) { if (!this->K) { throw std::runtime_error ("currently 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(); for (auto && mat: this->materials) { mat->initialise(); } // 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::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 auto CellBase::get_adaptor() -> 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 f6849a2..ac991ab 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,511 +1,511 @@ /** * @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 #include namespace muSpectre { /** * Cell adaptors implement the matrix-vector multiplication and * allow the system to be used like a sparse matrix in * conjugate-gradient-type solvers */ template class CellAdaptor; /** * Base class for cells that is not templated and therefore can be * in solvers that see cells as runtime-polymorphic objects. This * allows the use of standard * (i.e. spectral-method-implementation-agnostic) solvers, as for * instance the scipy solvers */ class Cell { public: //! sparse matrix emulation using Adaptor = CellAdaptor; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = Eigen::Matrix; //! dynamic matrix type for setting strains using Matrix_t = Eigen::Matrix; //! ref to constant vector using ConstVector_ref = Eigen::Map; //! output vector reference for solvers using Vector_ref = Eigen::Map; //! Default constructor Cell() = default; //! Copy constructor Cell(const Cell &other) = default; //! Move constructor Cell(Cell &&other) = default; //! Destructor virtual ~Cell() = default; //! Copy assignment operator Cell& operator=(const Cell &other) = default; //! Move assignment operator Cell& operator=(Cell &&other) = default; //! for handling double initialisations right bool is_initialised() const {return this->initialised;} //! returns the number of degrees of freedom in the cell virtual Dim_t get_nb_dof() const = 0; //! return the communicator object virtual const Communicator & get_communicator() const = 0; /** * formulation is hard set by the choice of the projection class */ virtual const Formulation & get_formulation() const = 0; /** * returns the material dimension of the problem */ virtual Dim_t get_material_dim() const = 0; /** * 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) */ virtual std::array get_strain_shape() const = 0; /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() = 0; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const = 0; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() = 0; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() = 0; /** * applies the projection operator in-place on the input vector */ virtual void apply_projection(Eigen::Ref vec) = 0; /** * freezes all the history variables of the materials */ virtual void save_history_variables() = 0; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) = 0; /** * set uniform strain (typically used to initialise problems */ virtual void set_uniform_strain(const Eigen::Ref &) = 0; //! get a sparse matrix view on the cell virtual Adaptor get_adaptor() = 0; protected: bool initialised{false}; //!< to handle double initialisation right private: }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase: public Cell { public: using Parent = Cell; 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; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = typename Parent::Vector_t; //! ref to constant vector using ConstVector_ref = typename Parent::ConstVector_ref; //! output vector reference for solvers using Vector_ref = typename Parent::Vector_ref; //! sparse matrix emulation using Adaptor = Parent::Adaptor; //! 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; + CellBase(CellBase &&other); //! 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); /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() override; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const override; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() override; /** * evaluates and returns the stress and stiffness for the currently set strain */ virtual std::array evaluate_stress_tangent() override; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc. (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness (Eigen::Ref delF) override; //! return the template param DimM (required for polymorphic use of `Cell` Dim_t get_material_dim() const override final {return DimM;} /** * 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; /** * applies the projection operator in-place on the input vector */ void apply_projection(Eigen::Ref vec) override final; /** * set uniform strain (typically used to initialise problems */ void set_uniform_strain(const Eigen::Ref &) override; /** * 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 */ Vector_ref directional_stiffness_vec(const Eigen::Ref & 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); /** * 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() override final; 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_domain_lengths() const {return this->domain_lengths;} /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const override final { return this->projection->get_formulation();} /** * 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() override; //! returns the number of degrees of freedom in the cell Dim_t get_nb_dof() const override {return this->size()*ipow(DimS, 2);}; //! return the communicator object virtual const Communicator & get_communicator() const override { 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 & 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 - const Formulation form; //!< formulation for solution + private: }; /** * lightweight resource handle wrapping a `muSpectre::Cell` 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.get_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. evaluate_projected_directional_stiffness(rhs); } }; } } #endif /* CELL_BASE_H */