Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/solvers/epic.cpp
a240eb50bf7efeatures/netcdf-mpi

Last Change
Raw File

epic.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 "epic.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
EPICSolver::EPICSolver(ContactSolver& csolver, EPSolver& epsolver,
Real tolerance, Real relaxation)
: csolver(csolver), epsolver(epsolver), tolerance(tolerance),
relaxation(relaxation) {
auto& model = csolver.getModel();
surface.wrap(csolver.getSurface());
pressure.resize(model.getTraction().getNbPoints());
#define SET_VIEWS(_, __, type) \
case type: { \
setViews<type>(); \
break; \
}
switch (model.getType()) {
BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
}
#undef SET_VIEWS
}
/* -------------------------------------------------------------------------- */
void EPICSolver::solve(std::vector<Real> load) {
UInt n = 0, nmax = 1000;
Real error = 0.;
const GridBase<Real> initial_surface(surface);
Real normalizing_factor = std::sqrt(initial_surface.var());
GridBase<Real> previous_residual(*residual_disp);
GridBase<Real> relaxation_direction(*residual_disp);
previous_residual = 0;
do {
csolver.solve(load); // solve contact problem
*pressure_inc -= pressure; // compute pressure increment
epsolver.solve(); // compute residual displacement
relaxation_direction = *residual_disp;
relaxation_direction -= previous_residual;
relaxation_direction *= relaxation;
error = computeError(*residual_disp, previous_residual, normalizing_factor);
previous_residual += relaxation_direction;
surface = initial_surface;
surface -= previous_residual; // update surface
std::cout << "error = " << std::scientific << error << std::fixed
<< std::endl;
} while (error > tolerance && n++ < nmax);
// computing full pressure
pressure += *pressure_inc;
// setting the model pressure to full converged pressure
*pressure_inc = pressure;
// updating plastic state
epsolver.updateState();
}
/* -------------------------------------------------------------------------- */
Real EPICSolver::computeError(const GridBase<Real>& current,
const GridBase<Real>& prev, Real factor) const {
GridBase<Real> error(current);
error -= prev;
return std::sqrt(error.dot(error)) / factor;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
c4science · Help