Page MenuHomec4science
Files F60250937

residual.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Apr 28, 15:44

residual.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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>
/* -------------------------------------------------------------------------- */
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);
// Initializing grids
for (auto grid : {&strain, &stress, &residual, &tmp}) {
*grid = allocateGrid<type, false, Real>(model->getDiscretization(),
trait::voigt);
}
// Registering fields
model->registerField("stress", stress);
model->registerField("strain", strain);
model->registerField("residual", residual);
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.getInternalData());
hardening.template computePlasticIncrement<false>(*residual, *strain, inc);
this->updateFilter(*residual);
this->model->applyElasticity(*residual, *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.template computeStress<false>(*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.template computeStress<true>(*stress, *strain, inc);
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.template computePlasticIncrement<false>(*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