diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index a587d36..2b83ff9 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,109 +1,112 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "logger.hh"
#include "statistics.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
void wrapStatistics(py::module& mod) {
auto name = makeDimensionName("Statistics", dim);
py::class_>(mod, name.c_str())
.def_static("computePowerSpectrum",
&Statistics::computePowerSpectrum,
- py::return_value_policy::move)
- .def_static("computeAutocorrelation",
- &Statistics::computeAutocorrelation,
- py::return_value_policy::move)
- .def_static("computeMoments", &Statistics::computeMoments)
+ py::return_value_policy::move, "Compute PSD of surface")
+ .def_static(
+ "computeAutocorrelation", &Statistics::computeAutocorrelation,
+ py::return_value_policy::move, "Compute autocorrelation of surface")
+ .def_static("computeMoments", &Statistics::computeMoments,
+ "Compute spectral moments")
.def_static("computeSpectralRMSSlope",
- &Statistics::computeSpectralRMSSlope)
- .def_static("computeRMSHeights", &Statistics::computeRMSHeights)
+ &Statistics::computeSpectralRMSSlope,
+ "Compute hrms' in Fourier space")
+ .def_static("computeRMSHeights", &Statistics::computeRMSHeights,
+ "Compute hrms")
.def_static("contact", &Statistics::contact, "tractions"_a,
"perimeter"_a = 0,
"Compute the (corrected) contact area. Permieter is the "
"total contact perimeter in number of segments.");
}
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod) {
// Exposing logging facility
py::enum_(mod, "LogLevel")
.value("debug", LogLevel::debug)
.value("info", LogLevel::info)
.value("warning", LogLevel::warning)
.value("error", LogLevel::error);
mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); });
py::class_(mod, "Logger")
.def(py::init<>())
.def("get",
[](Logger& logger, LogLevel level) -> Logger& {
logger.get(level);
return logger;
})
.def("__lshift__", [](Logger& logger, std::string msg) -> Logger& {
auto level = logger.getWishLevel();
if (level >= Logger::getCurrentLevel() and
not(mpi::rank() != 0 and level == LogLevel::info)) {
py::print(msg, "file"_a = py::module::import("sys").attr("stderr"));
}
return logger;
});
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
mod.def("to_voigt",
[](const Grid& field) {
- Logger().get(LogLevel::warning)
- << "tamaas.to_voigt deprecated. Use tamaas.compute.to_voigt\n";
+ TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()");
+
if (field.getNbComponents() == 9) {
Grid voigt(field.sizes(), 6);
Loop::loop(
[] CUDA_LAMBDA(MatrixProxy in,
SymMatrixProxy out) {
out.symmetrize(in);
},
range>(field),
range>(voigt));
return voigt;
} else
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 87fb854..4f160c7 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,431 +1,451 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "adhesion_functional.hh"
#include "functional.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_extensions.hh"
#include "model_factory.hh"
#include "numpy.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
struct model_operator_accessor {
Model& m;
decltype(auto) get(std::string name) {
return m.getIntegralOperator(std::move(name));
}
};
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_,
functional::wrap::PyFunctional>
func(mod, "Functional");
func.def(py::init<>())
.def("computeF", &functional::Functional::computeF)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
- .def("setParameters",
- [](functional::AdhesionFunctional& f,
- const std::map& m) {
- TAMAAS_DEPRECATE("setParameters()", "the parameters property");
- f.setParameters(m);
- });
+ .def("setParameters", [](functional::AdhesionFunctional& f,
+ const std::map& m) {
+ TAMAAS_DEPRECATE("setParameters()", "the parameters property");
+ f.setParameters(m);
+ });
py::class_(
mod, "ExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
}
template
std::unique_ptr> instanciateFromNumpy(numpy& num) {
std::unique_ptr> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique>>(num);
return result;
case 3:
result = std::make_unique>>(num);
return result;
case 4:
result = std::make_unique>>(num);
return result;
default:
TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
"array to be the number of components");
}
}
/// Wrap IntegralOperator
void wrapIntegralOperator(py::module& mod) {
py::class_(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy input, numpy output) {
TAMAAS_DEPRECATE("apply()", "the () operator");
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"),
py::return_value_policy::reference)
.def(TAMAAS_DEPRECATE_ACCESSOR(getKind, IntegralOperator, "kind"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getType, IntegralOperator, "type"))
.def("__call__",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def("updateFromModel", &IntegralOperator::updateFromModel)
.def_property_readonly("kind", &IntegralOperator::getKind)
.def_property_readonly("model", &IntegralOperator::getModel)
.def_property_readonly("type", &IntegralOperator::getType);
- py::enum_(mod, "integration_method")
- .value("linear", integration_method::linear)
- .value("cutoff", integration_method::cutoff);
+ py::enum_(mod, "integration_method",
+ "Integration method used for the computation "
+ "of volumetric Fourier operators")
+ .value(
+ "linear", integration_method::linear,
+ "No approximation error, O(N_1·N_2·N_3) complexity, may cause float "
+ "overflow/underflow")
+ .value("cutoff", integration_method::cutoff,
+ "Approximation, O(sqrt(N_1^2+N_2^2)·N_3^2) complexity, no "
+ "overflow/underflow risk");
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet)
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, BEEngine, "model"),
py::return_value_policy::reference)
.def_property_readonly("model", &BEEngine::getModel);
}
template
void wrapModelTypeTrait(py::module& mod) {
using trait = model_type_traits;
py::class_(mod, trait::repr)
.def_property_readonly_static("dimension",
- [](py::object) { return trait::dimension; })
+ [](py::object) { return trait::dimension; },
+ "Dimension of computational domain")
.def_property_readonly_static(
- "components", [](py::object) { return trait::components; })
+ "components", [](py::object) { return trait::components; },
+ "Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
- [](py::object) { return trait::boundary_dimension; })
- .def_property_readonly_static("voigt",
- [](py::object) { return trait::voigt; })
+ [](py::object) { return trait::boundary_dimension; },
+ "Dimension of boundary of computational domain")
+ .def_property_readonly_static(
+ "voigt", [](py::object) { return trait::voigt; },
+ "Number of components of symmetrical tensor fields")
.def_property_readonly_static("indices",
[](py::object) { return trait::indices; });
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
- .value("basic_1d", model_type::basic_1d)
- .value("basic_2d", model_type::basic_2d)
- .value("surface_1d", model_type::surface_1d)
- .value("surface_2d", model_type::surface_2d)
- .value("volume_1d", model_type::volume_1d)
- .value("volume_2d", model_type::volume_2d);
+ .value("basic_1d", model_type::basic_1d,
+ "Normal contact with 1D interface")
+ .value("basic_2d", model_type::basic_2d,
+ "Normal contact with 2D interface")
+ .value("surface_1d", model_type::surface_1d,
+ "Normal & tangential contact with 1D interface")
+ .value("surface_2d", model_type::surface_2d,
+ "Normal & tangential contact with 2D interface")
+ .value("volume_1d", model_type::volume_1d,
+ "Contact with volumetric representation and 1D interface")
+ .value("volume_2d", model_type::volume_2d,
+ "Contact with volumetric representation and 2D interface");
auto trait_mod = mod.def_submodule("_type_traits");
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
py::class_(mod, "_model_operator_acessor")
.def(py::init())
.def("__getitem__",
[](model_operator_accessor& acc, std::string name) {
try {
return acc.get(name);
} catch (std::out_of_range&) {
throw py::key_error(name);
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](model_operator_accessor& acc, std::string key) {
const auto ops = acc.m.getIntegralOperators();
return std::find(ops.begin(), ops.end(), key) != ops.end();
})
.def("__iter__",
[](const model_operator_accessor& acc) {
const auto& ops = acc.m.getIntegralOperatorsMap();
return py::make_key_iterator(ops.cbegin(), ops.cend());
},
py::keep_alive<0, 1>());
py::class_(mod, "Model")
.def_property_readonly("type", &Model::getType)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus,
"Young's modulus")
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio,
"Poisson's ratio")
.def_property_readonly("mu", &Model::getShearModulus, "Shear modulus")
.def_property_readonly("E_star", &Model::getHertzModulus,
"Contact (Hertz) modulus")
- .def_property_readonly("be_engine", &Model::getBEEngine)
+ .def_property_readonly("be_engine", &Model::getBEEngine,
+ "Boundary element engine")
.def("setElasticity",
[](Model& m, Real E, Real nu) {
TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties");
m.setElasticity(E, nu);
},
"E"_a, "nu"_a)
.def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model,
"boundary_system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model,
"boundary_shape"))
- .def("solveNeumann", &Model::solveNeumann)
- .def("solveDirichlet", &Model::solveDirichlet)
- .def("dump", &Model::dump)
- .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
+ .def("solveNeumann", &Model::solveNeumann,
+ "Solve surface tractions -> displacements")
+ .def("solveDirichlet", &Model::solveDirichlet,
+ "Solve surface displacemnts -> tractions")
+ .def("dump", &Model::dump, "Write model data to registered dumpers")
+ .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>(),
+ "Register a dumper")
.def("getBEEngine",
[](Model& m) -> decltype(m.getBEEngine()) {
TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property");
return m.getBEEngine();
},
py::return_value_policy::reference_internal)
.def("getIntegralOperator",
[](const Model& m, std::string name) {
TAMAAS_DEPRECATE("getIntegralOperator()",
"the operators property");
return m.getIntegralOperator(std::move(name));
},
"operator_name"_a, py::return_value_policy::reference_internal)
.def("registerField",
[](Model& m, std::string name, numpy field) {
TAMAAS_DEPRECATE("registerField()", "the [] operator");
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def("getField",
[](const Model& m, std::string name) -> decltype(m.getField(name)) {
TAMAAS_DEPRECATE("getField()", "the [] operator");
return m.getField(std::move(name));
},
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields",
[](const Model& m) {
TAMAAS_DEPRECATE("getFields()", "list(model)");
return m.getFields();
},
"Return fields list")
.def("applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
},
"Apply Hooke's law")
// Python magic functions
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def("__getitem__",
[](const Model& m, std::string key) -> decltype(m[key]) {
try {
return m[key];
} catch (std::out_of_range&) {
throw py::key_error(key);
}
},
py::return_value_policy::reference_internal, "Get field")
.def("__setitem__",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
py::keep_alive<1, 3>(), "Register new field")
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), std::move(key)) !=
fields.end();
},
py::keep_alive<0, 1>(), "Test field existence")
.def("__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>(), "Iterator on fields")
.def_property_readonly(
"operators", [](Model& m) { return model_operator_accessor{m}; },
"Returns a dict-like object allowing access to the model's "
"integral "
"operators")
// More python-like access to model properties
.def_property_readonly("shape", &Model::getDiscretization,
"Discretization (local in MPI environment)")
.def_property_readonly("global_shape", &Model::getGlobalDiscretization,
"Global discretization (in MPI environement)")
.def_property_readonly("boundary_shape",
&Model::getBoundaryDiscretization,
"Number of points on boundary")
.def_property_readonly("system_size", &Model::getSystemSize,
"Size of physical domain")
.def_property_readonly("boundary_system_size",
&Model::getBoundarySystemSize,
"Physical size of surface")
.def_property_readonly("traction",
(const GridBase& (Model::*)() const) &
Model::getTraction,
"Surface traction field")
.def_property_readonly("displacement",
(const GridBase& (Model::*)() const) &
Model::getDisplacement,
"Displacement field");
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a, "Dump model")
.def("__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; },
"Dump model");
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static("createModel", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "global_discretization"_a,
"Create a new model of a given type, physical size and "
"*global* discretization")
.def_static("createResidual", &ModelFactory::createResidual, "model"_a,
"sigma_y"_a, "hardening"_a = 0.,
"Create an isotropic linear hardening residual")
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a,
"Register Boussinesq and Mindlin operators to model");
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
"cutoff"_a = 1e-12)
.def_property("yield_stress", &Residual::getYieldStress,
&Residual::setYieldStress)
.def_property("hardening_modulus", &Residual::getHardeningModulus,
&Residual::setHardeningModulus)
.def_property_readonly("model", &Residual::getModel);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 42ca273..eddb3ff 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,188 +1,189 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#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;
void solve() override { PYBIND11_OVERLOAD_PURE(void, EPSolver, solve); }
void updateState() override {
PYBIND11_OVERLOAD(void, EPSolver, updateState);
}
};
/* -------------------------------------------------------------------------- */
void wrapSolvers(py::module& mod) {
py::class_(mod, "ContactSolver")
.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("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("addFunctionalTerm", &ContactSolver::addFunctionalTerm)
+ .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())
.def("solve", py::overload_cast(&ContactSolver::solve),
"target_normal_pressure"_a,
py::call_guard());
py::class_ pkr(mod, "PolonskyKeerRey");
// 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);
py::class_(mod, "KatoSaturated")
.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);
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");
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")
.def(py::init());
py::class_(mod, "EPSolver")
.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&) {
return std::make_unique(res);
}),
"residual"_a, "model"_a, py::keep_alive<1, 2>(),
"Deprecated constructor: only here for compatibility reasons");
py::class_(mod, "EPICSolver")
.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())
.def("acceleratedSolve",
[](EPICSolver& solver, Real pressure) {
return solver.acceleratedSolve(std::vector{pressure});
},
"normal_pressure"_a,
py::call_guard());
}
} // namespace wrap
} // namespace tamaas
diff --git a/src/core/surface_complex.hh b/src/core/surface_complex.hh
deleted file mode 100644
index c9f136b..0000000
--- a/src/core/surface_complex.hh
+++ /dev/null
@@ -1,159 +0,0 @@
-/**
- * @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 .
- *
- */
-/* -------------------------------------------------------------------------- */
-#ifndef SURFACE_COMPLEX_H
-#define SURFACE_COMPLEX_H
-/* -------------------------------------------------------------------------- */
-#include "grid_hermitian.hh"
-#include "surface.hh"
-#include
-
-namespace tamaas {
-
-/**
- * @deprecated should use GridHermitian instead
- * @brief Square 2D complex surface
- */
-template
-class SurfaceComplex : public GridHermitian {
-
-public:
- SurfaceComplex(UInt a, Real L);
- using GridHermitian::GridHermitian;
- using GridHermitian::operator=;
- using GridHermitian::operator();
-
- //! make a surface real by taking norm as real part and imaginary 0
- void makeItRealBySquare();
- //! make it real by summing the imaginary and real parts
- void makeItRealBySum();
- //! make a surface real by taking norm as real part and imaginary 0
- void makeItRealByAbs();
- //! get real part
- Surface real() const;
- //! get imaginary part
- Surface imag() const;
-
- UInt size() const { return this->n[0]; }
- UInt getL() const { return 1.; }
-
- void setGridSize(UInt s);
-
- virtual const SurfaceComplex& getWrappedNumpy() { return *this; };
-
-#define SURF_OPERATOR(op) \
- inline void operator op(const Surface& other); \
- inline void operator op(const SurfaceComplex& other)
-
- SURF_OPERATOR(+=);
- SURF_OPERATOR(*=);
- SURF_OPERATOR(-=);
- SURF_OPERATOR(/=);
-#undef SURF_OPERATOR
-};
-
-/* -------------------------------------------------------------------------- */
-
-template
-void SurfaceComplex::setGridSize(UInt s) {
- this->n[0] = s;
- this->n[1] = s / 2 + 1;
- this->resize(this->n);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template
-Surface SurfaceComplex::real() const {
- Surface res(this->size(), this->getL());
- UInt n = this->size();
- for (UInt i = 0; i < n; ++i)
- for (UInt j = 0; j < n; ++j) {
- res(i, j) = (*this)(i, j).real();
- }
- return res;
-}
-/* -------------------------------------------------------------------------- */
-
-template
-Surface SurfaceComplex::imag() const {
- Surface res(this->size(), this->getL());
- UInt n = this->size();
- for (UInt i = 0; i < n; ++i)
- for (UInt j = 0; j < n; ++j) {
- res(i, j) = (*this)(i, j).imag();
- }
- return res;
-}
-
-/* -------------------------------------------------------------------------- */
-template
-void SurfaceComplex::makeItRealBySquare() {
- Loop::loop([] CUDA_LAMBDA(complex & x) { x *= thrust::conj(x); }, *this);
-}
-
-/* -------------------------------------------------------------------------- */
-template
-void SurfaceComplex::makeItRealByAbs() {
- Loop::loop([] CUDA_LAMBDA(complex & x) { x = thrust::abs(x); }, *this);
-}
-
-/* -------------------------------------------------------------------------- */
-template
-void SurfaceComplex::makeItRealBySum() {
- Loop::loop([] CUDA_LAMBDA(complex & x) { x = x.real() + x.imag(); },
- *this);
-}
-
-/* -------------------------------------------------------------------------- */
-template
-SurfaceComplex::SurfaceComplex(UInt a, Real /*L*/) : GridHermitian() {
- this->n = GridHermitian::hermitianDimensions(std::array{a, a});
- this->resize(this->n);
-}
-
-/* -------------------------------------------------------------------------- */
-#define SURF_OP_IMPL(op) \
- template \
- inline void SurfaceComplex::operator op(const Surface& other) { \
- TAMAAS_ASSERT(this->n[0] == other.sizes()[0] && \
- this->n[1] == other.sizes()[1] / 2 + 1, \
- "surface size does not match"); \
- _Pragma("omp parallel for") for (UInt i = 0; i < this->n[0]; \
- i++) for (UInt j = 0; j < this->n[1]; \
- j++) (*this)(i, j)op other(i, \
- j); \
- } \
- template \
- inline void SurfaceComplex::operator op(const SurfaceComplex& other) { \
- Grid, 2>::operator op(other); \
- }
-
-SURF_OP_IMPL(+=)
-SURF_OP_IMPL(*=)
-SURF_OP_IMPL(-=)
-SURF_OP_IMPL(/=)
-#undef SURF_OP_IMPL
-
-} // namespace tamaas
-
-#endif /* SURFACE_COMPLEX_H */
diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index ee230f4..a615a12 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,143 +1,140 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "kelvin.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
template
Kelvin::Kelvin(Model* model) : VolumePotential(model) {
setIntegrationMethod(integration_method::linear, 1e-12);
}
template
void Kelvin::setIntegrationMethod(integration_method method,
Real cutoff) {
this->method = method;
this->cutoff = cutoff;
Logger logger;
if (this->method == integration_method::linear) {
logger.get(LogLevel::debug)
<< TAMAAS_DEBUG_MSG("Setting linear integration method");
this->initialize(dtrait::template source_components,
dtrait::template out_components,
this->model->getDiscretization()[0]);
}
else {
logger.get(LogLevel::debug) << TAMAAS_DEBUG_MSG(
"Setting cutoff integration method (cutoff " << this->cutoff << ')');
this->initialize(dtrait::template source_components,
dtrait::template out_components, 1);
}
auto max_q = Loop::reduce(
[] CUDA_LAMBDA(VectorProxy qv) {
return qv.l2norm();
},
range>(
this->wavevectors));
if (this->method == integration_method::linear and
not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
logger.get(LogLevel::warning)
<< "Probable overflow of integral computation (consider "
"changing integration method to integration_method::cutoff or "
"compiling with real_type='long double')\n";
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template
void Kelvin::applyIf(GridBase& source,
GridBase& out,
filter_t pred) const {
KelvinInfluence kelvin(this->model->getShearModulus(),
this->model->getPoissonRatio());
parent::transformSource(source, pred);
// Reset buffer values
for (auto&& layer : this->out_buffer)
layer = 0;
if (method == integration_method::linear) {
linearIntegral(out, kelvin);
} else {
cutoffIntegral(out, kelvin);
}
}
template
void Kelvin::linearIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper helper;
helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
this->model->getSystemSize().front(), kelvin);
helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
helper.makeFundamentalGreatAgain(this->out_buffer);
parent::transformOutput(
[](auto&& out_buffer, auto layer) ->
typename parent::BufferType& { return out_buffer[layer]; },
out);
}
template
void Kelvin::cutoffIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper helper;
auto func = [&](auto&& out_buffer, auto layer) ->
typename parent::BufferType& {
auto&& out = out_buffer.front();
out = 0;
helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors,
this->model->getSystemSize().front(), cutoff, kelvin);
helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
helper.makeFundamentalGreatAgain(out);
return out;
};
parent::transformOutput(func, out);
}
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
/* Template instanciation */
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
template class Kelvin;
template class Kelvin;
template class Kelvin;
-/* --------------------------------------------------------------------------
- */
+/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp
index d8f89e5..08f54f6 100644
--- a/src/surface/isopowerlaw.cpp
+++ b/src/surface/isopowerlaw.cpp
@@ -1,101 +1,101 @@
/**
* @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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "isopowerlaw.hh"
#include