Page MenuHomec4science
Files F60244331

bem_grid_condat.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, Apr 28, 14:32

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), q_N(tmp, 3);
q_N = q;
q_N *= q_step;
//#pragma omp parallel for
for (auto it = gaps.begin(3); it < gaps.end(3); ++it) {
VectorProxy<Real> g(&(*it), 3);
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__

Event Timeline

c4science · Help