diff --git a/examples/new_adhesion.py b/examples/new_adhesion.py index bb3a09e..1785ac3 100644 --- a/examples/new_adhesion.py +++ b/examples/new_adhesion.py @@ -1,80 +1,103 @@ import tamaas as tm import sys import matplotlib.pyplot as plt import numpy as np +import argparse + +class AdhesionPython(tm.Functional): + def __init__(self, engine, rho, gamma): + tm.Functional.__init__(self, engine) + self.rho = rho + self.gamma = gamma + + def computeF(self, gap, pressure): + return -self.gamma * np.sum(np.exp(-gap / self.rho)) + + def computeGradF(self, gap, gradient): + gradient += self.gamma * np.exp(-gap / self.rho) / self.rho + + +parser = argparse.ArgumentParser() +parser.add_argument('--local-functional', dest="py_adh", action="store_true", + help="use the adhesion functional written in python") + +args = parser.parse_args() # Initialize threads and fftw tm.initialize() # Surface size n = 1024 # Surface generator sg = tm.SurfaceGeneratorFilter2D() sg.setSizes([n, n]) sg.setRandomSeed(0) # Spectrum spectrum = tm.Isopowerlaw2D() sg.setFilter(spectrum) # Parameters params = { "Q0": 16, "Q1": 16, "Q2": 64, "Hurst": 0.8 } helper = tm.ParamHelper(spectrum) helper.set(params) # Generating surface surface = sg.buildSurface() #surface /= tm.SurfaceStatistics.computeSpectralRMSSlope(surface) surface /= n #print(spectrum.rmsSlopes()) #print(tm.SurfaceStatistics.computeRMSSlope(surface)) plt.imshow(surface) # Creating model model = tm.ModelFactory.createModel(tm.model_type_basic_2d, [1., 1.], [n, n]) # Solver solver = tm.PolonskyKeerRey(model, surface, 1e-12, tm.PolonskyKeerRey.gap, tm.PolonskyKeerRey.gap) - -adhesion = tm.ExponentialAdhesionFunctional(model.getBEEngine(), surface) - adhesion_params = { - "rho":2e-3, - "surface_energy":2e-5 + "rho":2e-3, + "surface_energy":2e-5 } -adhesion.setParameters(adhesion_params) +if args.py_adh: + adhesion = AdhesionPython(model.getBEEngine(), + adhesion_params["rho"], + adhesion_params["surface_energy"]) +else: + adhesion = tm.ExponentialAdhesionFunctional(model.getBEEngine(), surface) + adhesion.setParameters(adhesion_params) solver.addFunctionalTerm(adhesion) # Solve for target pressure g_target = 5e-2 solver.solve(g_target) tractions = model.getTraction() plt.figure() plt.imshow(tractions) plt.colorbar() plt.figure() zones = np.zeros_like(tractions) tol = 1e-6 zones[tractions > tol] = 1 zones[tractions < -tol] = -1 plt.imshow(zones, cmap='Greys') # Cleanup threads tm.finalize() plt.show() diff --git a/python/wrap/bem.cpp b/python/wrap/bem.cpp index 6033b83..815c14c 100644 --- a/python/wrap/bem.cpp +++ b/python/wrap/bem.cpp @@ -1,106 +1,111 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "wrap.hh" #include "bem_fft_base.hh" #include "bem_polonski.hh" #include "bem_kato.hh" #include "bem_gigipol.hh" #include "bem_grid.hh" #include "bem_grid_polonski.hh" +#include "bem_grid_condat.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { using namespace py::literals; /* -------------------------------------------------------------------------- */ void wrapBEM(py::module& mod) { py::class_(mod, "BemFFTBase") .def("computeDisplacementsFromTractions", &BemFFTBase::computeDisplacementsFromTractions) .def("computeTractionsFromDisplacements", &BemFFTBase::computeTractionsFromDisplacements) .def("applyInfluenceFunctions", &BemFFTBase::applyInfluenceFunctions) .def("applyInverseInfluenceFunctions", &BemFFTBase::applyInverseInfluenceFunctions) .def("computeSpectralInfluenceOverDisplacements", &BemFFTBase::computeSpectralInfluenceOverDisplacement) .def("computeTrueDisplacements", &BemFFTBase::computeTrueDisplacements) .def("computeGaps", &BemFFTBase::computeGaps) .def("getSpectralInfluenceOverDisplacement", &BemFFTBase::getSpectralInfluenceOverDisplacement) .def("getTractions", &BemFFTBase::getTractions) .def("getDisplacements", &BemFFTBase::getDisplacements) .def("getTrueDisplacements", &BemFFTBase::getTrueDisplacements) .def("getGap", &BemFFTBase::getGap) .def("setEffectiveModulus", &BemFFTBase::setEffectiveModulus, "E_star"_a) .def("getEffectiveModulus", &BemFFTBase::getEffectiveModulus) .def("setMaxIterations", &BemFFTBase::setMaxIterations) .def("getMaxIterations", &BemFFTBase::getMaxIterations) .def("setDumpFreq", &BemFFTBase::setDumpFreq) .def("getDumpFreq", &BemFFTBase::getDumpFreq); py::class_(mod, "BemPolonski") .def(py::init&>()) .def("computeEquilibrium", &BemPolonski::computeEquilibrium) .def("computeEquilibriuminit", &BemPolonski::computeEquilibriuminit) .def("computeMeanGapsInContact", &BemPolonski::computeMeanGapsInContact) .def("setMaxPressure", &BemPolonski::setMaxPressure) .def("getNbIterations", &BemPolonski::getNbIterations); py::class_(mod, "BemKato") .def(py::init&>()) .def("computeEquilibrium", &BemKato::computeEquilibrium) .def("setMaxPressure", &BemKato::setMaxPressure); py::class_(mod, "BemGigiPol") .def(py::init&>()) .def("computeEquilibrium", &BemGigipol::computeEquilibrium) .def("computeEquilibrium2", &BemGigipol::computeEquilibrium2) .def("computeEquilibriuminit", &BemGigipol::computeEquilibriuminit) .def("computeEquilibrium2init", &BemGigipol::computeEquilibrium2init) .def("computeTractionsFromDisplacements", &BemGigipol::computeTractionsFromDisplacements); py::class_(mod, "BemGrid") .def("computeDisplacementsFromTractions", &BemGrid::computeDisplacementsFromTractions) .def("computeSurfaceNormals", &BemGrid::computeSurfaceNormals) .def("computeInfluence", &BemGrid::computeInfluence) .def("setElasticity", &BemGrid::setElasticity) .def("getTractions", &BemGrid::getTractions) .def("getDisplacements", &BemGrid::getDisplacements) .def("getSurfaceNormals", &BemGrid::getSurfaceNormals); py::class_(mod, "BemGridPolonski") .def(py::init&>()); + + py::class_(mod, "BemGridCondat") + .def(py::init&>()) + .def("computeEquilibrium", &BemGridCondat::computeEquilibrium); } } // namespace wrap __END_TAMAAS__ diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index 52d4239..0768f4f 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,136 +1,137 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "wrap.hh" #include "model.hh" #include "model_factory.hh" #include "functional.hh" #include "adhesion_functional.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace functional { namespace wrap { class PyFunctional : public Functional { public: using Functional::Functional; // Overriding pure virtual functions Real computeF(GridBase& variable, GridBase& dual) const override { PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual); } void computeGradF(GridBase& variable, GridBase& gradient) const override { PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient); } }; } } namespace wrap { using namespace py::literals; void wrapFunctionals(py::module& mod) { py::class_ func( mod, "Functional"); - func.def("computeF", &functional::Functional::computeF) + func.def(py::init()) + .def("computeF", &functional::Functional::computeF) .def("computeGradF", &functional::Functional::computeGradF); py::class_ adh(mod, "AdhesionFunctional", func); adh.def_property("parameters", &functional::AdhesionFunctional::getParameters, &functional::AdhesionFunctional::setParameters) // legacy wrapper code .def("setParameters", &functional::AdhesionFunctional::setParameters); py::class_( mod, "ExponentialAdhesionFunctional", adh) .def(py::init&>(), "be_engine"_a, "surface"_a); py::class_( mod, "MaugisAdhesionFunctional", adh) .def(py::init&>(), "be_engine"_a, "surface"_a); py::class_( mod, "SquaredExponentialAdhesionFunctional", adh) .def(py::init&>(), "be_engine"_a, "surface"_a); } void wrapBEEngine(py::module& mod) { py::class_(mod, "BEEngine") .def("solveNeumann", &BEEngine::solveNeumann) .def("solveDirichlet", &BEEngine::solveDirichlet) .def("getModel", &BEEngine::getModel, py::return_value_policy::reference); } void wrapModelClass(py::module& mod) { py::enum_(mod, "model_type") .value("basic_1d", model_type::basic_1d) .value("basic_2d", model_type::basic_2d) .value("surface_1d", model_type::surface_1d) .value("surface_2d", model_type::surface_2d) .value("volume_1d", model_type::volume_1d) .value("volume_2d", model_type::volume_2d); py::class_(mod, "Model") .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a) .def("getHertzModulus", &Model::getHertzModulus) .def("getYoungModulus", &Model::getYoungModulus) .def("getPoissonRatio", &Model::getPoissonRatio) .def("getTraction", (GridBase & (Model::*)()) & Model::getTraction) .def("getDisplacement", (GridBase & (Model::*)()) & Model::getDisplacement) .def("getSystemSize", &Model::getSystemSize) .def("getDiscretization", &Model::getDiscretization) .def("solveNeumann", &Model::solveNeumann) .def("solveDirichlet", &Model::solveDirichlet) .def("getBEEngine", &Model::getBEEngine, py::return_value_policy::reference); } void wrapModelFactory(py::module& mod) { py::class_(mod, "ModelFactory") .def_static("createModel", &ModelFactory::createModel, "model_type"_a, "system_size"_a, "discretization"_a); } void wrapModel(py::module& mod) { wrapBEEngine(mod); wrapModelClass(mod); wrapModelFactory(mod); wrapFunctionals(mod); } } // namespace wrap __END_TAMAAS__