Page MenuHomec4science
Files F69942022

residual.cpp
No OneTemporary
Actions

File Metadata

Created
Thu, Jul 4, 08:05

residual.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "residual.hh"
#include "grid_view.hh"
#include "model_factory.hh"
#include "model_type.hh"
#include <list>
#include <memory>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <model_type type>
ResidualTemplate<type>::ResidualTemplate(Model* model, Real sigma_0, Real h)
: Residual(model), hardening(model, sigma_0, h) {
// Registering operators for residual and displacement calculation
ModelFactory::registerVolumeOperators(*model);
const auto n = model->getDiscretization();
const auto nc = trait::voigt;
strain = model->request<type, false, Real>("strain", n, nc);
stress = model->request<type, false, Real>("stress", n, nc);
residual = model->request<type, false, Real>("residual", n, nc);
tmp = allocateGrid<type, false, Real>(n, nc);
plastic_filter = [this](UInt layer) {
#if 0
return true;
#else
return this->plastic_layers.find(layer) != this->plastic_layers.end();
#endif
};
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::computeResidual(GridBase<Real>& strain_increment) {
Grid<Real, dim> inc(strain->sizes(), strain->getNbComponents(),
strain_increment.view());
hardening.computeEigenStress(*residual, *strain, inc);
this->updateFilter(*residual);
integralOperator("mindlin_gradient")
.applyIf(*residual, *residual, plastic_filter);
integralOperator("boussinesq_gradient")
.apply(this->model->getTraction(), *tmp);
*residual -= strain_increment;
*residual += *tmp;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::applyTangent(
GridBase<Real>& output, GridBase<Real>& input,
GridBase<Real>& current_strain_increment) {
auto& inc = dynamic_cast<Grid<Real, dim>&>(current_strain_increment);
auto& out = dynamic_cast<Grid<Real, dim>&>(output);
auto& in = dynamic_cast<Grid<Real, dim>&>(input);
hardening.applyTangentIncrement(out, in, *strain, inc);
this->updateFilter(out);
integralOperator("mindlin_gradient").applyIf(out, out, plastic_filter);
out -= in;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::computeStress(GridBase<Real>& strain_increment) {
auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
hardening.computeStress(*stress, *strain, inc);
// integralOperator("boussinesq_gradient")
// .apply(this->model->getTraction(), tmp);
// this->model->applyElasticity(tmp, tmp);
// stress += tmp;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::updateState(
GridBase<Real>& converged_strain_increment) {
auto& inc = dynamic_cast<Grid<Real, dim>&>(converged_strain_increment);
hardening.computeStress(*stress, *strain, inc);
hardening.update();
this->model->applyElasticity(*residual, hardening.getPlasticStrain());
updateFilter(*residual);
*strain += converged_strain_increment;
// Computing displacements
integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
plastic_filter);
Grid<Real, dim> disp_tmp(this->model->getDiscretization(), trait::components);
integralOperator("boussinesq").apply(this->model->getTraction(), disp_tmp);
this->model->getDisplacement() += disp_tmp;
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::computeResidualDisplacement(
GridBase<Real>& strain_increment) {
auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
hardening.computeInelasticDeformationIncrement(*tmp, *strain, inc);
updateFilter(*tmp);
this->model->applyElasticity(*residual, *tmp);
integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
plastic_filter);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::updateFilter(
Grid<Real, dim>& plastic_strain_increment) {
// plastic_layers.clear();
for (UInt i : Loop::range(plastic_strain_increment.sizes().front())) {
auto slice = make_view(plastic_strain_increment, i);
if (slice.dot(slice) / slice.globalDataSize() > 1e-14)
plastic_layers.insert(i);
}
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::setIntegrationMethod(integration_method method,
Real cutoff) {
auto& mindlin = dynamic_cast<Mindlin<type, 1>&>(
*this->model->getIntegralOperator("mindlin"));
auto& mindlin_g = dynamic_cast<Mindlin<type, 2>&>(
*this->model->getIntegralOperator("mindlin_gradient"));
mindlin.setIntegrationMethod(method, cutoff);
mindlin_g.setIntegrationMethod(method, cutoff);
}
/* -------------------------------------------------------------------------- */
template class ResidualTemplate<model_type::volume_2d>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas

Event Timeline

c4science · Help