Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/bem_grid_condat.cpp
0e70f46e9d04maxwell_viscoelastic…

Last Change
Raw File

bem_grid_condat.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 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 "bem_grid_condat.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BemGridCondat::BemGridCondat(Surface<Real> & surface):
BemGridKato(surface),
grad_step(1.8),
q_step(0.1)
{}
/* -------------------------------------------------------------------------- */
BemGridCondat::~BemGridCondat() {}
/* -------------------------------------------------------------------------- */
void BemGridCondat::updateGradient(Grid<Real, 1> &q_vec) {
Real tmp[3] = {0};
VectorProxy<Real>
q(q_vec.getInternalData(), 3),
g(nullptr, 3),
q_N(tmp, 3);
const UInt N = gaps.getNbPoints();
q_N = q;
q_N *= q_step;
#pragma omp parallel for firstprivate(g)
for (UInt i = 0 ; i < N ; i++) {
g.setPointer(gaps.getInternalData() + i * g.dataSize());
g *= grad_step / lipschitz;
g += q;
}
}
/* -------------------------------------------------------------------------- */
void BemGridCondat::computeEquilibrium(Real tol, const Grid<Real, 1> & p_target) {
this->computeInfluence();
Real f = 0;
UInt n_iter = 0;
TAMAAS_ASSERT(p_target.dataSize() == 3, "Applied pressure must have only 3 components");
TAMAAS_ASSERT(q_step > 0 && q_step < 1, "q step should be in ]0, 1[");
TAMAAS_ASSERT(grad_step * lipschitz < 2 * (1 - q_step), "grad step * lipschitz should be less than 2 q steps");
this->tractions.uniformSetComponents(p_target);
Real q_mem[3] = {0};
Real e_mem[3] = {0};
Real pn_mem[3] = {0};
Real pn1_mem[3] = {0};
Real zero_mem[3] = {0};
VectorProxy<Real> q(q_mem, 3), e(e_mem, 3), pn(pn_mem, 3), pn1(pn1_mem, 3), zero(zero_mem, 3);
Real norm =
p_target(0)*p_target(0)
+ p_target(1)*p_target(1)
+ p_target(2)*p_target(2);
do {
computeDisplacementsFromTractions();
computeGaps(displacements);
updateGradient(q);
beta = computeLambda();
BemGrid::computeMeanValueError(tractions, zero, pn);
tractions -= gaps;
BemGrid::projectOnFrictionCone(tractions, mu);
BemGrid::computeMeanValueError(tractions, zero, pn1);
pn1 *= 2;
q -= p_target;
q += pn1;
q -= pn;
lambda -= beta;
f = computeStoppingCriterion();
printFriendlyMessage(n_iter, f);
} while ((f > tol || e.l2norm() / std::sqrt(norm) > tol) && n_iter++ < this->max_iter);
}
__END_TAMAAS__
c4science · Help