diff --git a/INFOS.py b/INFOS.py
index 0b9b760..28ba298 100755
--- a/INFOS.py
+++ b/INFOS.py
@@ -1,54 +1,55 @@
#!/usr/bin/env python3
"""Defines the information to be used throughout the builds."""
# -*- mode:python; coding: utf-8 -*-
from collections import namedtuple
import versioneer
TamaasInfo = namedtuple('TamaasInfo',
['version',
'authors',
'maintainer',
'email',
'copyright',
'description',
'license'])
TAMAAS_INFOS = TamaasInfo(
version=versioneer.get_version(),
authors=([
u'Lucas Frérot',
'Guillaume Anciaux',
'Valentine Rey',
'Son Pham-Ba',
u'Jean-François Molinari'
]),
maintainer=u'Lucas Frérot',
email='lucas.frerot@imtek.uni-freiburg.de',
copyright=(
- u"Copyright (©) 2016-2022 EPFL "
- u"(École Polytechnique Fédérale de Lausanne), "
- u"Laboratory (LSMS - Laboratoire de Simulation en "
- u"Mécanique des Solides)"
+ "Copyright (©) 2016-2020 EPFL "
+ "(École Polytechnique Fédérale de Lausanne), "
+ "Laboratory (LSMS - Laboratoire de Simulation en "
+ "Mécanique des Solides)\\n"
+ "Copyright (©) 2020-2022 Lucas Frérot"
),
description='A high-performance library for periodic rough surface contact',
license="SPDX-License-Identifier: AGPL-3.0-or-later",
)
def main():
import argparse
parser = argparse.ArgumentParser(description="Print Tamaas info")
parser.add_argument("--version", action="store_true")
args = parser.parse_args()
if args.version:
print(TAMAAS_INFOS.version)
if __name__ == "__main__":
main()
diff --git a/python/cast.hh b/python/cast.hh
index e5cc85c..f99534c 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,196 +1,197 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 CAST_HH
#define CAST_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "grid_base.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace pybind11 {
// Format descriptor necessary for correct wrap of tamaas complex type
template
struct format_descriptor<
tamaas::complex, detail::enable_if_t::value>> {
static constexpr const char c = format_descriptor::c;
static constexpr const char value[3] = {'Z', c, '\0'};
static std::string format() { return std::string(value); }
};
#ifndef PYBIND11_CPP17
template
constexpr const char format_descriptor<
tamaas::complex,
detail::enable_if_t::value>>::value[3];
#endif
namespace detail {
// declare tamaas complex as a complex type for pybind11
template
struct is_complex> : std::true_type {};
template
struct is_fmt_numeric,
detail::enable_if_t::value>>
: std::true_type {
static constexpr int index = is_fmt_numeric::index + 3;
};
static inline handle policy_switch(return_value_policy policy, handle parent) {
switch (policy) {
case return_value_policy::copy:
case return_value_policy::move:
return handle();
case return_value_policy::automatic_reference: // happens in python-derived
// classes
case return_value_policy::reference:
return none();
case return_value_policy::reference_internal:
return parent;
default:
TAMAAS_EXCEPTION("Policy is not handled");
}
}
template
handle grid_to_python(const tamaas::Grid& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes(dim);
std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
if (grid.getNbComponents() != 1)
sizes.push_back(grid.getNbComponents());
return array(sizes, grid.getInternalData(), parent).release();
}
template
handle grid_to_python(const tamaas::GridBase& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes = {grid.dataSize()};
return array(sizes, grid.getInternalData(), parent).release();
}
/**
* Type caster for grid classes
* inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
* and pybind11/eigen.h
*/
template class G, typename T,
tamaas::UInt dim>
struct type_caster> {
using type = G;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridWrap"));
/**
* Conversion part 1 (Python->C++): convert a PyObject into a grid
* instance or return false upon failure. The second argument
* indicates whether implicit conversions should be applied.
*/
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf)
return false;
value.move(tamaas::wrap::GridNumpy>(buf));
return true;
}
/**
* Conversion part 2 (C++ -> Python): convert a grid instance into
* a Python object. The second and third arguments are used to
* indicate the return value policy and parent object (for
* ``return_value_policy::reference_internal``) and are generally
* ignored by implicit casters.
*/
static handle cast(const type& src, return_value_policy policy,
handle parent) {
return grid_to_python(src, policy, parent);
}
};
/**
* Type caster for GridBase classes
*/
template
struct type_caster> {
using type = tamaas::GridBase;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridBaseWrap"));
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
if (!buf)
return false;
value.move(tamaas::wrap::GridBaseNumpy(buf));
return true;
}
static handle cast(const type& src, return_value_policy policy,
handle parent) {
#define GRID_BASE_CASE(unused1, unused2, dim) \
case dim: { \
const auto* conv = dynamic_cast*>(&src); \
if (conv) \
return grid_to_python(*conv, policy, parent); \
else \
return grid_to_python(src, policy, parent); \
}
switch (src.getDimension()) {
BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
default:
return grid_to_python(src, policy, parent);
}
#undef GRID_BASE_CASE
}
};
} // namespace detail
} // namespace pybind11
#endif // CAST_HH
diff --git a/python/numpy.hh b/python/numpy.hh
index 0e8f591..d4e6aa1 100644
--- a/python/numpy.hh
+++ b/python/numpy.hh
@@ -1,92 +1,93 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 NUMPY_HH
#define NUMPY_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
#include
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// Type alias for usual numpy array type
template
using numpy = py::array_t;
/// GridBase helper class wrapping numpy array
template
class GridBaseNumpy : public GridBase {
public:
GridBaseNumpy(numpy& buffer) : GridBase() {
this->nb_components = 1;
this->data.wrapMemory(buffer.mutable_data(), buffer.size());
}
};
/// Grid helper class wrapping numpy array
template
class GridNumpy : public Parent {
public:
GridNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
const UInt ndim = buffer.ndim();
if (ndim > Parent::dimension + 1 || ndim < Parent::dimension)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected grid dimensions");
if (ndim == Parent::dimension + 1)
this->nb_components = array_shape[ndim - 1];
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
/// Surface helper class wrapping numpy array
template
class SurfaceNumpy : public Parent {
public:
SurfaceNumpy(numpy& buffer) : Parent() {
auto* array_shape = buffer.shape();
if (buffer.ndim() != 2)
TAMAAS_EXCEPTION(
"Numpy array dimension do not match expected surface dimensions");
std::copy_n(array_shape, Parent::dimension, this->n.begin());
this->computeStrides();
this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
}
};
} // namespace wrap
} // namespace tamaas
#endif // NUMPY_HH
diff --git a/python/tamaas_module.cpp b/python/tamaas_module.cpp
index e2445b3..554dfe7 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,86 +1,87 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "tamaas.hh"
#include "tamaas_info.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace py = pybind11;
namespace detail {
template
struct dtype_helper {
static const py::dtype dtype;
};
template <>
const py::dtype dtype_helper::dtype("=f8");
template <>
const py::dtype dtype_helper::dtype("=f16");
} // namespace detail
/// Creating the tamaas python module
PYBIND11_MODULE(_tamaas, mod) {
mod.doc() = "Compiled component of Tamaas";
// Wrapping the base methods
mod.def("initialize", &initialize, py::arg("num_threads") = 0,
"Initialize tamaas with desired number of threads. Automatically "
"called upon import of the tamaas module, but can be manually called "
"to set the desired number of threads.");
mod.def("finalize", []() {
PyErr_WarnEx(PyExc_DeprecationWarning,
"finalize() is deprecated, it is now automatically called "
"at the end of the program", 1);
});
// Default dtype of numpy arrays
mod.attr("dtype") = detail::dtype_helper::dtype;
// Wrapping release information
auto info = py::class_(mod, "TamaasInfo");
info.attr("version") = TamaasInfo::version;
info.attr("build_type") = TamaasInfo::build_type;
info.attr("branch") = TamaasInfo::branch;
info.attr("commit") = TamaasInfo::commit;
info.attr("diff") = TamaasInfo::diff;
info.attr("remotes") = TamaasInfo::remotes;
info.attr("has_mpi") = TamaasInfo::has_mpi;
info.attr("backend") = TamaasInfo::backend;
// Wrapping tamaas components
wrap::wrapCore(mod);
wrap::wrapPercolation(mod);
wrap::wrapSurface(mod);
wrap::wrapModel(mod);
wrap::wrapSolvers(mod);
wrap::wrapCompute(mod);
wrap::wrapMPI(mod);
/// Wrapping test features
wrap::wrapTestFeatures(mod);
}
} // namespace tamaas
diff --git a/python/wrap.hh b/python/wrap.hh
index c6ac47d..4a18d8f 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,76 +1,77 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 WRAP_HH
#define WRAP_HH
/* -------------------------------------------------------------------------- */
#include "cast.hh"
#include "numpy.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
namespace py = pybind11;
/// For naming classes templated with dimension
inline std::string makeDimensionName(const std::string& name, UInt dim) {
std::stringstream str;
str << name << dim << "D";
return str.str();
}
/* -------------------------------------------------------------------------- */
/* Prototypes for wrapping of tamaas components */
/* -------------------------------------------------------------------------- */
void wrapCore(py::module& mod);
void wrapPercolation(py::module& mod);
void wrapSurface(py::module& mod);
void wrapModel(py::module& mod);
void wrapSolvers(py::module& mod);
void wrapTestFeatures(py::module& mod);
void wrapCompute(py::module& mod);
void wrapMPI(py::module& mod);
/* -------------------------------------------------------------------------- */
/* Deprecation macro */
/* -------------------------------------------------------------------------- */
#define TAMAAS_DEPRECATE(olds, news) \
do { \
PyErr_WarnEx(PyExc_DeprecationWarning, \
olds " is deprecated, use " news " instead.", 1); \
} while (0)
#define TAMAAS_DEPRECATE_ACCESSOR(acc, type, property) \
#acc, [](const type& m) -> decltype(m.acc()) { \
TAMAAS_DEPRECATE(#acc "()", "the " property " property"); \
return m.acc(); \
}
} // namespace wrap
} // namespace tamaas
#endif // WRAP_HH
diff --git a/python/wrap/compute.cpp b/python/wrap/compute.cpp
index 666e9d5..fee8395 100644
--- a/python/wrap/compute.cpp
+++ b/python/wrap/compute.cpp
@@ -1,90 +1,91 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "computes.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
void wrapCompute(py::module& mod) {
auto compute_mod = mod.def_submodule("compute");
compute_mod.def(
"eigenvalues",
[](model_type type, Grid& eigs, const Grid& field) {
eigenvalues(type, eigs, field);
},
"model_type"_a, "eigenvalues_out"_a, "field"_a,
"Compute eigenvalues of a tensor field");
compute_mod.def(
"von_mises",
[](model_type type, Grid& vm, const Grid& field) {
vonMises(type, vm, field);
},
"model_type"_a, "von_mises"_a, "field"_a,
"Compute the Von Mises invariant of a tensor field");
compute_mod.def(
"deviatoric",
[](model_type type, Grid& dev, const Grid& field) {
deviatoric(type, dev, field);
},
"model_type"_a, "deviatoric"_a, "field"_a,
"Compute the deviatoric part of a tensor field");
compute_mod.def(
"to_voigt",
[](const Grid& field) {
if (field.getNbComponents() != 9)
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
Grid voigt(field.sizes(), 6);
Loop::loop(
[] CUDA_LAMBDA(MatrixProxy in,
SymMatrixProxy out) { out.symmetrize(in); },
range>(field),
range>(voigt));
return voigt;
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
compute_mod.def(
"from_voigt",
[](const Grid& field) {
if (field.getNbComponents() != 6)
TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
Grid full(field.sizes(), 9);
Loop::loop([] CUDA_LAMBDA(
SymMatrixProxy in,
MatrixProxy out) { out.fromSymmetric(in); },
range>(field),
range>(full));
return full;
},
"Convert a 3D tensor field to voigt notation",
py::return_value_policy::copy);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index e61c995..05a5a10 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,120 +1,121 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "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, "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,
"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.")
.def_static("graphArea", &Statistics::graphArea, "zdisplacement"_a,
"Compute area defined by function graph.");
}
/* -------------------------------------------------------------------------- */
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); });
mod.def("get_log_level", Logger::getCurrentLevel);
py::class_(mod, "Logger")
.def(py::init<>())
.def(
"get",
[](Logger& logger, LogLevel level) -> Logger& {
logger.get(level);
return logger;
},
"Get a logger object for a log level")
.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;
},
"Log a message");
wrapStatistics<1>(mod);
wrapStatistics<2>(mod);
mod.def(
"to_voigt",
[](const Grid& field) {
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;
}
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 3296dd5..ead4a54 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,516 +1,517 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "model.hh"
#include "adhesion_functional.hh"
#include "elastic_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(const std::string& name) {
return m.getIntegralOperator(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,
"Compute functional value")
.def("computeGradF", &functional::Functional::computeGradF,
"Compute functional gradient");
py::class_(
mod, "ElasticFunctionalPressure", func)
.def(py::init&>());
py::class_(mod, "ElasticFunctionalGap",
func)
.def(py::init&>());
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters,
"Parameters dictionary")
.def("setParameters", [](functional::AdhesionFunctional& f,
const std::map& m) {
TAMAAS_DEPRECATE("setParameters()", "the parameters property");
f.setParameters(m);
});
py::class_(
mod, "ExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-g/ρ)")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh,
"Cohesive zone potential F = H(g - ρ)·γ/ρ")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-0.5·(g/ρ)²)")
.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);
},
"Apply the integral operator")
.def("updateFromModel", &IntegralOperator::updateFromModel,
"Resets internal persistent variables from the model")
.def_property_readonly("kind", &IntegralOperator::getKind)
.def_property_readonly("model", &IntegralOperator::getModel)
.def_property_readonly("type", &IntegralOperator::getType)
.def_property_readonly("shape", &IntegralOperator::matvecShape)
.def(
"matvec",
[](const IntegralOperator& op, numpy X) -> GridBase {
GridBaseNumpy x(X);
auto y = op.matvec(x);
return y;
},
py::return_value_policy::move);
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₁·N₂·N₃) time complexity, may cause "
"float overflow/underflow")
.value("cutoff", integration_method::cutoff,
"Approximation, O(sqrt(N₁²+N₂²)·N₃²) time 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", [](const py::object&) { return trait::dimension; },
"Dimension of computational domain")
.def_property_readonly_static(
"components", [](const py::object&) { return trait::components; },
"Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
[](const py::object&) { return trait::boundary_dimension; },
"Dimension of boundary of computational domain")
.def_property_readonly_static(
"voigt", [](const py::object&) { return trait::voigt; },
"Number of components of symmetrical tensor fields")
.def_property_readonly_static(
"indices", [](const py::object&) { return trait::indices; });
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
.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(py::init(py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel)))
.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,
"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,
"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(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(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(), 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(
"__copy__",
[](const Model&) {
throw std::runtime_error("__copy__ not implemented");
},
"Shallow copy of model. Not implemented.")
.def(
"__deepcopy__",
[](const Model& m, py::dict) { return ModelFactory::createModel(m); },
"Deep copy of model.")
.def_property_readonly("boundary_fields", &Model::getBoundaryFields)
.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",
py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel),
"model_type"_a, "system_size"_a, "global_discretization"_a,
R"-(Create a new model of a given type, physical size and *global*
discretization.
:param model_type: the type of desired model
:param system_size: the physical size of the domain in each direction
:param global_discretization: number of points in each direction)-")
.def_static("createModel",
py::overload_cast(&ModelFactory::createModel),
"model"_a, "Create a deep copy of a model.")
.def_static("createResidual", &ModelFactory::createResidual, "model"_a,
"sigma_y"_a, "hardening"_a = 0.,
R"-(Create an isotropic linear hardening residual.
:param model: the model on which to define the residual
:param sigma_y: the (von Mises) yield stress
:param hardening: the hardening modulus)-")
.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/model_extensions.hh b/python/wrap/model_extensions.hh
index 9ff457d..09debbd 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,147 +1,148 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "functional.hh"
#include "model.hh"
#include "model_dumper.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace functional {
namespace wrap {
/// Class for extension of Functional in python
class PyFunctional : public Functional {
public:
using Functional::Functional;
// Overriding pure virtual functions
Real computeF(GridBase& variable, GridBase& dual) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
}
void computeGradF(GridBase& variable,
GridBase& gradient) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
}
};
} // namespace wrap
} // namespace functional
/* -------------------------------------------------------------------------- */
namespace wrap {
/* -------------------------------------------------------------------------- */
class PyModelDumper : public ModelDumper {
public:
using ModelDumper::ModelDumper;
void dump(const Model& model) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, ModelDumper, dump, model);
}
};
/* -------------------------------------------------------------------------- */
class PyResidual : public Residual {
public:
using Residual::Residual;
void computeResidual(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidual, strain_increment);
}
void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, applyTangent, output, input,
current_strain_increment);
}
void computeStress(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeStress, strain_increment);
}
void updateState(GridBase& converged_strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, updateState,
converged_strain_increment);
}
const GridBase& getVector() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getVector);
}
const GridBase& getPlasticStrain() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getPlasticStrain);
}
const GridBase& getStress() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getStress);
}
void computeResidualDisplacement(GridBase& strain_increment) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidualDisplacement,
strain_increment);
}
void setIntegrationMethod(integration_method method, Real cutoff) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setIntegrationMethod, method,
cutoff);
}
Real getYieldStress() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Residual, getYieldStress);
}
Real getHardeningModulus() const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(Real, Residual, getHardeningModulus);
}
void setYieldStress(Real val) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setYieldStress, val);
}
void setHardeningModulus(Real val) override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Residual, setHardeningModulus, val);
}
};
/* -------------------------------------------------------------------------- */
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/mpi.cpp b/python/wrap/mpi.cpp
index d0c2a83..b36ad30 100644
--- a/python/wrap/mpi.cpp
+++ b/python/wrap/mpi.cpp
@@ -1,100 +1,101 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "mpi_interface.hh"
#include "partitioner.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
void wrapMPI(py::module& mod) {
auto mpi_mod = mod.def_submodule("mpi");
py::class_(mpi_mod, "sequential")
.def(py::init<>())
.def("__enter__",
[](mpi::sequential& /*obj*/) { mpi::sequential::enter(); })
.def("__exit__", [](mpi::sequential& /*obj*/, py::object /*unused*/,
py::object /*unused*/,
py::object /*unused*/) { mpi::sequential::exit(); });
mpi_mod.def(
"local_shape",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_size(global);
case 2:
return Partitioner<2>::local_size(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a, "Gives the local size of a 1D/2D global shape");
mpi_mod.def(
"global_shape",
[](std::vector local) {
switch (local.size()) {
case 1:
return Partitioner<1>::global_size(local);
case 2:
return Partitioner<2>::global_size(local);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"local_shape"_a, "Gives the global shape of a 1D/2D local shape");
mpi_mod.def(
"local_offset",
[](std::vector global) {
switch (global.size()) {
case 1:
return Partitioner<1>::local_offset(global);
case 2:
return Partitioner<2>::local_offset(global);
default:
TAMAAS_EXCEPTION("Please provide a 1D/2D shape");
}
},
"global_shape"_a, "Gives the local offset of a 1D/2D global shape");
mpi_mod.def("gather", &Partitioner<2>::gather,
py::return_value_policy::move, "Gather 2D surfaces");
mpi_mod.def("scatter", &Partitioner<2>::scatter,
py::return_value_policy::move, "Scatter 2D surfaces");
mpi_mod.def(
"size", []() { return mpi::size(); },
"Returns the number of MPI processes");
mpi_mod.def(
"rank", []() { return mpi::rank(); },
"Returns the rank of the local process");
}
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/python/wrap/percolation.cpp b/python/wrap/percolation.cpp
index 210f17c..fd1aa79 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,92 +1,93 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "flood_fill.hh"
#include "wrap.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
/* -------------------------------------------------------------------------- */
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
void wrapCluster(py::module& mod) {
auto name = makeDimensionName("Cluster", dim);
py::class_>(mod, name.c_str())
.def(py::init<>())
.def_property_readonly("area", &Cluster::getArea, "Area of cluster")
.def_property_readonly("points", &Cluster::getPoints,
"Get list of points of cluster")
.def_property_readonly("perimeter", &Cluster::getPerimeter,
"Get perimeter of cluster")
.def_property_readonly("bounding_box", &Cluster::boundingBox,
"Compute the bounding box of a cluster")
.def(TAMAAS_DEPRECATE_ACCESSOR(getArea, Cluster, "area"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPoints, Cluster, "points"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPerimeter, Cluster, "perimeter"))
.def("__str__", [](const Cluster& cluster) {
std::stringstream sstr;
sstr << '{';
auto&& points = cluster.getPoints();
auto print_point = [&sstr](const auto& point) {
sstr << '(';
for (UInt i = 0; i < point.size() - 1; ++i)
sstr << point[i] << ", ";
sstr << point.back() << ')';
};
auto point_it = points.begin();
for (UInt p = 0; p < points.size() - 1; ++p) {
print_point(*(point_it++));
sstr << ", ";
}
print_point(points.back());
sstr << "}";
return sstr.str();
});
}
void wrapPercolation(py::module& mod) {
wrapCluster<1>(mod);
wrapCluster<2>(mod);
wrapCluster<3>(mod);
py::class_(mod, "FloodFill")
.def_static("getSegments", &FloodFill::getSegments,
"Return a list of segments from boolean map", "contact"_a)
.def_static("getClusters", &FloodFill::getClusters,
"Return a list of clusters from boolean map", "contact"_a,
"diagonal"_a)
.def_static("getVolumes", &FloodFill::getVolumes,
"Return a list of volume clusters", "map"_a, "diagonal"_a);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/solvers.cpp b/python/wrap/solvers.cpp
index 769718d..3f9e3c0 100644
--- a/python/wrap/solvers.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,238 +1,239 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "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_readonly("functional", &ContactSolver::getFunctional)
.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);
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&) {
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",
"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!");
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index 07ec191..2f8a0f4 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,177 +1,178 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * Copyright (©) 2020-2022 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 "isopowerlaw.hh"
#include "regularized_powerlaw.hh"
#include "surface_generator_filter.hh"
#include "surface_generator_random_phase.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace wrap {
using namespace py::literals;
/* -------------------------------------------------------------------------- */
template
class PyFilter : public Filter {
public:
using Filter::Filter;
// Overriding pure virtual functions
void
computeFilter(GridHermitian& filter_coefficients) const override {
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_OVERLOAD_PURE(void, Filter, computeFilter,
filter_coefficients);
}
};
template
void wrapFilter(py::module& mod) {
auto name = makeDimensionName("Filter", dim);
py::class_, std::shared_ptr>, PyFilter>(
mod, name.c_str(), "Mother class for Fourier filter objects")
.def(py::init<>())
.def("computeFilter",
(void (Filter::*)(GridHermitian&) const) &
Filter::computeFilter,
"Compute the Fourier coefficient of the surface");
}
/* -------------------------------------------------------------------------- */
template
void wrapIsopowerlaw(py::module& mod) {
std::string name = makeDimensionName("Isopowerlaw", dim);
py::class_, Filter, std::shared_ptr>>(
mod, name.c_str(), "Isotropic powerlaw spectrum with a rolloff plateau")
.def(py::init<>())
.def_property("q0", &Isopowerlaw::getQ0, &Isopowerlaw::setQ0,
"Long wavelength cutoff")
.def_property("q1", &Isopowerlaw::getQ1, &Isopowerlaw::setQ1,
"Rolloff wavelength")
.def_property("q2", &Isopowerlaw::getQ2, &Isopowerlaw::setQ2,
"Short wavelength cutoff")
.def_property("hurst", &Isopowerlaw::getHurst,
&Isopowerlaw