Page MenuHomec4science
Files F65647172

residual.cpp
No OneTemporary
Actions

File Metadata

Created
Wed, Jun 5, 06:39

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 {
/* -------------------------------------------------------------------------- */
Residual::Residual(Model& model, std::shared_ptr<Material> material)
: model(model), material(material) {
// Registering operators for residual and displacement calculation
ModelFactory::registerVolumeOperators(model);
const auto n = model.getDiscretization();
const auto nc = model_type_traits<type>::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
};
}
/* -------------------------------------------------------------------------- */
void Residual::computeResidual(GridBase<Real>& strain_increment) {
// Grid<Real, dim> inc{stress->sizes(), stress->getNbComponents(),
// strain_increment.view()};
material->computeEigenStress(*residual, *strain, strain_increment);
this->updateFilter(*residual);
integralOperator("mindlin_gradient")
.applyIf(*residual, *residual, plastic_filter);
integralOperator("boussinesq_gradient").apply(model.getTraction(), *tmp);
*residual -= strain_increment;
*residual += *tmp;
}
/* -------------------------------------------------------------------------- */
void Residual::computeResidualDisplacement(GridBase<Real>& strain_increment) {
material->computeEigenStress(*tmp, *strain, strain_increment);
updateFilter(*tmp);
integralOperator("mindlin").applyIf(*tmp, model.getDisplacement(),
plastic_filter);
}
/* -------------------------------------------------------------------------- */
void Residual::updateState(GridBase<Real>& converged_strain_increment) {
material->computeStress(*stress, *strain, converged_strain_increment);
material->computeEigenStress(*residual, *strain, converged_strain_increment);
material->update();
updateFilter(*residual);
*strain += converged_strain_increment;
// Computing displacements
integralOperator("mindlin").applyIf(*residual, model.getDisplacement(),
plastic_filter);
Grid<Real, dim> disp_tmp(model.getDiscretization(),
model_type_traits<type>::components);
integralOperator("boussinesq").apply(model.getTraction(), disp_tmp);
model.getDisplacement() += disp_tmp;
}
/* -------------------------------------------------------------------------- */
void Residual::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);
}
}
} // namespace tamaas

Event Timeline

c4science · Help