Diffusion tamaas (maxwell_viscoelastic-branch) Browse

model.cpp
5fab848001d9maxwell_viscoelastic…

model.cpp
View Options

			/**
			* @file
			* @section LICENSE
			*
			* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
			* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
			*
			* This program is free software: you can redistribute it and/or modify
			* it under the terms of the GNU Affero General Public License as published
			* by the Free Software Foundation, either version 3 of the License, or
			* (at your option) any later version.
			*
			* This program is distributed in the hope that it will be useful,
			* but WITHOUT ANY WARRANTY; without even the implied warranty of
			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
			* GNU Affero General Public License for more details.
			*
			* You should have received a copy of the GNU Affero General Public License
			* along with this program. If not, see <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(std::string name) {
			return m.getIntegralOperator(std::move(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)
			.def("computeGradF", &functional::Functional::computeGradF);

			py::class_<functional::AdhesionFunctional> adh(mod, "AdhesionFunctional",
			func);
			adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
			&functional::AdhesionFunctional::setParameters)
			.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)
			.def(py::init<const GridBase<Real>&>(), "surface"_a);
			py::class_<functional::MaugisAdhesionFunctional>(
			mod, "MaugisAdhesionFunctional", adh)
			.def(py::init<const GridBase<Real>&>(), "surface"_a);
			py::class_<functional::SquaredExponentialAdhesionFunctional>(
			mod, "SquaredExponentialAdhesionFunctional", adh)
			.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);
			})
			.def("updateFromModel", &IntegralOperator::updateFromModel)

			.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")
			.value("linear", integration_method::linear)
			.value("cutoff", integration_method::cutoff);
			}

			/// 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; })
			.def_property_readonly_static(
			"components", [](py::object) { return trait::components; })
			.def_property_readonly_static(
			"boundary_dimension",
			[](py::object) { return trait::boundary_dimension; })
			.def_property_readonly_static("voigt",
			[](py::object) { return trait::voigt; })
			.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)
			.value("basic_2d", model_type::basic_2d)
			.value("surface_1d", model_type::surface_1d)
			.value("surface_2d", model_type::surface_2d)
			.value("volume_1d", model_type::volume_1d)
			.value("volume_2d", model_type::volume_2d);

			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)

			.def("setElasticity",
			[](Model& m, Real E, Real nu) {
			TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties");
			m.setElasticity(E, nu);
			},
			"E"_a, "nu"_a)
			.def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"),
			py::return_value_policy::reference_internal)
			.def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
			py::return_value_policy::reference_internal)
			.def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model,
			"boundary_system_size"))
			.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model,
			"boundary_shape"))

			.def("solveNeumann", &Model::solveNeumann)
			.def("solveDirichlet", &Model::solveDirichlet)
			.def("dump", &Model::dump)
			.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())

			.def("getBEEngine",
			[](Model& m) -> decltype(m.getBEEngine()) {
			TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property");
			return m.getBEEngine();
			},
			py::return_value_policy::reference_internal)
			.def("getIntegralOperator",
			[](const Model& m, std::string name) {
			TAMAAS_DEPRECATE("getIntegralOperator()",
			"the operators property");
			return m.getIntegralOperator(std::move(name));
			},
			"operator_name"_a, py::return_value_policy::reference_internal)
			.def("registerField",
			[](Model& m, std::string name, numpy<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(std::move(name));
			},
			"field_name"_a, py::return_value_policy::reference_internal)
			.def("getFields",
			[](const Model& m) {
			TAMAAS_DEPRECATE("getFields()", "list(model)");
			return m.getFields();
			},
			"Return fields list")

			.def("applyElasticity",
			[](Model& model, numpy<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(), std::move(key)) !=
			fields.end();
			},
			py::keep_alive<0, 1>(), "Test field existence")
			.def("__iter__",
			[](const Model& m) {
			const auto& fields = m.getFieldsMap();
			return py::make_key_iterator(fields.cbegin(), fields.cend());
			},
			py::keep_alive<0, 1>(), "Iterator on fields")
			.def_property_readonly(
			"operators", [](Model& m) { return model_operator_accessor{m}; },
			"Returns a dict-like object allowing access to the model's "
			"integral "
			"operators")
			// More python-like access to model properties
			.def_property_readonly("shape", &Model::getDiscretization,
			"Discretization (local in MPI environment)")
			.def_property_readonly("global_shape", &Model::getGlobalDiscretization,
			"Global discretization (in MPI environement)")
			.def_property_readonly("boundary_shape",
			&Model::getBoundaryDiscretization,
			"Number of points on boundary")
			.def_property_readonly("system_size", &Model::getSystemSize,
			"Size of physical domain")
			.def_property_readonly("boundary_system_size",
			&Model::getBoundarySystemSize,
			"Physical size of surface")
			.def_property_readonly("traction",
			(const GridBase<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

tamaas/python/wrap/model.cpp5fab848001d9maxwell_viscoelastic…

model.cppView Options

tamaas/python/wrap/model.cpp
5fab848001d9maxwell_viscoelastic…

model.cpp
View Options