diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 6cbd075..39eb938 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,246 +1,252 @@
/*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
+#include "anderson.hh"
#include "beck_teboulle.hh"
#include "condat.hh"
#include "contact_solver.hh"
#include "dfsane_solver.hh"
#include "ep_solver.hh"
#include "epic.hh"
#include "kato.hh"
#include "kato_saturated.hh"
#include "polonsky_keer_rey.hh"
#include "polonsky_keer_tan.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
/* -------------------------------------------------------------------------- */
using namespace py::literals;
class PyEPSolver : public EPSolver {
public:
using EPSolver::EPSolver;
// NOLINTNEXTLINE(readability-else-after-return)
void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
void updateState() override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD(void, EPSolver, updateState);
}
};
class PyContactSolver : public ContactSolver {
public:
using ContactSolver::ContactSolver;
Real solve(std::vector load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
Real solve(Real load) override {
PYBIND11_OVERLOAD(Real, ContactSolver, solve, load);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_(mod, "ContactSolver")
.def(py::init&, Real>(),
py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def(
"setMaxIterations",
[](ContactSolver& m, UInt iter) {
TAMAAS_DEPRECATE("setMaxIterations()", "the max_iter property");
m.setMaxIterations(iter);
},
"max_iter"_a)
.def(
"setDumpFrequency",
[](ContactSolver& m, UInt freq) {
TAMAAS_DEPRECATE("setDumpFrequency()", "the dump_freq property");
m.setDumpFrequency(freq);
},
"dump_freq"_a)
.def_property("tolerance", &ContactSolver::getTolerance,
&ContactSolver::setTolerance, "Solver tolerance")
.def_property("max_iter", &ContactSolver::getMaxIterations,
&ContactSolver::setMaxIterations,
"Maximum number of iterations")
.def_property("dump_freq", &ContactSolver::getDumpFrequency,
&ContactSolver::setDumpFrequency,
"Frequency of displaying solver info")
.def_property_readonly("model", &ContactSolver::getModel)
.def_property_readonly("surface", &ContactSolver::getSurface)
.def_property("functional", &ContactSolver::getFunctional,
&ContactSolver::setFunctional)
.def("addFunctionalTerm", &ContactSolver::addFunctionalTerm,
"Add a term to the contact functional to minimize")
.def("solve", py::overload_cast>(&ContactSolver::solve),
"target_force"_a,
py::call_guard(),
"Solve the contact for a mean traction/gap vector")
.def("solve", py::overload_cast(&ContactSolver::solve),
"target_normal_pressure"_a,
py::call_guard(),
"Solve the contact for a mean normal pressure/gap");
py::class_ pkr(
mod, "PolonskyKeerRey",
"Main solver class for normal elastic contact problems. Its functional "
"can be customized to add an adhesion term, and its primal variable can "
"be set to either the gap or the pressure.");
// Need to export enum values before defining PKR constructor
py::enum_(pkr, "type")
.value("gap", PolonskyKeerRey::gap)
.value("pressure", PolonskyKeerRey::pressure)
.export_values();
pkr.def(py::init&, Real, PolonskyKeerRey::type,
PolonskyKeerRey::type>(),
"model"_a, "surface"_a, "tolerance"_a,
"primal_type"_a = PolonskyKeerRey::type::pressure,
"constraint_type"_a = PolonskyKeerRey::type::pressure,
py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def("computeError", &PolonskyKeerRey::computeError)
.def("setIntegralOperator", &PolonskyKeerRey::setIntegralOperator);
py::class_(
mod, "KatoSaturated",
"Solver for pseudo-plasticity problems where the normal pressure is "
"constrained above by a saturation pressure \"pmax\"")
.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "pmax"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def_property("pmax", &KatoSaturated::getPMax, &KatoSaturated::setPMax,
"Saturation normal pressure");
py::class_ kato(mod, "Kato");
kato.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &Kato::solve, "p0"_a, "proj_iter"_a = 50)
.def("solveRelaxed", &Kato::solveRelaxed, "g0"_a)
.def("solveRegularized", &Kato::solveRegularized, "p0"_a, "r"_a = 0.01)
.def("computeCost", &Kato::computeCost, "use_tresca"_a = false);
py::class_ bt(mod, "BeckTeboulle");
bt.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &BeckTeboulle::solve, "p0"_a)
.def("computeCost", &BeckTeboulle::computeCost);
py::class_ cd(
mod, "Condat",
"Main solver for frictional contact problems. It has no restraint on the "
"material properties or friction coefficient values, but solves an "
"associated version of the Coulomb friction law, which differs from the "
"traditional Coulomb friction in that the normal and tangential slip "
"components are coupled.");
cd.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &Condat::solve, "p0"_a, "grad_step"_a = 0.9)
.def("computeCost", &Condat::computeCost);
py::class_ pkt(mod, "PolonskyKeerTan");
pkt.def(py::init&, Real, Real>(), "model"_a,
"surface"_a, "tolerance"_a, "mu"_a, py::keep_alive<1, 2>(),
py::keep_alive<1, 3>())
.def("solve", &PolonskyKeerTan::solve, "p0"_a)
.def("solveTresca", &PolonskyKeerTan::solveTresca, "p0"_a)
.def("computeCost", &PolonskyKeerTan::computeCost,
"use_tresca"_a = false);
py::class_(
mod, "_tolerance_manager",
"Manager object for the tolereance of nonlinear plasticity solvers. "
"Decreases the solver tolerance by geometric progression.")
.def(py::init(), "start_tol"_a, "end_tol"_a, "rate"_a);
py::class_(
mod, "EPSolver", "Mother class for nonlinear plasticity solvers")
.def(py::init(), "residual"_a, py::keep_alive<1, 2>())
.def("solve", &EPSolver::solve)
.def("getStrainIncrement", &EPSolver::getStrainIncrement,
py::return_value_policy::reference_internal)
.def("getResidual", &EPSolver::getResidual,
py::return_value_policy::reference_internal)
.def("updateState", &EPSolver::updateState)
.def_property("tolerance", &EPSolver::getTolerance,
&EPSolver::setTolerance)
.def("setToleranceManager", &EPSolver::setToleranceManager)
.def("beforeSolve", &EPSolver::beforeSolve);
py::class_(mod, "_DFSANESolver")
.def(py::init(), "residual"_a, py::keep_alive<1, 2>())
.def(py::init([](Residual& res, Model&) {
TAMAAS_DEPRECATE("Solver(residual, model)", "Solver(residual)");
return std::make_unique(res);
}),
"residual"_a, "model"_a, py::keep_alive<1, 2>());
py::class_(
mod, "EPICSolver",
"Main solver class for elastic-plastic contact problems")
.def(py::init(),
"contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10,
"relaxation"_a = 0.3, py::keep_alive<1, 2>(), py::keep_alive<1, 3>())
.def(
"solve",
[](EPICSolver& solver, Real pressure) {
return solver.solve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard(),
"Solves the EP contact with a relaxed fixed-point scheme. Adjust "
"the relaxation parameter to help convergence.")
.def(
"acceleratedSolve",
[](EPICSolver& solver, Real pressure) {
return solver.acceleratedSolve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard(),
"Solves the EP contact with an accelerated fixed-point scheme. May "
"not converge!")
.def_property("tolerance", &EPICSolver::getTolerance,
&EPICSolver::setTolerance)
.def_property("relaxation", &EPICSolver::getRelaxation,
&EPICSolver::setRelaxation)
.def_property_readonly("model", &EPICSolver::getModel);
+
+ py::class_(mod, "AndersonMixing")
+ .def(py::init(),
+ "contact_solver"_a, "elasto_plastic_solver"_a, "tolerance"_a = 1e-10,
+ "memory"_a = 5, py::keep_alive<1, 2>(), py::keep_alive<1, 3>());
}
} // namespace wrap
} // namespace tamaas
diff --git a/src/SConscript b/src/SConscript
index 2ed290c..fdd85ae 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,147 +1,148 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
# 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 .
from SCons.Script import Import, Export
def prepend(env, path, list):
return [env.File(x, path) for x in list]
Import('main_env')
env = main_env.Clone()
env.AddMethod(prepend, 'PrependDir')
# Core
core_list = """
fft_engine.cpp
grid.cpp
grid_hermitian.cpp
statistics.cpp
tamaas.cpp
loop.cpp
computes.cpp
logger.cpp
mpi_interface.cpp
""".split()
core_list = env.PrependDir('core', core_list)
if env['use_fftw']:
core_list += ['core/fftw/fftw_engine.cpp']
if env['use_mpi']:
core_list += ['core/fftw/mpi/fftw_mpi_engine.cpp']
if env['use_cuda']:
core_list += ['core/cuda/cufft_engine.cpp']
info_file = env.Substfile('tamaas_info.cpp', 'tamaas_info.cpp.in')
env.Depends(info_file, Value(env.subst(env['SUBST_DICT'].values())))
core_list.append(info_file)
# Lib roughcontact
generator_list = """
surface_generator.cpp
surface_generator_filter.cpp
surface_generator_random_phase.cpp
isopowerlaw.cpp
regularized_powerlaw.cpp
""".split()
generator_list = env.PrependDir('surface', generator_list)
# Lib PERCOLATION
percolation_list = """
flood_fill.cpp
""".split()
percolation_list = env.PrependDir('percolation', percolation_list)
# Model
model_list = """
model.cpp
model_factory.cpp
model_type.cpp
model_template.cpp
integral_operator.cpp
be_engine.cpp
westergaard.cpp
dcfft.cpp
elastic_functional.cpp
meta_functional.cpp
adhesion_functional.cpp
volume_potential.cpp
kelvin.cpp
mindlin.cpp
boussinesq.cpp
hooke.cpp
residual.cpp
materials/isotropic_hardening.cpp
materials/mfront_material.cpp
integration/element.cpp
""".split()
model_list = env.PrependDir('model', model_list)
# Solvers
solvers_list = """
contact_solver.cpp
polonsky_keer_rey.cpp
kato_saturated.cpp
kato.cpp
beck_teboulle.cpp
condat.cpp
polonsky_keer_tan.cpp
ep_solver.cpp
dfsane_solver.cpp
epic.cpp
+anderson.cpp
""".split()
solvers_list = env.PrependDir('solvers', solvers_list)
# Assembling total list
rough_contact_list = \
core_list + generator_list + percolation_list + model_list + solvers_list
# Adding extra warnings for Tamaas base lib
env.AppendUnique(CXXFLAGS=['-Wextra'])
# Allowing libTamaas.so to find libs in the same directory
env.AppendUnique(RPATH=["'$$$$ORIGIN'"])
# Build static library for packaging
if env['build_static_lib']:
env.AppendUnique(CXXFLAGS='-fPIC')
libTamaas = env.StaticLibrary('Tamaas', rough_contact_list)
# Build shared library (default)
else:
libTamaas = env.SharedLibrary('Tamaas', rough_contact_list,
SHLIBVERSION=env['version'])
# Specify install target to install lib
lib_prefix = env.Dir('lib', env['prefix'])
lib_install = env.InstallVersionedLib(target=lib_prefix,
source=libTamaas)
# Defining alias targets
main_env.Alias('build-cpp', libTamaas)
main_env.Alias('install-lib', lib_install)
# Export target for use in python builds
Export('libTamaas')
diff --git a/src/solvers/anderson.cpp b/src/solvers/anderson.cpp
new file mode 100644
index 0000000..529d1eb
--- /dev/null
+++ b/src/solvers/anderson.cpp
@@ -0,0 +1,192 @@
+/*
+ * 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 .
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "anderson.hh"
+#include
+#include
+/* -------------------------------------------------------------------------- */
+namespace tamaas {
+
+/* -------------------------------------------------------------------------- */
+AndersonMixing::AndersonMixing(ContactSolver& csolver, EPSolver& epsolver,
+ Real tolerance, UInt memory)
+ : EPICSolver(csolver, epsolver, tolerance, 0.5), M(memory) {}
+
+/* -------------------------------------------------------------------------- */
+Real AndersonMixing::solve(const std::vector& load) {
+ UInt n = 0, nmax = 1000;
+ Real error = 0.;
+ const GridBase initial_surface(surface);
+ Real normalizing_factor = std::sqrt(initial_surface.var());
+ GridBase x(*residual_disp), x_prev(*residual_disp);
+ memory_t memory, residual_memory;
+
+ do {
+ {
+ // Compute a residual vector for this iteration
+ GridBase F(*residual_disp);
+ fixedPoint(F, x, initial_surface, load);
+ F -= x; // F = g(x) - x
+ memory.push_front(x);
+ residual_memory.push_front(std::move(F));
+ }
+
+ // Keep M previous iterations: memory has M + 1 vectors
+ if (memory.size() > M + 1)
+ memory.pop_back();
+ if (residual_memory.size() > M + 1)
+ residual_memory.pop_back();
+
+ x = mixingUpdate(std::move(x), computeGamma(residual_memory), memory,
+ residual_memory, relaxation);
+ error = computeError(x, x_prev, normalizing_factor);
+ x_prev = x;
+
+ Logger().get(LogLevel::info) << "[Anderson] error = " << std::scientific
+ << error << std::fixed << std::endl;
+ } while (error > tolerance && n++ < nmax);
+
+ std::cout << "Converged in " << n << " iterations\n";
+
+ // computing full pressure
+ pressure += *pressure_inc;
+ // setting the model pressure to full converged pressure
+ *pressure_inc = pressure;
+ // updating plastic state
+ epsolver.updateState();
+ return error;
+}
+
+/* -------------------------------------------------------------------------- */
+GridBase AndersonMixing::mixingUpdate(GridBase x,
+ std::vector gamma,
+ const memory_t& memory,
+ const memory_t& residual_memory,
+ Real relaxation) {
+ Loop::loop([beta = relaxation](Real& x, Real Fl) { x += beta * Fl; }, x,
+ residual_memory.front());
+
+ for (UInt m = 1; m < memory.size(); ++m) {
+ Loop::loop(
+ [gm = gamma[m - 1], beta = relaxation](Real& x, Real xmp1, Real xm,
+ Real Fmp1, Real Fm) {
+ // Eyert (10.1006/jcph.1996.0059) p282 eq (7.7)
+ x -= gm * ((xmp1 - xm) + beta * (Fmp1 - Fm));
+ },
+ x, memory[m - 1], memory[m], residual_memory[m - 1],
+ residual_memory[m]);
+ }
+
+ return x;
+}
+
+/* -------------------------------------------------------------------------- */
+/// Crout(e) (Antoine... et Daniel...) factorization from
+/// https://en.wikipedia.org/wiki/Crout_matrix_decomposition
+auto factorize(Grid A) {
+ TAMAAS_ASSERT(A.sizes()[0] == A.sizes()[1], "Matrix is not square");
+ const auto N = A.sizes()[0];
+ auto LU =
+ std::make_pair(Grid{A.sizes(), 1}, Grid{A.sizes(), 1});
+ auto& L = LU.first;
+ auto& U = LU.second;
+
+ for (UInt i = 0; i < N; ++i) {
+ U(i, i) = 1;
+ }
+
+ for (UInt j = 0; j < N; ++j) {
+ for (UInt i = j; i < N; ++i) {
+ double sum = 0;
+ for (UInt k = 0; k < j; ++k) {
+ sum += L(i, k) * U(k, j);
+ }
+
+ L(i, j) = A(i, j) - sum;
+ }
+
+ for (UInt i = j; i < N; ++i) {
+ double sum = 0;
+ for (UInt k = 0; k < j; ++k) {
+ sum += L(j, k) * U(k, i);
+ }
+
+ U(j, i) = (A(j, i) - sum) / L(j, j);
+ }
+ }
+ return LU;
+}
+
+/* -------------------------------------------------------------------------- */
+auto LUsubstitute(std::pair, Grid> LU, Grid b) {
+ const auto& L = LU.first;
+ const auto& U = LU.second;
+ const auto N = b.sizes()[0];
+
+ // Forward substitution
+ for (UInt m = 0; m < N; ++m) {
+ for (UInt i = 0; i < m; ++i)
+ b(m) -= L(m, i) * b(i);
+ b(m) /= L(m, m);
+ }
+
+ // Backward substitution
+ for (UInt l = 0; l < N; ++l) {
+ auto m = N - 1 - l;
+ for (UInt i = m + 1; i < N; ++i)
+ b(m) -= U(m, i) * b(i);
+ b(m) /= U(m, m);
+ }
+
+ return b;
+}
+
+/* -------------------------------------------------------------------------- */
+std::vector AndersonMixing::computeGamma(const memory_t& residual) {
+ auto N = residual.size() - 1;
+ Grid A({N, N}, 1);
+ Grid b({N}, 1);
+ const auto& Fl = residual.front();
+
+ // Fill RHS vector
+ for (UInt n = 1; n < residual.size(); ++n) {
+ b(n - 1) = residual[n - 1].dot(Fl) - residual[n].dot(Fl);
+ }
+
+ /// Regularization from Eyert
+ auto w0 = 0;
+ // Fill matrix
+ for (UInt n = 1; n < residual.size(); ++n) {
+ for (UInt m = 1; m < residual.size(); ++m) {
+ A(n - 1, m - 1) = residual[n - 1].dot(residual[m - 1]) -
+ residual[n - 1].dot(residual[m]) -
+ residual[n].dot(residual[m - 1]) +
+ residual[n].dot(residual[m]);
+ A(n - 1, m - 1) *= (1 + w0 * w0 * (n == m));
+ }
+ }
+
+ auto x = LUsubstitute(factorize(A), b);
+ return std::vector(x.begin(), x.end());
+}
+
+} // namespace tamaas
diff --git a/src/solvers/anderson.hh b/src/solvers/anderson.hh
new file mode 100644
index 0000000..87f415d
--- /dev/null
+++ b/src/solvers/anderson.hh
@@ -0,0 +1,53 @@
+/*
+ * 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 .
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef ANDERSON_HH
+#define ANDERSON_HH
+/* -------------------------------------------------------------------------- */
+#include "epic.hh"
+#include
+/* -------------------------------------------------------------------------- */
+namespace tamaas {
+
+class AndersonMixing : public EPICSolver {
+ using memory_t = std::deque>;
+
+public:
+ AndersonMixing(ContactSolver& csolver, EPSolver& epsolver,
+ Real tolerance = 1e-10, UInt memory = 5);
+
+ Real solve(const std::vector& load) override;
+
+private:
+ static std::vector computeGamma(const memory_t& residual);
+ static GridBase mixingUpdate(GridBase x, std::vector gamma,
+ const memory_t& memory,
+ const memory_t& residual_memory,
+ Real relaxation);
+
+private:
+ UInt M;
+};
+
+} // namespace tamaas
+/* -------------------------------------------------------------------------- */
+#endif
diff --git a/src/solvers/epic.hh b/src/solvers/epic.hh
index 4b6f9c4..1f13d7e 100644
--- a/src/solvers/epic.hh
+++ b/src/solvers/epic.hh
@@ -1,88 +1,88 @@
/*
* 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef EPIC_HH
#define EPIC_HH
/* -------------------------------------------------------------------------- */
#include "contact_solver.hh"
#include "ep_solver.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
class EPICSolver {
public:
/// Constructor
EPICSolver(ContactSolver& csolver, EPSolver& epsolver, Real tolerance = 1e-10,
Real relaxation = 0.3);
- Real solve(const std::vector& load);
+ virtual Real solve(const std::vector& load);
Real acceleratedSolve(const std::vector& load);
Real computeError(const GridBase& current, const GridBase& prev,
Real factor) const;
void fixedPoint(GridBase& result, const GridBase& x,
const GridBase& initial_surface,
std::vector load);
template
void setViews();
TAMAAS_ACCESSOR(tolerance, Real, Tolerance);
TAMAAS_ACCESSOR(relaxation, Real, Relaxation);
const Model& getModel() const { return csolver.getModel(); }
protected:
GridBase surface; ///< corrected surface
GridBase pressure; ///< current pressure
std::unique_ptr> residual_disp; ///< plastic residual disp
std::unique_ptr> pressure_inc; ///< pressure increment
ContactSolver& csolver;
EPSolver& epsolver;
Real tolerance, relaxation;
};
/* -------------------------------------------------------------------------- */
/* Template impl. */
/* -------------------------------------------------------------------------- */
template
void EPICSolver::setViews() {
constexpr UInt dim = model_type_traits::dimension;
constexpr UInt bdim = model_type_traits::boundary_dimension;
constexpr UInt comp = model_type_traits::components;
pressure_inc =
std::unique_ptr>{new GridView(
csolver.getModel().getTraction(), {}, comp - 1)};
residual_disp =
std::unique_ptr>{new GridView(
csolver.getModel().getDisplacement(),
model_type_traits::indices, comp - 1)};
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // EPIC_HH
diff --git a/tests/test_epic.py b/tests/test_epic.py
index cb3f87f..ed24b09 100644
--- a/tests/test_epic.py
+++ b/tests/test_epic.py
@@ -1,78 +1,77 @@
# -*- coding: utf-8 -*-
#
# 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 .
from __future__ import division, print_function
import pytest
from conftest import pkridfn
import tamaas as tm
from numpy.linalg import norm
from tamaas.nonlinear_solvers import DFSANESolver
@pytest.fixture(scope="module",
- params=[tm.PolonskyKeerRey.pressure],
- ids=pkridfn)
+ params=[tm.EPICSolver, tm.AndersonMixing],
+ ids=["EPIC", "Anderson"])
def model_setup(hertz_coarse, request):
models, solvers = [], []
for mtype in [tm.model_type.basic_2d, tm.model_type.volume_2d]:
dim = tm.type_traits[mtype].dimension
n = [hertz_coarse.n] * dim
d = [hertz_coarse.domain_size] * dim
if dim == 3:
n[0] = 2
d[0] = 0.001
model = tm.ModelFactory.createModel(mtype, d, n)
model.E, model.nu = hertz_coarse.e_star, 0
- solver = tm.PolonskyKeerRey(model, hertz_coarse.surface, 1e-12,
- request.param, request.param)
+ solver = tm.PolonskyKeerRey(model, hertz_coarse.surface, 1e-12)
models.append(model)
solvers.append(solver)
mat = tm.materials.IsotropicHardening(model,
sigma_y=1e2 * model.E,
hardening=0)
epsolver = DFSANESolver(tm.Residual(model, mat))
- solvers[1] = tm.EPICSolver(solvers[1], epsolver, 1e-12)
+ solvers[1] = request.param(solvers[1], epsolver, 1e-12)
for solver in solvers:
solver.solve(hertz_coarse.load)
return models, hertz_coarse
@pytest.mark.xfail(tm.TamaasInfo.backend == 'tbb',
reason='TBB reductions are unstable?')
def test_energy(model_setup):
"""Test the full elastic-plastic solve step in elasticity"""
# We use computed values from basic PKR to test
(model_elastic, model_plastic), hertz = model_setup
error = norm(
(model_plastic.traction[..., 2] - model_elastic.traction) *
(model_plastic.displacement[0, ..., 2] - model_elastic.displacement))
error /= norm(hertz.pressure * hertz.displacement)
assert error < 1e-16