Page MenuHomec4science
Files F62839547

residual.cpp
No OneTemporary
Actions

File Metadata

Created
Thu, May 16, 00:30

residual.cpp
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "residual.hh"
#include "grid_view.hh"
#include <list>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <UInt dim, UInt comp = dim>
void symmetrize(Grid<Real, dim>& grad) {
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));
}
/* -------------------------------------------------------------------------- */
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 calculation
model->registerIntegralOperator<Mindlin<type, 2>>("mindlin_gradient");
model->registerIntegralOperator<Boussinesq<type, 1>>("boussinesq_gradient");
// Registering operator for displacement calculation
model->registerIntegralOperator<Mindlin<type, 1>>("mindlin");
model->registerIntegralOperator<Boussinesq<type, 0>>("boussinesq");
// Initializing grids
for (auto grid : {&strain, &stress, &residual, &tmp}) {
grid->setNbComponents(trait::components * trait::components);
grid->resize(model->getDiscretization());
}
// Registering fields
model->registerField("stress", getStress());
model->registerField("plastic_strain", getPlasticStrain());
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