Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/python/wrap/model.cpp
2b24d9d43e38maxwell_viscoelastic…

Last Change
Raw File

model.cpp
View Options

/*
* 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)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
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::Functional, std::shared_ptr<functional::Functional>,
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_<functional::AdhesionFunctional> adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters,
"Parameters dictionary")
.def("setParameters", [](functional::AdhesionFunctional& f,
const std::map<std::string, Real>& m) {
TAMAAS_DEPRECATE("setParameters()", "the parameters property");
f.setParameters(m);
});
py::class_<functional::ExponentialAdhesionFunctional>(
mod, "ExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-g/ρ)")
.def(py::init<const GridBase<Real>&>(), "surface"_a);
py::class_<functional::MaugisAdhesionFunctional>(
mod, "MaugisAdhesionFunctional", adh,
"Cohesive zone potential F = H(g - ρ)·γ/ρ")
.def(py::init<const GridBase<Real>&>(), "surface"_a);
py::class_<functional::SquaredExponentialAdhesionFunctional>(
mod, "SquaredExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-0.5·(g/ρ)²)")
.def(py::init<const GridBase<Real>&>(), "surface"_a);
}
template <typename T>
std::unique_ptr<GridBase<T>> instanciateFromNumpy(numpy<T>& num) {
std::unique_ptr<GridBase<T>> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique<GridNumpy<Grid<T, 1>>>(num);
return result;
case 3:
result = std::make_unique<GridNumpy<Grid<T, 2>>>(num);
return result;
case 4:
result = std::make_unique<GridNumpy<Grid<T, 3>>>(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_<IntegralOperator>(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy<Real> input, numpy<Real> 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<Real> input, numpy<Real> 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);
py::enum_<integration_method>(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_<BEEngine>(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 <model_type type>
void wrapModelTypeTrait(py::module& mod) {
using trait = model_type_traits<type>;
py::class_<trait>(mod, trait::repr)
.def_property_readonly_static("dimension",
[](py::object) { return trait::dimension; },
"Dimension of computational domain")
.def_property_readonly_static(
"components", [](py::object) { return trait::components; },
"Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
[](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_<model_type>(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<model_type::basic_1d>(trait_mod);
wrapModelTypeTrait<model_type::surface_1d>(trait_mod);
wrapModelTypeTrait<model_type::volume_1d>(trait_mod);
wrapModelTypeTrait<model_type::basic_2d>(trait_mod);
wrapModelTypeTrait<model_type::surface_2d>(trait_mod);
wrapModelTypeTrait<model_type::volume_2d>(trait_mod);
py::class_<model_operator_accessor>(mod, "_model_operator_acessor")
.def(py::init<Model&>())
.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_<Model>(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,
"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<Real> 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<Real> stress, numpy<Real> 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<Real> 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_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<Real>& (Model::*)() const) &
Model::getTraction,
"Surface traction field")
.def_property_readonly("displacement",
(const GridBase<Real>& (Model::*)() const) &
Model::getDisplacement,
"Displacement field");
py::class_<ModelDumper, PyModelDumper, std::shared_ptr<ModelDumper>>(
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_<ModelFactory>(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_<Residual, PyResidual>(mod, "Residual")
.def(py::init<Model*>())
.def("computeResidual",
[](Residual& res, numpy<Real>& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy<Real>& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy<Real>& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy<Real>& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy<Real>& output, numpy<Real>& input,
numpy<Real>& 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
c4science · Help