model.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, Apr 26, 08:01

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 "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 <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::ElasticFunctionalPressure>(
	mod, "ElasticFunctionalPressure", func)
	.def(py::init<const IntegralOperator&, const GridBase<Real>&>());
	py::class_<functional::ElasticFunctionalGap>(mod, "ElasticFunctionalGap",
	func)
	.def(py::init<const IntegralOperator&, const GridBase<Real>&>());

	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)
	.def_property_readonly("shape", &IntegralOperator::matvecShape)
	.def(
	"matvec",
	[](const IntegralOperator& op, numpy<Real> X) -> GridBase<Real> {
	GridBaseNumpy<Real> x(X);
	auto y = op.matvec(x);
	return y;
	},
	py::return_value_policy::move);

	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(py::init(py::overload_cast<model_type, const std::vector<Real>&,
	const std::vector<UInt>&>(
	&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<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(
	"__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<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",
	py::overload_cast<model_type, const std::vector<Real>&,
	const std::vector<UInt>&>(
	&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<const Model&>(&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_<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

model.cppNo OneTemporaryActions

File Metadata

model.cppView Options

Event Timeline

model.cpp
No OneTemporary
Actions

model.cpp
View Options