diff --git a/doc/sphinx/source/model.rst b/doc/sphinx/source/model.rst index 4be5543..dc2fc75 100644 --- a/doc/sphinx/source/model.rst +++ b/doc/sphinx/source/model.rst @@ -1,103 +1,103 @@ Model and integral operators ============================ The class :cpp:class:`tamaas::Model` (and its counterpart :py:class:`Model `) is both a central class in Tamaas and one of the simplest. It mostly serves as holder for material properties, fields and integral operators, and apart from a linear elastic behavior does not perform any computation on its own. Model types ----------- :cpp:class:`tamaas::Model` has a concrete subclass :cpp:class:`tamaas::ModelTemplate` which implements the model function for a given :cpp:type:`tamaas::model_type`: :cpp:enumerator:`tamaas::basic_2d` Model type used in normal frictionless contact: traction and displacement are 2D fields with only one component. :cpp:enumerator:`tamaas::surface_2d` Model type used in frictional contact: traction and displacement are 2D fields with three components. :cpp:enumerator:`tamaas::volume_2d` Model type used in elastoplastic contact: tractions are the same as with :cpp:enumerator:`tamaas::surface_2d` but the displacement is a 3D field. The enumeration values suffixed with ``_1d`` are the one dimensional (line contact) counterparts of the above model types. The domain physical dimension and number of components are encoded in the class :cpp:class:`tamaas::model_type_traits`. Model creation and basic functionality -------------------------------------- The instanciation of a :py:class:`Model ` is done with the :py:class:`ModelFactory ` class and its :py:func:`createModel ` function:: physical_size = [1., 1.] discretization = [512, 512] model = tm.ModelFactory.createModel(tm.model_type.basic_2d, physical_size, discretization) .. warning:: For models of type ``volume_*d``, the first component of the ``physical_size`` and ``discretization`` arrays corresponds to the depth dimension (:math:`z` in most cases). For example:: tm.ModelFactory.createModel(tm.model_type.basic_2d, [0.3, 1, 1], [64, 81, 81]) creates a model of depth 0.3 and surface size 1\ :superscript:`2`, while the number of points is 64 in depth and 81\ :sup:`2` on the surface. This is done for data contiguity reasons, as we do discrete Fourier transforms in the horizontal plane. .. note:: If ran in an MPI context, the method :py:meth:`createModel ` expects the *global* system sizes and discretization of the model. The properties ``E`` and ``nu`` can be used to set the Young's modulus and Poisson ratio respectively:: model.E = 1 model.nu = 0.3 Fields can be easlily accessed with the ``[]`` operator, similar to Python's dictionaries:: surface_traction = model['traction'] To know what fields are available, you can call the :py:func:`getFields ` method. You can add new fields to a model object with :py:func:`registerField `, which is convenient for dumping. A model can also be used to compute stresses from a strain field:: import numpy as np - strain = np.zeros(model.getDiscretization() + [6]) # Mandel--Voigt notation + strain = np.zeros(model.shape + [6]) # Mandel--Voigt notation stress = np.zeros_like(strain) model.applyElasticity(stress, strain) Model dumpers ------------- The submodule `tamaas.dumpers` contains a number of classes to save model data into different formats: :py:class:`UVWDumper ` Dumps a model to `VTK `_ format. Requires the `UVW `_ python package which you can install with pip:: pip install uvw This dumper is made for visualization with VTK based software like `Paraview `_. :py:class:`NumpyDumper ` Dumps a model to a compressed Numpy file. :py:class:`H5Dumper ` Dumps a model to a compressed `HDF5 `_ file. Requires the `h5py `_ package. The dumpers are initialzed with a basename and the fields that you wish to write to file (optionally you can set ``all_fields`` to ``True`` to dump all fields in the model). By default, each write operation creates a new file in a separate directory (e.g. :py:class:`UVWDumper ` creates a ``paraview`` directory). To write to a specific file you can use the `dump_to_file` method. Here is a usage example:: from tamaas.dumpers import UVWDumper, H5Dumper # Create dumper uvw_dumper = UVWDumper('rough_contact_example', 'stress', 'plastic_strain') # Dump model uvw_dumper << model # Or alternatively model.addDumper(H5Dumper('rough_contact_archive', all_fields=True)) model.addDumper(uvw_dumper) model.dump() The last ``model.dump()`` call will call both dumpers. The resulting files will have the following hierachy:: ./paraview/rough_contact_example_0000.vtr ./paraview/rough_contact_example_0001.vtr ./hdf5/rough_contact_archive_0000.h5 .. note:: Currently, only :py:class:`H5Dumper ` supports parallel output with MPI. diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py index 2d29c6e..b45baec 100644 --- a/python/tamaas/dumpers/__init__.py +++ b/python/tamaas/dumpers/__init__.py @@ -1,342 +1,342 @@ # -*- mode:python; coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 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 . """ Dumpers for the class tamaas.Model """ from __future__ import print_function from sys import stderr from os import makedirs import os.path import numpy as np from .. import ModelDumper, model_type, mpi, type_traits from ._helper import step_dump, directory_dump, local_slice, _is_surface_field def _get_attributes(model): "Get model attributes" return { 'model_type': str(model.type), 'system_size': model.system_size, 'discretization': model.global_shape, } class FieldDumper(ModelDumper): """Abstract dumper for python classes using fields""" postfix = '' extension = '' name_format = "{basename}{postfix}.{extension}" def __init__(self, basename, *fields, **kwargs): """Construct with desired fields""" super(FieldDumper, self).__init__() self.basename = basename self.fields = list(fields) self.all_fields = kwargs.get('all_fields', False) def add_field(self, field): """Add another field to the dump""" if field not in self.fields: self.fields.append(field) def dump_to_file(self, file_descriptor, model): """Dump to a file (name or handle)""" def get_fields(self, model): """Get the desired fields""" if not self.all_fields: requested_fields = self.fields else: - requested_fields = model.getFields() + requested_fields = list(model) - return {field: model.getField(field) for field in requested_fields} + return {field: model[field] for field in requested_fields} def dump(self, model): "Dump model" self.dump_to_file(self.file_path, model) @property def file_path(self): """Get the default filename""" return self.name_format.format(basename=self.basename, postfix=self.postfix, extension=self.extension) @directory_dump('numpys') @step_dump class NumpyDumper(FieldDumper): """Dumper to compressed numpy files""" extension = 'npz' def dump_to_file(self, file_descriptor, model): """Saving to compressed multi-field Numpy format""" if mpi.size() > 1: raise RuntimeError("NumpyDumper does not function " "at all in parallel") np.savez_compressed(file_descriptor, attrs=_get_attributes(model), **self.get_fields(model)) try: import h5py @directory_dump('hdf5') @step_dump class H5Dumper(FieldDumper): """Dumper to HDF5 file format""" extension = 'h5' def dump_to_file(self, file_descriptor, model): """Saving to HDF5 with metadata about the model""" # Setup for MPI if not h5py.get_config().mpi and mpi.size() > 1: raise RuntimeError("HDF5 does not have MPI support") if mpi.size() > 1: from mpi4py import MPI # noqa mpi_args = dict(driver='mpio', comm=MPI.COMM_WORLD) comp_args = {} # compression does not work in parallel else: mpi_args = {} comp_args = dict(compression='gzip', compression_opts=7) with h5py.File(file_descriptor, 'w', **mpi_args) as handle: # Writing data for name, field in self.get_fields(model).items(): shape = list(field.shape) if mpi.size() > 1: xdim = 0 if _is_surface_field(field, model) else 1 shape[xdim] = MPI.COMM_WORLD.allreduce(shape[xdim]) dset = handle.create_dataset(name, shape, field.dtype, **comp_args) dset[local_slice(field, model)] = field # Writing metadata for name, attr in _get_attributes(model).items(): handle.attrs[name] = attr except ImportError: pass try: import uvw # noqa import uvw.parallel @directory_dump('paraview') @step_dump class UVWDumper(FieldDumper): """Dumper to VTK files for elasto-plastic calculations""" extension = 'vtr' forbidden_fields = ['traction', 'gap'] def dump_to_file(self, file_descriptor, model): """Dump displacements, plastic deformations and stresses""" if mpi.size() > 1: raise RuntimeError("UVWDumper does not function " "properly in parallel") - bdim = len(model.getBoundaryDiscretization()) + bdim = len(model.boundary_shape) # Local MPI size - lsize = model.getDiscretization() - gsize = mpi.global_shape(model.getBoundaryDiscretization()) + lsize = model.shape + gsize = mpi.global_shape(model.boundary_shape) gshape = gsize if len(lsize) > bdim: - gshape = [model.getDiscretization()[0]] + gshape + gshape = [model.shape[0]] + gshape # Space coordinates coordinates = [np.linspace(0, L, N, endpoint=False) - for L, N in zip(model.getSystemSize(), gshape)] + for L, N in zip(model.system_size, gshape)] # If model has subsurfce domain, z-coordinate is always first dimension_indices = np.arange(bdim) if len(lsize) > bdim: dimension_indices += 1 dimension_indices = np.concatenate((dimension_indices, [0])) coordinates[0] = \ - np.linspace(0, model.getSystemSize()[0], gshape[0]) + np.linspace(0, model.system_size[0], gshape[0]) offset = np.zeros_like(dimension_indices) offset[0] = mpi.local_offset(gsize) rectgrid = uvw.RectilinearGrid if mpi.size() == 1 \ else uvw.parallel.PRectilinearGrid # Creating rectilinear grid with correct order for components coordlist = [coordinates[i][o:o+lsize[i]] for i, o in zip(dimension_indices, offset)] grid = rectgrid( file_descriptor, coordlist, compression=True, offsets=offset, ) # Iterator over fields we want to dump fields_it = filter(lambda t: t[0] not in self.forbidden_fields, self.get_fields(model).items()) # We make fields periodic for visualization for name, field in fields_it: array = uvw.DataArray(field, dimension_indices, name) grid.addPointData(array) grid.write() @directory_dump('paraview') class UVWGroupDumper(FieldDumper): "Dumper to ParaViewData files" extension = 'pvd' def __init__(self, basename, *fields, **kwargs): """Construct with desired fields""" super(UVWGroupDumper, self).__init__(basename, *fields, **kwargs) subdir = os.path.join('paraview', basename + '-VTR') if not os.path.exists(subdir): makedirs(subdir) self.uvw_dumper = UVWDumper( os.path.join(basename + '-VTR', basename), *fields, **kwargs ) self.group = uvw.ParaViewData(self.file_path, compression=True) def dump_to_file(self, file_descriptor, model): self.group.addFile( self.uvw_dumper.file_path.replace('paraview/', ''), timestep=self.uvw_dumper.count ) self.group.write() self.uvw_dumper.dump(model) except ImportError as error: print(error, file=stderr) try: from netCDF4 import Dataset @directory_dump('netcdf') class NetCDFDumper(FieldDumper): """Dumper to netCDF4 files""" extension = "nc" boundary_fields = ['traction', 'gap'] def _file_setup(self, grp, model): grp.createDimension('frame', None) # Local dimensions - model_dim = len(model.getDiscretization()) + model_dim = len(model.shape) voigt_dim = type_traits[model.type].voigt self._vec = grp.createDimension('spatial', model_dim) self._tens = grp.createDimension('Voigt', voigt_dim) - self.model_info = model.getDiscretization(), model.type + self.model_info = model.shape, model.type # Create boundary dimensions for label, size, length in zip( "xy", - model.getBoundaryDiscretization(), - model.getBoundarySystemSize() + model.boundary_shape, + model.boundary_system_size ): grp.createDimension(label, size) coord = grp.createVariable(label, 'f8', (label,)) coord[:] = np.linspace(0, length, size, endpoint=False) self._create_variables( grp, model, lambda f: _is_surface_field(f[1], model), - model.getBoundaryDiscretization(), "xy" + model.boundary_shape, "xy" ) # Create volume dimension if model.type in {model_type.volume_1d, model_type.volume_2d}: - size = model.getDiscretization()[0] + size = model.shape[0] grp.createDimension("z", size) coord = grp.createVariable("z", 'f8', ("z",)) - coord[:] = np.linspace(0, model.getSystemSize()[0], size) + coord[:] = np.linspace(0, model.system_size[0], size) self._create_variables( grp, model, lambda f: not _is_surface_field(f[1], model), - model.getDiscretization(), "zxy" + model.shape, "zxy" ) self.has_setup = True def dump_to_file(self, file_descriptor, model): if mpi.size() > 1: raise RuntimeError("NetCDFDumper does not function " "properly in parallel") mode = 'a' if os.path.isfile(file_descriptor) \ and getattr(self, 'has_setup', False) else 'w' with Dataset(file_descriptor, mode, format='NETCDF4_CLASSIC') as rootgrp: if rootgrp.dimensions == {}: self._file_setup(rootgrp, model) - if self.model_info != (model.getDiscretization(), model.type): + if self.model_info != (model.shape, model.type): raise Exception("Unexpected model {}".format(model)) self._dump_generic(rootgrp, model) def _create_variables(self, grp, model, predicate, shape, dimensions): field_dim = len(shape) fields = list(filter(predicate, self.get_fields(model).items())) dim_labels = list(dimensions[:field_dim]) for label, data in fields: local_dim = [] # If we have an extra component if data.ndim > field_dim: if data.shape[-1] == self._tens.size: local_dim = [self._tens.name] elif data.shape[-1] == self._vec.size: local_dim = [self._vec.name] grp.createVariable(label, 'f8', ['frame'] + dim_labels + local_dim, zlib=True) def _dump_generic(self, grp, model): fields = self.get_fields(model).items() new_frame = grp.dimensions['frame'].size for label, data in fields: var = grp[label] slice_in_global = (new_frame,) + local_slice(data, model) var[slice_in_global] = np.array(data, dtype=np.double) except ImportError: pass diff --git a/python/tamaas/dumpers/_helper.py b/python/tamaas/dumpers/_helper.py index 66bbbcb..08795a7 100644 --- a/python/tamaas/dumpers/_helper.py +++ b/python/tamaas/dumpers/_helper.py @@ -1,141 +1,141 @@ # -*- mode:python; coding: utf-8 -*- # @file # @section LICENSE # # Copyright (©) 2016-2020 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 . """ Helper functions for dumpers """ from functools import wraps from collections import defaultdict import os import numpy as np from .. import model_type, type_traits, mpi __all__ = ["step_dump", "directory_dump"] def _is_surface_field(field, model): - bn = model.getBoundaryDiscretization() + bn = model.boundary_shape return list(field.shape[:len(bn)]) == bn def local_slice(field, model): - n = model.getDiscretization() - bn = model.getBoundaryDiscretization() + n = model.shape + bn = model.boundary_shape gshape = mpi.global_shape(bn) offsets = np.zeros_like(gshape) offsets[0] = mpi.local_offset(gshape) if not _is_surface_field(field, model) and len(n) > len(bn): gshape = [n[0]] + gshape offsets = np.concatenate(([0], offsets)) return tuple((slice(offset, offset+size, None) for offset, size in zip(offsets, field.shape))) def step_dump(cls): """ Decorator for dumper with counter for steps """ orig_init = cls.__init__ orig_dump = cls.dump @wraps(cls.__init__) def __init__(obj, *args, **kwargs): orig_init(obj, *args, **kwargs) obj.count = 0 def postfix(obj): return "_{:04d}".format(obj.count) @wraps(cls.dump) def dump(obj, *args, **kwargs): orig_dump(obj, *args, **kwargs) obj.count += 1 cls.__init__ = __init__ cls.dump = dump cls.postfix = property(postfix) return cls def directory_dump(directory=""): "Decorator for dumper in a directory" def actual_decorator(cls): orig_dump = cls.dump orig_filepath = cls.file_path.fget @wraps(cls.dump) def dump(obj, *args, **kwargs): if not os.path.exists(directory) and mpi.rank() == 0: os.mkdir(directory) if not os.path.isdir(directory) and mpi.rank() == 0: raise Exception('{} is not a directory'.format(directory)) orig_dump(obj, *args, **kwargs) @wraps(cls.file_path.fget) def file_path(obj): return os.path.join(directory, orig_filepath(obj)) cls.dump = dump cls.file_path = property(file_path) return cls return actual_decorator def hdf5toVTK(inpath, outname): "Convert HDF5 dump of a model to VTK" import h5py as h5 # noqa from .. import ModelFactory # noqa from . import UVWDumper # noqa type_translate = { 'model_type.basic_1d': model_type.basic_1d, 'model_type.basic_2d': model_type.basic_2d, 'model_type.surface_1d': model_type.surface_1d, 'model_type.surface_2d': model_type.surface_2d, 'model_type.volume_1d': model_type.volume_1d, 'model_type.volume_2d': model_type.volume_2d, } with h5.File(inpath, 'r') as h5file: model_t = h5file.attrs['model_type'] system_size = list(h5file.attrs['system_size']) n = list(h5file.attrs['discretization']) model = ModelFactory.createModel(type_translate[model_t], system_size, n) fields = [] for field in h5file: model.registerField(field, h5file[field][:]) fields.append(field) uvw_dumper = UVWDumper(outname, *fields) uvw_dumper << model diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index d3a1914..cc8dfc8 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,372 +1,377 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 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 . * */ /* -------------------------------------------------------------------------- */ #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 /* -------------------------------------------------------------------------- */ 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::wrap::PyFunctional> func(mod, "Functional"); 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 + // TODO [legacy] remove wrapper code .def("setParameters", &functional::AdhesionFunctional::setParameters); py::class_( mod, "ExponentialAdhesionFunctional", adh) .def(py::init&>(), "surface"_a); py::class_( mod, "MaugisAdhesionFunctional", adh) .def(py::init&>(), "surface"_a); py::class_( mod, "SquaredExponentialAdhesionFunctional", adh) .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") // TODO [legacy] remove .def("apply", [](IntegralOperator& op, numpy input, numpy output) { auto in = instanciateFromNumpy(input); auto out = instanciateFromNumpy(output); op.apply(*in, *out); }) .def("getModel", &IntegralOperator::getModel, py::return_value_policy::reference) .def("getKind", &IntegralOperator::getKind) .def("getType", &IntegralOperator::getType) .def("__call__", [](IntegralOperator& op, numpy input, numpy output) { auto in = instanciateFromNumpy(input); auto out = instanciateFromNumpy(output); op.apply(*in, *out); }) .def("updateFromModel", &IntegralOperator::updateFromModel) .def_property_readonly("model", &IntegralOperator::getModel) .def_property_readonly("type", &IntegralOperator::getType); py::enum_(mod, "integration_method") .value("linear", integration_method::linear) .value("cutoff", integration_method::cutoff); } /// Wrap BEEngine classes 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) .def("registerNeumann", &BEEngine::registerNeumann) - .def("registerDirichlet", &BEEngine::registerDirichlet); + .def("registerDirichlet", &BEEngine::registerDirichlet) + // TODO [legacy] remove + .def("getModel", &BEEngine::getModel, 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", [](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_(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(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_property_readonly("type", &Model::getType) .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a) .def_property("E", &Model::getYoungModulus, &Model::setYoungModulus) .def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio) .def("getHertzModulus", &Model::getHertzModulus) .def("getYoungModulus", &Model::getYoungModulus) .def("getShearModulus", &Model::getShearModulus) .def("getPoissonRatio", &Model::getPoissonRatio) .def("getTraction", (GridBase & (Model::*)()) & Model::getTraction, py::return_value_policy::reference_internal) .def("getDisplacement", (GridBase & (Model::*)()) & Model::getDisplacement, py::return_value_policy::reference_internal) .def("getSystemSize", &Model::getSystemSize) .def("getDiscretization", &Model::getDiscretization) .def("getBoundarySystemSize", &Model::getBoundarySystemSize) .def("getBoundaryDiscretization", &Model::getBoundaryDiscretization) .def("solveNeumann", &Model::solveNeumann) .def("solveDirichlet", &Model::solveDirichlet) .def("dump", &Model::dump) .def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>()) .def("setDumper", [](Model&, const ModelDumper&) { TAMAAS_EXCEPTION("setDumper() is not a member of Model; use " "addDumper() instead"); }) .def("getBEEngine", &Model::getBEEngine, py::return_value_policy::reference_internal) .def("getIntegralOperator", &Model::getIntegralOperator, "operator_name"_a, py::return_value_policy::reference_internal) .def("registerField", [](Model& m, std::string name, numpy field) { auto f = instanciateFromNumpy(field); m.registerField(name, std::move(f)); }, "field_name"_a, "field"_a, py::keep_alive<1, 3>()) .def("getField", (GridBase & (Model::*)(const std::string&)) & Model::getField, "field_name"_a, py::return_value_policy::reference_internal) .def("getFields", &Model::getFields) .def("applyElasticity", [](Model& model, numpy stress, numpy strain) { auto out = instanciateFromNumpy(stress); auto in = instanciateFromNumpy(strain); model.applyElasticity(*out, *in); }) // Python magic functions .def("__repr__", [](const Model& m) { std::stringstream ss; ss << m; return ss.str(); }) .def("__getitem__", [](const Model& m, std::string key) -> const GridBase& { try { return m.getField(key); } catch (std::out_of_range&) { throw py::key_error(key); } }, py::return_value_policy::reference_internal) .def("__setitem__", [](Model& m, std::string name, numpy field) { auto f = instanciateFromNumpy(field); m.registerField(name, std::move(f)); }, py::keep_alive<1, 3>()) .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>()) .def("__iter__", [](const Model& m) { const auto& fields = m.getFieldsMap(); return py::make_key_iterator(fields.cbegin(), fields.cend()); }, py::keep_alive<0, 1>()) .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) .def_property_readonly("global_shape", &Model::getGlobalDiscretization) .def_property_readonly("boundary_shape", &Model::getBoundaryDiscretization) - .def_property_readonly("system_size", &Model::getSystemSize); + .def_property_readonly("system_size", &Model::getSystemSize) + .def_property_readonly("boundary_system_size", + &Model::getBoundarySystemSize); py::class_>( mod, "ModelDumper") .def(py::init<>()) .def("dump", &ModelDumper::dump, "model"_a) .def("__lshift__", [](ModelDumper& dumper, Model& model) { dumper << model; }); } /// Wrap factory for models void wrapModelFactory(py::module& mod) { py::class_(mod, "ModelFactory") .def_static("createModel", &ModelFactory::createModel, "model_type"_a, "system_size"_a, "global_discretization"_a) .def_static("createResidual", &ModelFactory::createResidual, "model"_a, "sigma_y"_a, "hardening"_a = 0.) .def_static("registerVolumeOperators", &ModelFactory::registerVolumeOperators, "model"_a); } /// 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/surface.cpp b/python/wrap/surface.cpp index fb95ed6..cbce617 100644 --- a/python/wrap/surface.cpp +++ b/python/wrap/surface.cpp @@ -1,216 +1,217 @@ /** * @file * @section LICENSE * * Copyright (©) 2016-2020 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 . * */ /* -------------------------------------------------------------------------- */ #include "isopowerlaw.hh" #include "regularized_powerlaw.hh" #include "surface_generator_filter.hh" #include "surface_generator_filter_fft.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 { 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()) .def(py::init<>()) .def("computeFilter", (void (Filter::*)(GridHermitian&) const) & Filter::computeFilter); } /* -------------------------------------------------------------------------- */ template void wrapIsopowerlaw(py::module& mod) { std::string name = makeDimensionName("Isopowerlaw", dim); py::class_, Filter, std::shared_ptr>>( mod, name.c_str()) .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::setHurst, "Hurst exponent") .def("rmsHeights", &Isopowerlaw::rmsHeights, "Theoretical RMS of heights") .def("moments", &Isopowerlaw::moments, "Theoretical first 3 moments of spectrum") .def("alpha", &Isopowerlaw::alpha, "Nayak's bandwidth parameter") .def("rmsSlopes", &Isopowerlaw::rmsSlopes, "Theoretical RMS of slopes") - // legacy wrapper code + // TODO [legacy] remove wrapper code .def("getQ0", &Isopowerlaw::getQ0, py::return_value_policy::reference) .def("getQ1", &Isopowerlaw::getQ1, py::return_value_policy::reference) .def("getQ2", &Isopowerlaw::getQ2, py::return_value_policy::reference) .def("setQ0", &Isopowerlaw::setQ0, py::return_value_policy::reference) .def("setQ1", &Isopowerlaw::setQ1, py::return_value_policy::reference) .def("setQ2", &Isopowerlaw::setQ2, py::return_value_policy::reference) .def("setHurst", &Isopowerlaw::setHurst) .def("getHurst", &Isopowerlaw::getHurst, py::return_value_policy::reference); name = makeDimensionName("RegularizedPowerlaw", dim); py::class_, Filter, std::shared_ptr>>(mod, name.c_str()) .def(py::init<>()) .def_property("q1", &RegularizedPowerlaw::getQ1, &RegularizedPowerlaw::setQ1, "Long wavelength cutoff") .def_property("q2", &RegularizedPowerlaw::getQ2, &RegularizedPowerlaw::setQ2, "Short wavelength cutoff") .def_property("hurst", &RegularizedPowerlaw::getHurst, &RegularizedPowerlaw::setHurst, "Hurst exponent"); } /* -------------------------------------------------------------------------- */ template void wrapSurfaceGenerators(py::module& mod) { std::string generator_name = makeDimensionName("SurfaceGenerator", dim); py::class_>(mod, generator_name.c_str()) .def("buildSurface", &SurfaceGenerator::buildSurface, py::return_value_policy::reference_internal) .def("setSizes", py::overload_cast>( &SurfaceGenerator::setSizes), "global_surface_shape"_a) .def_property("random_seed", &SurfaceGenerator::getRandomSeed, &SurfaceGenerator::setRandomSeed, "Random generator seed") .def_property("shape", &SurfaceGenerator::getSizes, py::overload_cast>( &SurfaceGenerator::setSizes), "Global shape of surfaces"); std::string filter_name = makeDimensionName("SurfaceGeneratorFilter", dim); py::class_, SurfaceGenerator>( mod, filter_name.c_str()) .def(py::init<>(), "Default constructor") .def(py::init>(), "Initialize with global surface shape") + // TODO [legacy] remove .def("setFilter", &SurfaceGeneratorFilter::setFilter, "Set PSD filter", "filter"_a) .def("setSpectrum", &SurfaceGeneratorFilter::setSpectrum, "Set PSD spectrum", "filter"_a) - // legacy wrapper code .def("setRandomSeed", &SurfaceGeneratorFilter::setRandomSeed, "[[Deprecated: use property]]") + .def_property("spectrum", &SurfaceGeneratorFilter::getSpectrum, &SurfaceGeneratorFilter::setSpectrum, "Power spectrum object"); std::string random_name = makeDimensionName("SurfaceGeneratorRandomPhase", dim); py::class_, SurfaceGeneratorFilter>( mod, random_name.c_str()) .def(py::init<>(), "Default constructor") .def(py::init>(), "Initialize with global surface shape"); } /* -------------------------------------------------------------------------- */ -// Legacy wrap +// TODO [legacy] remove wrap void wrapSurfaceGeneratorFilterFFT(py::module& mod) { py::class_(mod, "SurfaceGeneratorFilterFFT") .def(py::init<>()) .def("getQ0", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ0()); }) .def("getQ1", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ1()); }) .def("getQ2", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getQ2()); }) .def("getHurst", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getHurst()); }) .def("analyticalRMS", &SurfaceGeneratorFilterFFT::analyticalRMS) .def("buildSurface", &SurfaceGeneratorFilterFFT::buildSurface, py::return_value_policy::reference_internal) .def("getGridSize", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getGridSize()); }) .def("getRandomSeed", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getRandomSeed()); }) .def("getRMS", [](SurfaceGeneratorFilterFFT& f) { return smart_pointer(&f.getRMS()); }) .def("Init", &SurfaceGeneratorFilterFFT::Init); } /* -------------------------------------------------------------------------- */ void wrapSurface(py::module& mod) { wrapFilter<1>(mod); wrapFilter<2>(mod); wrapIsopowerlaw<1>(mod); wrapIsopowerlaw<2>(mod); wrapSurfaceGenerators<1>(mod); wrapSurfaceGenerators<2>(mod); // legacy wrap #if defined(TAMAAS_LEGACY_BEM) wrapSurfaceGeneratorFilterFFT(mod); #endif } } // namespace wrap /* -------------------------------------------------------------------------- */ } // namespace tamaas