diff --git a/src/model/elasto_plastic/residual.cpp b/src/model/elasto_plastic/residual.cpp index dc31807..6b123a6 100644 --- a/src/model/elasto_plastic/residual.cpp +++ b/src/model/elasto_plastic/residual.cpp @@ -1,122 +1,138 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "residual.hh" +#include "grid_view.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ template ResidualTemplate::ResidualTemplate(Model* model, Real sigma_0, Real h) : Residual(model), hardening(model, sigma_0, h) { // Registering operators for residual calculation model->registerIntegralOperator>("mindlin_gradient"); model->registerIntegralOperator>("boussinesq_gradient"); // Registering operator for displacement calculation model->registerIntegralOperator>("mindlin"); model->registerIntegralOperator>("boussinesq"); // Initializing grids for (auto grid : {&strain, &stress, &residual, &tmp}) { grid->setNbComponents(trait::components * trait::components); grid->resize(model->getDiscretization()); } } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeResidual(GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computePlasticIncrement(residual, strain, inc); + this->updateFilter(residual); this->model->applyElasticity(residual, residual); - integralOperator("mindlin_gradient").apply(residual, residual); + integralOperator("mindlin_gradient") + .applyIf(residual, residual, plastic_filter); integralOperator("boussinesq_gradient") .apply(this->model->getTraction(), tmp); residual -= strain_increment; residual += tmp; // Symmetrize gradient Loop::stridedLoop( [](MatrixProxy&& gradient) { - // Explicit loop here to avoid temporaries + // Explicit loop here to avoid temporaries for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { gradient(i, j) = gradient(j, i) = 0.5 * (gradient(i, j) + gradient(j, i)); } } }, residual); } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeStress(GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computeStress(stress, strain, inc); // integralOperator("boussinesq_gradient") // .apply(this->model->getTraction(), tmp); // this->model->applyElasticity(tmp, tmp); // stress += tmp; } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::updateState( GridBase& converged_strain_increment) { computeStress(converged_strain_increment); Grid& inc = dynamic_cast(converged_strain_increment); hardening.template computePlasticIncrement(residual, strain, inc); strain += converged_strain_increment; // Computing displacements - integralOperator("mindlin").apply(hardening.getPlasticStrain(), - this->model->getDisplacement()); + integralOperator("mindlin").applyIf(hardening.getPlasticStrain(), + this->model->getDisplacement(), + plastic_filter); Grid disp_tmp(this->model->getDiscretization(), trait::components); integralOperator("boussinesq").apply(this->model->getTraction(), disp_tmp); this->model->getDisplacement() += disp_tmp; } /* -------------------------------------------------------------------------- */ template void ResidualTemplate::computeResidualDisplacement( GridBase& strain_increment) { Grid& inc = dynamic_cast(strain_increment); hardening.template computePlasticIncrement(residual, strain, inc); this->model->applyElasticity(residual, residual); - integralOperator("mindlin").apply(residual, this->model->getDisplacement()); + integralOperator("mindlin").applyIf(residual, this->model->getDisplacement(), + plastic_filter); +} + +/* -------------------------------------------------------------------------- */ +template +void ResidualTemplate::updateFilter( + Grid& plastic_strain_increment) { + for (UInt i : Loop::range(plastic_strain_increment.sizes().front())) { + auto slice = make_view(plastic_strain_increment, i); + if (slice.dot(slice) / slice.dataSize() > 1e-14) + plastic_layers.insert(i); + } } /* -------------------------------------------------------------------------- */ template class ResidualTemplate; /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/elasto_plastic/residual.hh b/src/model/elasto_plastic/residual.hh index 24b3873..249bda3 100644 --- a/src/model/elasto_plastic/residual.hh +++ b/src/model/elasto_plastic/residual.hh @@ -1,122 +1,132 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef RESIDUAL_HH #define RESIDUAL_HH /* -------------------------------------------------------------------------- */ +#include "boussinesq.hh" #include "isotropic_hardening.hh" -#include "model_type.hh" #include "mindlin.hh" -#include "boussinesq.hh" +#include "model_type.hh" +/* -------------------------------------------------------------------------- */ +#include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ /** * @brief Residual manager */ class Residual { public: /// Constructor Residual(Model* model) : model(model) {} /// Destructor virtual ~Residual() = default; // Pure virtual methods public: /// Compute the residual vector for a given strain increment virtual void computeResidual(GridBase& strain_increment) = 0; /// Compute the stresses for a given strain increment virtual void computeStress(GridBase& strain_increment) = 0; /// Update the plastic state virtual void updateState(GridBase& converged_strain_increment) = 0; /// Get residual vector virtual const GridBase& getVector() const = 0; /// Get plastic strain virtual const GridBase& getPlasticStrain() const = 0; /// Get stresses virtual const GridBase& getStress() const = 0; /// Compute displacement virtual void computeResidualDisplacement(GridBase& strain_increment) = 0; // Accessors public: Model& getModel() { return *model; } protected: Model* model; }; /** * @brief Templated residual manager */ template class ResidualTemplate : public Residual { using trait = model_type_traits; static constexpr UInt dim = trait::dimension; public: - ResidualTemplate(Model * model, Real sigma_0, Real h); + ResidualTemplate(Model* model, Real sigma_0, Real h); // Implementation public: /// Compute the residual vector for a given strain increment void computeResidual(GridBase& strain_increment) override; /// Compute the stresses for a given strain increment void computeStress(GridBase& strain_increment) override; /// Update the plastic state void updateState(GridBase& converged_strain_increment) override; /// Compute displacement void computeResidualDisplacement(GridBase& strain_increment) override; // Accessors public: /// Get residual vector const GridBase& getVector() const override { return residual; } /// Get plastic strain const GridBase& getPlasticStrain() const override { return hardening.getPlasticStrain(); } /// Get stresses const GridBase& getStress() const override { return stress; } private: /// Convenience function const IntegralOperator& integralOperator(const std::string& name) const { return *this->model->getIntegralOperator(name); } + /// Add non-zero layers of plastic strain into the filter + void updateFilter(Grid& plastic_strain_increment); + protected: IsotropicHardening hardening; Grid strain, stress, residual, tmp; + + std::unordered_set plastic_layers; + std::function plastic_filter = [this](UInt layer) { + return plastic_layers.find(layer) != plastic_layers.end(); + }; }; /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif // RESIDUAL_HH diff --git a/src/model/integral_operator.hh b/src/model/integral_operator.hh index ce47a68..6405c64 100644 --- a/src/model/integral_operator.hh +++ b/src/model/integral_operator.hh @@ -1,103 +1,109 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef INTEGRAL_OPERATOR_HH #define INTEGRAL_OPERATOR_HH /* -------------------------------------------------------------------------- */ #include "grid_base.hh" #include "model_type.hh" -#include #include +#include #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ class Model; /** * @brief Generic class for integral operators */ class IntegralOperator { public: /// Kind of operator enum kind { neumann, dirichlet }; public: /// Constructor IntegralOperator(Model* model) : model(model) {} /// Virtual destructor for subclasses virtual ~IntegralOperator() = default; /// Apply operator on input virtual void apply(GridBase& input, GridBase& output) const = 0; + /// Apply operator on filtered input + virtual void applyIf(GridBase& input, GridBase& output, + std::function) const { + apply(input, output); + } + /// Update any data dependent on model parameters virtual void updateFromModel() = 0; /// Get model const Model& getModel() const { return *model; } /// Kind virtual kind getKind() const = 0; /// Type virtual model_type getType() const = 0; /// Norm virtual Real getOperatorNorm() { TAMAAS_EXCEPTION("operator does not implement norm"); } protected: Model* model = nullptr; }; /* -------------------------------------------------------------------------- */ /* Printing IntegralOperator::kind to string */ /* -------------------------------------------------------------------------- */ #define INTEGRAL_KINDS (neumann)(dirichlet) inline std::ostream& operator<<(std::ostream& o, const IntegralOperator::kind& val) { switch (val) { #define PRINT_INTEGRAL_KIND(r, data, kind) \ case data::kind: \ o << BOOST_PP_STRINGIZE(kind); \ break; BOOST_PP_SEQ_FOR_EACH(PRINT_INTEGRAL_KIND, IntegralOperator, INTEGRAL_KINDS); #undef PRINT_INTEGRAL_KIND } return o; } #undef INTEGRAL_KINDS /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif // INTEGRAL_OPERATOR_HH diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp index 4945000..2fe9c73 100644 --- a/src/model/kelvin.cpp +++ b/src/model/kelvin.cpp @@ -1,128 +1,136 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "kelvin.hh" #include "elasto_plastic_model.hh" #include "kelvin_helper.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Constructors */ /* -------------------------------------------------------------------------- */ template <> Kelvin::Kelvin(Model* model) : VolumePotential(model) { initialize(trait::components, trait::components); } template <> Kelvin::Kelvin(Model* model) : VolumePotential(model) { initialize(trait::components * trait::components, trait::components); } template <> Kelvin::Kelvin(Model* model) : VolumePotential(model) { initialize(trait::components * trait::components, trait::components * trait::components); } /* -------------------------------------------------------------------------- */ /* Operator implementation */ /* -------------------------------------------------------------------------- */ template -void Kelvin::apply(GridBase& source, - GridBase& out) const { +void Kelvin::applyIf(GridBase& source, GridBase& out, + filter_t pred) const { constexpr UInt derivative = order - 2; Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus(); influence::Kelvin kelvin(mu, nu); - auto apply = [this, &kelvin](UInt node, - decltype(this->source_buffers)& source_buffers, - decltype(this->disp_buffer)& out_buffer) { + auto apply = [this, &kelvin, &pred]( + UInt node, decltype(this->source_buffers)& source_buffers, + decltype(this->disp_buffer)& out_buffer) { constexpr UInt derivative = order - 2; const Real L = this->model->getSystemSize().front(); const UInt N = this->model->getDiscretization().front(); const Real dl = L / (N - 1); const Real xi = node * dl; detail::LinearElement element{ thrust::make_pair(node, xi), kelvin, source_buffers, out_buffer, this->wavevectors}; out_buffer = 0; + // Simple condition to ensure correct integration + auto filter_out_elem = [&pred](UInt e) { + return not pred(e) and not pred(e + 1); + }; + // Loop over elements for (UInt e : Loop::range(N - 1)) { + if (filter_out_elem(e)) + continue; + const thrust::pair node_positions{e * dl, (e + 1) * dl}; element.integrateElement(e, node_positions); } // Setting fundamental frequency to zero typename detail::KelvinTrait< influence::Kelvin>::out_t out_fundamental(out_buffer(0)); out_fundamental = 0; }; - this->fourierApply(apply, source, out); + this->fourierApplyIf(apply, source, out, pred); } /* -------------------------------------------------------------------------- */ template void Kelvin::initialize(UInt source_components, - UInt out_components) { + UInt out_components) { // Copy horizontal sizes std::array sizes; std::copy(this->model->getDiscretization().begin() + 1, this->model->getDiscretization().end(), sizes.begin()); auto hermitian_sizes = GridHermitian::hermitianDimensions( sizes); /// Initializing buffers this->source_buffers.resize(this->model->getDiscretization()[0]); std::for_each(this->source_buffers.begin(), this->source_buffers.end(), [&](typename parent::BufferType& buffer) { buffer.setNbComponents(source_components); buffer.resize(hermitian_sizes); }); this->disp_buffer.setNbComponents(out_components); this->disp_buffer.resize(hermitian_sizes); } /* -------------------------------------------------------------------------- */ /* Template instanciation */ /* -------------------------------------------------------------------------- */ template class Kelvin; template class Kelvin; template class Kelvin; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/model/kelvin.hh b/src/model/kelvin.hh index 7e5525b..f6535e3 100644 --- a/src/model/kelvin.hh +++ b/src/model/kelvin.hh @@ -1,61 +1,64 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef KELVIN_HH #define KELVIN_HH /* -------------------------------------------------------------------------- */ -#include "volume_potential.hh" -#include "model_type.hh" #include "grid_hermitian.hh" +#include "model_type.hh" +#include "volume_potential.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Kelvin tensor */ template class Kelvin : public VolumePotential { static_assert(type == model_type::volume_1d || type == model_type::volume_2d, "Only volume types are supported"); using trait = model_type_traits; using parent = VolumePotential; +protected: + using filter_t = typename parent::filter_t; + public: /// Constructor - Kelvin(Model * model); + Kelvin(Model* model); protected: void initialize(UInt source_components, UInt out_components); public: /// Apply the Kelvin-tensor_order operator - void apply(GridBase& source, GridBase& out) const override; + void applyIf(GridBase& source, GridBase& out, + filter_t pred) const override; }; +__END_TAMAAS__ - __END_TAMAAS__ - -#endif // KELVIN_HH +#endif // KELVIN_HH diff --git a/src/model/mindlin.cpp b/src/model/mindlin.cpp index 8d1bf65..dd72841 100644 --- a/src/model/mindlin.cpp +++ b/src/model/mindlin.cpp @@ -1,172 +1,185 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mindlin.hh" #include "elasto_plastic_model.hh" #include "influence.hh" #include "kelvin_helper.hh" /* -------------------------------------------------------------------------- */ -__BEGIN_TAMAAS__ +namespace tamaas { /* -------------------------------------------------------------------------- */ namespace detail { template class MindlinBoussinesqHelper { using boussinesq_t = T; using trait = model_type_traits; static constexpr UInt dim = trait::dimension; static constexpr UInt bdim = trait::boundary_dimension; public: MindlinBoussinesqHelper(const T& boussinesq, const influence::ElasticHelper& el, Real y_3, Real cutoff) : boussinesq(boussinesq), el(el), y_3(y_3), cutoff(cutoff) {} inline void applyIntegral(GridHermitian& out, GridHermitian& surface_gradients, - GridHermitian& source_stresses, + GridHermitian& source_stresses, const Grid& wavevectors) const { Loop::stridedLoop( [&](typename KelvinTrait::out_t&& out_local, MatrixProxy&& surface_gradient, MatrixProxy&& source_stress, VectorProxy&& q) { // Cutoff if (-q.l2norm() * std::abs(y_3) < std::log(cutoff)) return; - Vector normal{{{0, 0, 1}}}; - auto traction = el(surface_gradient) * - normal; // computing strains from gradu - traction -= source_stress * normal; // subtracting source stresses - out_local += boussinesq.applyU0(traction, q) * - influence::KelvinIntegrator<0>::g0(q.l2norm() * y_3); - out_local += (boussinesq.applyU1(traction, q) * - influence::KelvinIntegrator<0>::g1(q.l2norm() * y_3)); + Vector normal{{{0, 0, 1}}}; + auto traction = + el(surface_gradient) * normal; // computing strains from gradu + traction -= source_stress * normal; // subtracting source stresses + out_local += + boussinesq.applyU0(traction, q) * + influence::KelvinIntegrator<0>::g0(q.l2norm() * y_3); + out_local += + (boussinesq.applyU1(traction, q) * + influence::KelvinIntegrator<0>::g1(q.l2norm() * y_3)); }, out, surface_gradients, source_stresses, wavevectors); } inline void applyConstantTerm(typename KelvinTrait::out_t& out, Vector& source, const influence::ElasticHelper& el) const { out = 0; } protected: const T& boussinesq; const influence::ElasticHelper& el; const Real y_3, cutoff; }; template <> inline void MindlinBoussinesqHelper>:: applyConstantTerm( typename KelvinTrait>::out_t& out, Vector& source, const influence::ElasticHelper& el) const { out = 0; out(2, 0) = -source(0) / el.mu; out(2, 1) = -source(1) / el.mu; out(2, 2) = -source(2) / (el.lambda + 2 * el.mu); } } // namespace detail template -void Mindlin::apply(GridBase& source, - GridBase& out) const { +void Mindlin::applyIf(GridBase& source, GridBase& out, + filter_t pred) const { Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus(); constexpr UInt derivative = order - 2; influence::Kelvin kelvin(mu, nu); influence::Kelvin kelvin_strain(mu, nu); influence::Boussinesq boussinesq(mu, nu); influence::ElasticHelper elasticity(mu, nu); auto apply = [&](UInt node, decltype(this->source_buffers)& source_buffers, decltype(this->disp_buffer)& out_buffer) { const Real L = this->model->getSystemSize().front(); const UInt N = this->model->getDiscretization().front(); const Real dl = L / (N - 1); const Real xi = node * dl; // this should only be used when node == xi == 0 detail::LinearElement strain_element{ - thrust::make_pair(node, xi), kelvin_strain, source_buffers, surface_strains, - this->wavevectors}; + thrust::make_pair(node, xi), kelvin_strain, source_buffers, + surface_strains, this->wavevectors}; detail::LinearElement element{ thrust::make_pair(node, xi), kelvin, source_buffers, out_buffer, this->wavevectors}; out_buffer = 0; + // Simple condition to ensure correct integration + auto filter_out_elem = [&pred](UInt e) { + return not pred(e) and not pred(e + 1); + }; + /// Compute surface strains only once (dirty) if (node == 0) { surface_strains = 0; for (UInt e : Loop::range(N - 1)) { + if (filter_out_elem(e)) + continue; + const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl); strain_element.integrateElement(e, node_positions); } } // Apply Kelvin operator for (UInt e : Loop::range(N - 1)) { + if (filter_out_elem(e)) + continue; + const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl); element.integrateElement(e, node_positions); } // Correcting for the tractions on the surface detail::MindlinBoussinesqHelper helper( boussinesq, elasticity, xi, 1e-3); helper.applyIntegral(out_buffer, surface_strains, source_buffers.front(), this->wavevectors); typename detail::KelvinTrait::out_t out_fundamental( out_buffer(0)); - typename detail::KelvinTrait::source_t - source_fundamental(source_buffers[node](0)); + typename detail::KelvinTrait::source_t source_fundamental( + source_buffers[node](0)); Vector e3; e3 = 0; e3(trait::dimension - 1) = -1; auto in_fundamental = source_fundamental * e3; // Uniform shift in case of gradient computation helper.applyConstantTerm(out_fundamental, in_fundamental, elasticity); }; - this->fourierApply(apply, source, out); + this->fourierApplyIf(apply, source, out, pred); } /* -------------------------------------------------------------------------- */ /* Template instanciation */ /* -------------------------------------------------------------------------- */ template class Mindlin; template class Mindlin; /* -------------------------------------------------------------------------- */ -__END_TAMAAS__ +} // namespace tamaas diff --git a/src/model/mindlin.hh b/src/model/mindlin.hh index 0cdc323..513f8f8 100644 --- a/src/model/mindlin.hh +++ b/src/model/mindlin.hh @@ -1,64 +1,66 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef MINDLIN_HH #define MINDLIN_HH /* -------------------------------------------------------------------------- */ -#include "volume_potential.hh" -#include "model_type.hh" #include "grid_hermitian.hh" #include "kelvin.hh" +#include "model_type.hh" +#include "volume_potential.hh" /* -------------------------------------------------------------------------- */ -__BEGIN_TAMAAS__ +namespace tamaas { /** * @brief Mindlin tensor */ template class Mindlin : public Kelvin { static_assert(type == model_type::volume_1d || type == model_type::volume_2d, "Only volume types are supported"); using trait = model_type_traits; using parent = Kelvin; + using filter_t = typename parent::filter_t; public: /// Constructor Mindlin(Model* model) : parent(model) { surface_strains.setNbComponents(trait::components * trait::components); surface_strains.resize(this->disp_buffer.sizes()); } public: /// Apply the Mindlin-tensor_order operator - void apply(GridBase& source, GridBase& out) const override; + void applyIf(GridBase& source, GridBase& out, + filter_t pred) const override; protected: mutable GridHermitian surface_strains; }; - __END_TAMAAS__ +} // namespace tamaas -#endif // MINDLIN_HH +#endif // MINDLIN_HH diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh index 3035585..b94d685 100644 --- a/src/model/volume_potential.hh +++ b/src/model/volume_potential.hh @@ -1,105 +1,126 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef VOLUME_POTENTIAL_HH #define VOLUME_POTENTIAL_HH /* -------------------------------------------------------------------------- */ +#include "fft_plan_manager.hh" +#include "grid_hermitian.hh" +#include "grid_view.hh" #include "integral_operator.hh" #include "model_type.hh" -#include "grid_view.hh" -#include "grid_hermitian.hh" -#include "fft_plan_manager.hh" +/* -------------------------------------------------------------------------- */ +#include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Volume potential operator class. Applies the operators for computation * of displacements and strains due to residual/eigen strains */ template class VolumePotential : public IntegralOperator { using trait = model_type_traits; +protected: + using filter_t = std::function; + public: - VolumePotential(Model * model); + VolumePotential(Model* model); /// Update from model (does nothing) void updateFromModel() override {} /// Kind IntegralOperator::kind getKind() const override { return IntegralOperator::neumann; } /// Type model_type getType() const override; + /// Applying to all source + void apply(GridBase& input, GridBase& output) const override { + applyIf(input, output, [](UInt) { return true; }); + } + protected: /// Function to handle layer-by-layer Fourier treatment of operators template - void fourierApply(Func func, GridBase& in, GridBase& out) const; + void fourierApply(Func func, GridBase& in, GridBase& out) const { + fourierApplyIf(func, in, out, [](UInt) { return true; }); + } + + /// Function to handle layer-by-layer Fourier treatment of operators + template + void fourierApplyIf(Func func, GridBase& in, GridBase& out, + filter_t pred) const; protected: Grid wavevectors; using BufferType = GridHermitian; mutable std::vector source_buffers; mutable BufferType disp_buffer; }; /* -------------------------------------------------------------------------- */ /* Template implementation */ /* -------------------------------------------------------------------------- */ template template -void VolumePotential::fourierApply(Func func, GridBase& in, - GridBase& out) const { +void VolumePotential::fourierApplyIf(Func func, GridBase& in, + GridBase& out, + filter_t pred) const { constexpr UInt dim = trait::dimension; Grid& i = dynamic_cast(in); Grid& o = dynamic_cast(out); // TAMAAS_ASSERT(i.sizes().front() == o.sizes().front(), // "Number of layers does not match"); // ^^^^ not applicable for Boussinesq for (UInt layer : Loop::range(i.sizes().front())) { + if (not pred(layer)) + continue; + auto in_layer = make_view(i, layer); FFTPlanManager::get() .createPlan(in_layer, source_buffers[layer]) .forwardTransform(); } for (UInt layer : Loop::range(o.sizes().front())) { func(layer, source_buffers, disp_buffer); auto out_layer = make_view(o, layer); FFTPlanManager::get() .createPlan(out_layer, disp_buffer) .backwardTransform(); } } __END_TAMAAS__ #endif // VOLUME_POTENTIAL_HH