Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/solvers/kato.hh
9fe244aa2311features/declarative

Last Change
Raw File

kato.hh
View Options

/**
* @file
* LICENSE
*
* Copyright (©) 2016-2021 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef KATO_HH
#define KATO_HH
/* -------------------------------------------------------------------------- */
#include "contact_solver.hh"
#include "meta_functional.hh"
#include "model_type.hh"
#include "static_types.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
class Kato : public ContactSolver {
public:
/// Constructor
Kato(Model& model, const GridBase<Real>& surface, Real tolerance, Real mu);
public:
/// Solve
Real solve(GridBase<Real>& p0, UInt proj_iter);
/// Solve relaxed problem
Real solveRelaxed(GridBase<Real>& g0);
/// Solve regularized problem
Real solveRegularized(GridBase<Real>& p0, Real r);
/// Compute cost function
Real computeCost(bool use_tresca = false);
public:
/// Template for solve function
template <model_type type>
Real solveTmpl(GridBase<Real>& p0, UInt proj_iter);
/// Template for solveRelaxed function
template <model_type type>
Real solveRelaxedTmpl(GridBase<Real>& g0);
/// Template for solveRegularized function
template <model_type type>
Real solveRegularizedTmpl(GridBase<Real>& p0, Real r);
public:
/// Creates surfaceComp form surface
template <model_type type>
void initSurfaceWithComponents();
/// Compute gradient of functional
template <UInt comp>
void computeGradient(bool use_tresca = false);
/// Project pressure on friction cone
template <UInt comp>
void enforcePressureConstraints(GridBase<Real>& p0, UInt proj_iter);
/// Project on C
template <UInt comp>
void enforcePressureMean(GridBase<Real>& p0);
/// Project on D
template <UInt comp>
void enforcePressureCoulomb();
/// Project on D (Tresca)
template <UInt comp>
void enforcePressureTresca();
/// Comupte mean value of field
template <UInt comp>
Vector<Real, comp> computeMean(GridBase<Real>& field);
/// Add vector to each point of field
template <UInt comp>
void addUniform(GridBase<Real>& field, GridBase<Real>& vec);
/// Regularization function with factor r (0 -> unregugularized)
static __device__ __host__ Real regularize(Real x, Real r);
/// Compute grids of dual and primal variables
template <model_type type>
void computeValuesForCost(GridBase<Real>& lambda, GridBase<Real>& eta,
GridBase<Real>& p_N, GridBase<Real>& p_C);
/// Compute dual and primal variables with Tresca friction
template <model_type type>
void computeValuesForCostTresca(GridBase<Real>& lambda, GridBase<Real>& eta,
GridBase<Real>& p_N, GridBase<Real>& p_C);
/// Compute total displacement
template <UInt comp>
void computeFinalGap();
protected:
BEEngine& engine;
GridBase<Real>* gap = nullptr;
GridBase<Real>* pressure = nullptr;
std::unique_ptr<GridBase<Real>> surfaceComp = nullptr;
Real mu = 0;
UInt N = 0; // number of points
};
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::computeGradient(bool use_tresca) {
engine.solveNeumann(*pressure, *gap);
*gap -= *surfaceComp;
// Impose zero tangential displacement in non-sliding zone
UInt count_static = 0;
Vector<Real, comp> g_static;
using pvector = VectorProxy<Real, comp>;
const auto mu = this->mu;
if (!use_tresca) { // with Coulomb friction
count_static = Loop::reduce<operation::plus>(
[mu] CUDA_LAMBDA(VectorProxy<Real, comp> p) -> UInt {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact
return 1;
} else {
return 0;
}
},
range<pvector>(*pressure));
g_static = Loop::reduce<operation::plus>(
[mu] CUDA_LAMBDA(VectorProxy<Real, comp> g,
VectorProxy<Real, comp> p) -> Vector<Real, comp> {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu * p_N > p_T_norm) { // non-sliding contact
return g; // to compute mean of g_T
} else {
return 0;
}
},
range<pvector>(*gap), range<pvector>(*pressure));
} else { // with Tresca friction
count_static = Loop::reduce<operation::plus>(
[mu] CUDA_LAMBDA(VectorProxy<Real, comp> p) -> UInt {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_T_norm = p_T.l2norm();
Real p_N = p(comp - 1);
if (0.99 * mu > p_T_norm && p_N > 0) { // non-sliding contact
return 1;
} else {
return 0;
}
},
range<pvector>(*pressure));
g_static = Loop::reduce<operation::plus>(
[mu] CUDA_LAMBDA(VectorProxy<Real, comp> g,
VectorProxy<Real, comp> p) -> Vector<Real, comp> {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_norm = p_T.l2norm();
if (0.99 * mu > p_T_norm && p_N > 0) { // non-sliding contact
return g; // to compute mean of g_T
} else {
return 0;
}
},
range<pvector>(*gap), range<pvector>(*pressure));
}
if (count_static != 0) {
g_static /= count_static;
// g_static(comp - 1) = 0;
} else { // if no sticking zone, mean computed on sliding zone
count_static = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> p) -> UInt {
if (p(comp - 1) > 0) {
return 1;
} else {
return 0;
}
},
range<pvector>(*pressure));
Real g_N_static = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> g,
VectorProxy<Real, comp> p) -> Real {
if (p(comp - 1) > 0) {
return g(comp - 1);
} else {
return 0;
}
},
range<pvector>(*gap), range<pvector>(*pressure));
g_static(comp - 1) = g_N_static / count_static;
}
// Add frictionnal term to functional
if (!use_tresca) { // with Coulomb friction
Loop::loop(
[mu, g_static] CUDA_LAMBDA(VectorProxy<Real, comp> g) {
g -= g_static;
VectorProxy<Real, comp - 1> g_T(g(0));
Real g_T_norm = g_T.l2norm();
g(comp - 1) += mu * g_T_norm; // Frictionnal work
},
range<pvector>(*gap));
} else { // with Tresca friction
*gap -= g_static;
// Loop::stridedLoop(
// [this, g_static] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
// VectorProxy<Real, comp>&& p) {
// g -= g_static;
// VectorProxy<Real, comp - 1> g_T = g(0);
// Real g_T_norm = g_T.l2norm();
// Real p_N = p(comp - 1);
// if (p_N > 0)
// g(comp - 1) += mu/p_N * g_T_norm; // Frictionnal work
// },
// *gap, *pressure);
}
// Loop::stridedLoop(
// [this] CUDA_LAMBDA(VectorProxy<Real, comp>&& g, VectorProxy<Real, comp>&&
// p) {
// VectorProxy<Real, comp - 1> g_T = g(0);
// Real g_T_norm = g_T.l2norm();
// VectorProxy<Real, 1> g_N = g(comp - 1);
// VectorProxy<Real, comp - 1> p_T = p(0);
// Real p_T_norm = p_T.l2norm();
// VectorProxy<Real, 1> p_N = p(comp - 1);
// if (p_T_norm != 0) {
// g_T = p_T;
// g_T *= -g_T_norm / p_T_norm;
// }
// g_N += mu * g_T_norm; // Frictionnal work
// },
// *gap, *pressure);
}
/* -------------------------------------------------------------------------- */
/**
* Projects \f$\vec{p}\f$ on \f$\mathcal{C}\f$ and \f$\mathcal{D}\f$.
*/
// template <UInt comp>
// void Kato::enforcePressureConstraints(GridBase<Real>& p0, UInt proj_iter) {
// for (UInt i = 0; i < proj_iter; i++) {
// enforcePressureMean<comp>(p0);
// enforcePressureCoulomb<comp>();
// }
// }
/* -------------------------------------------------------------------------- */
/**
* Projects \f$\vec{p}\f$ on \f$\mathcal{C}\f$ and \f$\mathcal{D}\f$.
*/
template <UInt comp>
void Kato::enforcePressureConstraints(GridBase<Real>& p0, UInt proj_iter) {
for (UInt i = 0; i < proj_iter; i++) {
enforcePressureMean<comp>(p0);
enforcePressureCoulomb<comp>();
}
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::enforcePressureCoulomb() {
const auto mu = this->mu;
Loop::loop(
[mu] CUDA_LAMBDA(VectorProxy<Real, comp> p) {
VectorProxy<Real, comp - 1> p_T(p(0));
Real p_N = p(comp - 1);
Real p_T_sqrd = p_T.l2squared();
// Projection normale au cône de friction
bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N);
bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu);
if (cond2) {
p_T = 0;
p(comp - 1) = 0;
} else if (!cond1) {
Real p_T_norm = std::sqrt(p_T_sqrd);
Real k = (p_N + mu * p_T_norm) / (1 + mu * mu);
p_T *= k * mu / p_T_norm;
p(comp - 1) = k;
}
},
range<VectorProxy<Real, comp>>(*pressure));
}
/* -------------------------------------------------------------------------- */
/**
* Compute mean of the field taking each component separately.
*/
template <UInt comp>
Vector<Real, comp> Kato::computeMean(GridBase<Real>& field) {
Vector<Real, comp> mean = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(VectorProxy<Real, comp> f) -> Vector<Real, comp> {
return f;
},
range<VectorProxy<Real, comp>>(field));
mean /= N;
return mean;
}
/* -------------------------------------------------------------------------- */
template <UInt comp>
void Kato::addUniform(GridBase<Real>& field, GridBase<Real>& vec) {
VectorProxy<Real, comp> _vec(vec(0));
field += _vec;
}
} // namespace tamaas
#endif // KATO_HH
c4science · Help