Page MenuHomec4science
Files F61370530

residual.cpp
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 06:12

residual.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 <list>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <UInt dim, UInt comp = dim>
void symmetrize(Grid<Real, dim>& /*grad*/) {
#if 0
Loop::loop(
[](MatrixProxy<Real, comp, comp> gradient) {
// 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));
}
}
},
range<MatrixProxy<Real, comp, comp>>(grad));
#endif
}
/* -------------------------------------------------------------------------- */
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->setNbComponents(trait::voigt);
grid->resize(model->getDiscretization());
}
// Registering fields
model->registerField("stress", getStress());
model->registerField("plastic_strain", getPlasticStrain());
model->registerField("residual", getVector());
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 = dynamic_cast<decltype(inc)>(strain_increment);
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;
symmetrize(residual);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::applyTangent(
GridBase<Real>& output, GridBase<Real>& input,
GridBase<Real>& current_strain_increment) {
Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(current_strain_increment);
Grid<Real, dim>& out = dynamic_cast<decltype(out)>(output);
Grid<Real, dim>& in = dynamic_cast<decltype(in)>(input);
hardening.applyTangentIncrement(out, in, strain, inc);
this->updateFilter(out);
integralOperator("mindlin_gradient").applyIf(out, out, plastic_filter);
out -= in;
symmetrize(residual);
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ResidualTemplate<type>::computeStress(GridBase<Real>& strain_increment) {
Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(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) {
Grid<Real, dim>& inc =
dynamic_cast<decltype(inc)>(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) {
Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(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.dataSize() > 1e-14)
plastic_layers.insert(i);
}
}
/* -------------------------------------------------------------------------- */
template class ResidualTemplate<model_type::volume_2d>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas

Event Timeline

c4science · Help