diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f72d0b..f7c8f3c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,386 +1,394 @@ # Changelog All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html) for final versions and [PEP440](https://www.python.org/dev/peps/pep-0440/) in case intermediate versions need to be released (e.g. development version `2.2.3.dev1` or release candidates `2.2.3rc1`), or individual commits are packaged. ## Unreleased ### Added -- Added a `Material` class tree to decouple material behavior from residuals +- Added a `Material` class tree to decouple material behavior from residuals. + The `Material` class can be derived in Python for custom material laws. +- `Model` can set the integration method for volume operators with the + `setIntegrationMethod` method. ### Changed -- Removed `model_type` template for `Residual` -- Simplified interface for `Residual` +- Removed `model_type` template for `Residual`: they are only defined on + the `model_type::volume_2d` type. +- Simplified interface for `Residual`: it now takes a `Material` object at + construction which handles the computation of eigenstresses involved in the + residual equation ### Deprecated - `getStress()`, `setIntegrationMethod()` and `getVector()` in `Residual` are deprecated, use `model["stress"]`, `model.setIntegrationMethod()` and - `residual.vector` respectively instead + `residual.vector` respectively instead. +- `ModelFactor.createResidual` is deprecated. Instead, instanciate a `Material` + object and use the `Residual` constructor to create the residual object. ## v2.6.0 -- 2022-12-23 ### Added - Added `matvec` interface for the following operators: - `Hooke` - `Kelvin` - `Westergaard` - Documented error handling with MPI - Added `span` class - `FieldContainer` class that defines an interface to access fields registered with an object - Operators now expose their fields, which can be modified from Python - Added `HookeField` operator for heterogeneous elastic properties - Documented proper unit choice - Added `Statistics::computeFDRMSSlope()`, finite difference approximation of RMS of surface slopes - Kato saturated uses L2 norm on residual displacement instead of L∞ - Request system to avoid erasing model data on restart - Boundary field names are dumped ### Changed - Resizing a grid no longer sets values to zero - Wheels built for `manylinux2014` instead of `manylinux2010` - Cleaned up `Array` interface - Shared-memory compile options disabled by default - Unknown build variables raise a warning - Removed conditional branches from PKR - NetCDF dumps in `NETCDF4` format instead of `NETCDF4_CLASSIC` - VTK no longer writes FieldData, this was causing issues reading converted files ### Fixed - Fixed flood fill algorithm to properly account for PBCs ## v2.5.0 -- 2022-07-05 ### Added - Added `tamaas.utils.publications` to print a list of publications relevant to the parts of Tamaas used in a script. Call at the end of a script for an exhaustive list. - Added `tamaas.mpi.gather` and `tamaas.mpi.scatter` to help handling of 2D data in MPI contexts. - Added `getBoundaryFields()/boundary_fields` function/property to `Model`. - Added Python constructor to `Model`, which is an alias to `ModelFactory.createModel`. - Added convenience `tamaas.dumpers._helper.netCDFtoParaview` function to convert model sequence to Paraview data file. - Added dump of surface fields as `FieldData` in VTK files - Added `scipy.sparse.linalg.LinearOperator` interface to `Westergaard`. A `LinearOperator` can be created by calling `scipy.sparse.linalg.aslinearoperator` on the operator object. - Allowed sub-classing of `ContactSolver` in Python. - Exposed elastic contact functionals in Python to help with custom solver implementation. - Added an example of custom solver with penalty method using Scipy. - Added `graphArea` function in statistics to compute area of function graph - Added `tamaas.utils.radial_average` to radialy average a 2D field - Added option to not make a dump directory (e.g. `netcdf`, `hdf5` directories for dump files). This allows to specify an absolute path for dump files, and is useful if simulation workflow is managed with file targets, as do systems like GNU Make and Snakemake. ### Changed - Filter for boundary fields in dumper should be smarter. - `read` functions in dumpers are now `@classmethod` - Changed sign of plastic correction in KatoSaturated solver to match the EPICSolver sign. - Plastic correction of KatoSaturated is registered in the model as `KatoSaturated::residual_displacement`. - Changed to modern Python packaging (PEP517) with `setup.cfg` and `pyproject.toml`. A `setup.py` script is still provided for editable installations, see [gh#7953](https://github.com/pypa/pip/issues/7953). Once setuptools has stable support for [PEP621](https://peps.python.org/pep-0621/), `setup.cfg` will be merged into `pyproject.toml`. Note that source distributions that can be built with [`build`](https://pypi.org/project/build/) are not true source distributions, since they contain the compiled Python module of Tamaas. - Updated copyright to reflect Tamaas' development history since leaving EPFL ### Fixed - Fixed `tamaas.dumpers.NetCDFDumper` to dump in MPI context. ### Deprecated - SCons versions < 3 and Python < 3.6 are explicitely no longer supported. ## v2.4.0 -- 2022-03-22 ### Added - Added a `tamaas.utils` module. - Added `tamaas.utils.load_path` generator, which yields model objects for a sequence of applied loads. - Added `tamaas.utils.seeded_surfaces` generator, which yields surfaces for a sequence of random seeds. - Added `tamaas.utils.hertz_surface` which generates a parabolic (Hertzian) surface. - Added `tamaas.utils.log_context` context manager which locally sets the log level for Tamaas' logger. - Added `tamaas.compute.from_voigt` to convert Voigt/Mendel fields to dense tensor fields. - Added a deep-copy function to `ModelFactory` to copy `Model` objects. Use `model_copy = copy.deepcopy(model)` in Python to make a copy. Currently only copies registered fields, same as dumpers/readers. - Added a `read_sequence` method for dumpers to read model frames. - Automatic draft release on Zenodo. ### Changed - `*_ROOT` variables for vendored dependencies GoogleTest and pybind11 now default to empty strings, so that the vendored trees in `third-party/` are not selected by default. This is so system packages are given priority. Vendored dependency submodules will eventually be deprecated. ### Fixed - Fixed an issue with `scons -c` when GTest was missing. ## v2.3.1 -- 2021-11-08 ### Added - Now using `clang-format`, `clang-tidy` and `flake8` for linting - Pre-built Tamaas images can be pulled with `docker pull registry.gitlab.com/tamaas/tamaas` - Added a `--version` flag to the `tamaas` command ### Changed - The root `Dockerfile` now compiles Tamaas, so using Tamaas in Docker is easier. - The model attributes dumped to Numpy files are now written in a JSON-formatted string to avoid unsafe loading/unpickling of objects. - Removed the `build_doc` build option: now the doc targets are automatically added if the dependencies are met, and built if `scons doc` is called. - Removed the `use_googletest` build option: if tests are built and gtest is present, the corresponding tests will be built ### Fixed - The command `tamaas plot` gracefully fails if Matplotlib is missing - Better heuristic to guess which fields are defined on boundary in dumpers (still not perfect and may give false negatives) ## v2.3.0 -- 2021-06-15 ### Added - Added `read()` method to dumpers to create a model from a dump file - `getClusters()` can be called in MPI contact with partial contact maps - Added a JSON encoder class for models and a JSON dumper - CUDA compatibility is re-established, but has not been tested - Docstrings in the Python bindings for many classes/methods ### Changed - Tamaas version numbers are now managed by [versioneer](https://github.com/python-versioneer/python-versioneer). This means that Git tags prefixed with `v` (e.g. `v2.2.3`) carry meaning and determine the version. When no tag is set, versioneer uses the last tag, specifies the commit short hash and the distance to the last tag (e.g. `2.2.2+33.ge314b0e`). This version string is used in the compiled library, the `setup.py` script and the `__version__` variable in the python module. - Tamaas migrated to [GitLab](https://gitlab.com/tamaas/tamaas) - Continuous delivery has been implemented: - the `master` branch will now automatically build and publish Python wheels to `https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple`. These "nightly" builds can be installed with: pip install \ --extra-index-url https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple \ tamaas - version tags pushed to `master` will automatically publish the wheels to [PyPI](https://pypi.org/project/tamaas/) ### Deprecated - The `finalize()` function is now deprecated, since it is automatically called when the process terminates - Python versions 3.5 and below are not supported anymore ### Fixed - Fixed a host of dump read/write issues when model type was not `volume_*d`. Dumper tests are now streamlined and systematic. - Fixed a bug where `Model::solveDirichlet` would not compute correctly - Fixed a bug where `Statistics::contact` would not normalize by the global number of surface points ## v2.2.2 -- 2021-04-02 ### Added - Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try executing `tamaas surface -h` from the command-line! ### Changed - `CXXFLAGS` are now passed to the linker - Added this changelog - Using absolute paths for environmental variables when running `scons test` - Reorganized documentation layout - Gave the build system a facelift (docs are now generated directly with SCons instead of a Makefile) ### Deprecated - Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with a Python 2 wheel. - The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command ### Fixed - `UVWDumper` no longer imports `mpi4py` in sequential - Compiling with different Thrust/FFTW backends ## v2.2.1 -- 2021-03-02 ### Added - Output registered fields and dumpers in `print(model)` - Added `operator[]` to the C++ model class (for fields) - Added `traction` and `displacement` properties to Python model bindings - Added `operators` property to Python model bindings, which provides a dict-like access to registered operators - Added `shape` and `spectrum` to properties to Python surface generator bindings - Surface generator constructor accepts surface global shape as argument - Choice of FFTW thread model ### Changed - Tests use `/tmp` for temporary files - Updated dependency versions (Thrust, Pybind11) ### Deprecated - Most `get___()` and `set___()` in Python bindings have been deprecated. They will generate a `DeprecationWarning`. ### Removed - All legacy code ## v2.2.0 -- 2020-12-31 ### Added - More accurate function for computation of contact area - Function to compute deviatoric of tensor fields - MPI implementation - Convenience `hdf5toVTK` function - Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give shape information ### Changed - Preprocessor defined macros are prefixed with `TAMAAS_` - Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt` ### Fixed - Warning about deprecated constructors with recent GCC versions - Wrong computation of grid strides - Wrong computation of grid sizes in views ## v2.1.4 -- 2020-08-07 ### Added - Possibility to generate a static `libTamaas` - C++ implementation of DFSANE solver - Allowing compilation without OpenMP ### Changed - NetCDF dumper writes frames to a single file ### Fixed - Compatibility with SCons+Python 3 ## v2.1.3 -- 2020-07-27 ### Added - Version number to `TamaasInfo` ### Changed - Prepending root directory when generating archive ## v2.1.2 -- 2020-07-24 This release changes some core internals related to discrete Fourier transforms for future MPI support. ### Added - Caching `CXXFLAGS` in SCons build - SCons shortcut to create code archive - Test of the elastic-plastic contact solver - Paraview data dumper (`.pvd` files) - Compression for UVW dumper - `__contains__` and `__iter__` Python bindings of model - Warning message of possible overflow in Kelvin ### Changed - Simplified `tamaas_info.cpp`, particularly the diff part - Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are re-used as much as possible with different data with the same shape. This is in view of future MPI developments - Redirecting I/O streams in solve functions so they can be used from Python (e.g. in Jupyter notebooks) - Calling `initialize()` and `finalize()` is no longer necessary ### Fixed - Convergence issue with non-linear solvers - Memory error in volume potentials ## v2.1.1 -- 2020-04-22 ### Added - SCons shortcut to run tests ### Fixed - Correct `RPATH` for shared libraries - Issues with SCons commands introduced in v2.1.0 - Tests with Python 2.7 ## v2.1.0 -- 2020-04-17 ### Added - SCons shortcuts to build/install Tamaas and its components - Selection of integration method for Kelvin operator - Compilation option to remove the legacy part of Tamaas - NetCDF dumper ### Fixed - Link bug with clang - NaNs in Kato saturated solver ## v2.0.0 -- 2019-11-11 First public release. Contains relatively mature elastic-plastic contact code. diff --git a/doc/sphinx/source/contact.rst b/doc/sphinx/source/contact.rst index 9f1e773..d9b4d5d 100644 --- a/doc/sphinx/source/contact.rst +++ b/doc/sphinx/source/contact.rst @@ -1,175 +1,176 @@ Solving contact =============== The resolution of contact problems is done with classes that inherit from :cpp:class:`tamaas::ContactSolver`. These usually take as argument a :cpp:class:`tamaas::Model` object, a surface described by a :cpp:class:`tamaas::Grid` or a 2D numpy array, and a tolerance. We will see the specificities of the different solver objects below. Elastic contact --------------- The most common case is normal elastic contact, and is most efficiently solved with :cpp:class:`tamaas::PolonskyKeerRey`. The advantage of this class is that it combines two algorithms into one. By default, it considers that the contact pressure field is the unknown, and tries to minimize the complementary energy of the system under the constraint that the mean pressure should be equal to the value supplied by the user, for example:: # ... solver = tm.PolonskyKeerRey(model, surface, 1e-12) solver.solve(1e-2) Here the average pressure is ``1e-2``. The solver can also be instanciated by specifying the the constraint should be on the mean gap instead of mean pressure:: solver = tm.PolonskyKeerRey(model, surface, 1e-12, constraint_type=tm.PolonskyKeerRey.gap) solver.solve(1e-2) The contact problem is now solved for a mean gap of ``1e-2``. Note that the choice of constraint affects the performance of the algorithm. Computing solutions for loading sequence ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The module :py:mod:`tamaas.utils` defines a convenience function :py:func:`load_path `, which generates solution models for a sequence of loads. This allows lazy evalution and reduces boiler-plate:: from tamaas.utils import load_path loads = np.linspace(0.01, 0.1, 10) for model in load_path(solver, loads): ... # do some computation on model, e.g. compute contact clusters The function :py:func:`load_path ` accepts the following optional arguments: ``verbose`` Prints solver output (i.e. iteration, cost function and error) ``callback`` A function to call after each solve, before the next load step. Useful for dumping model in generator expressions. Contact with adhesion --------------------- The second algorithm hidden in :cpp:class:`tamaas::PolonskyKeerRey` considers the **gap** to be the unknown. The constaint on the mean value can be chosen as above:: solver = tm.PolonskyKeerRey(model, surface, 1e-12, primal_type=tm.PolonskyKeerRey.gap, constraint_type=tm.PolonskyKeerRey.gap) solver.solve(1e-2) The advantage of this formulation is to be able to solve adhesion problems (`Rey et al. `_). This is done by adding a term to the potential energy functional that the solver tries to minimize:: adhesion_params = { "rho": 2e-3, "surface_energy": 2e-5 } adhesion = tm.ExponentialAdhesionFunctional(surface) adhesion.setParameters(adhesion) solver.addFunctionalterm(adhesion) solver.solve(1e-2) Custom classes can be used in place of the example term here. One has to inherit from :cpp:class:`tamaas::Functional`:: import numpy class AdhesionPython(tm.Functional): """ Functional class that extends a C++ class and implements the virtual methods """ def __init__(self, rho, gamma): super().__init__(self) self.rho = rho self.gamma = gamma def computeF(self, gap, pressure): return -self.gamma * numpy.sum(np.exp(-gap / self.rho)) def computeGradF(self, gap, gradient): gradient += self.gamma * numpy.exp(-gap / self.rho) / self.rho This example is actually equivalent to :cpp:class:`tamaas::functional::ExponentialAdhesionFunctional`. Tangential contact ------------------ For frictional contact, several solver classes are available. Among them, :cpp:class:`tamaas::Condat` is able to solve a coupled normal/tangential contact problem regardless of the material properties. It however solves an associated version of the Coulomb friction law. In general, the Coulomb friction used in contact makes the problem ill-posed. Solving a tangential contact problem is not much different from normal contact:: mu = 0.3 # Friction coefficient solver = tm.Condat(model, surface, 1e-12, mu) solver.max_iter = 5000 # The default of 1000 may be too little solver.solve([1e-2, 0, 1e-2]) # 3D components of applied mean pressure Elasto-plastic contact ---------------------- For elastic-plastic contact, one needs three different solvers: an elastic contact solver like the ones described above, a non-linear solver and a coupling solver. The non-linear solvers available in Tamaas are implemented in python and inherit from the C++ class :cpp:class:`tamaas::EPSolver`. They make use of the non-linear solvers available in scipy: :py:class:`DFSANESolver ` Implements an interface to :py:func:`scipy.optimize.root` wiht the DFSANE method. :py:class:`NewtonKrylovSolver ` Implements an interface to :py:func:`scipy.optimize.newton_krylov`. These solvers require a residual vector to cancel, the creation of which is handled with :cpp:class:`tamaas::ModelFactory`. Then, an instance of :cpp:class:`tamaas::EPICSolver` is needed to couple the contact and non-linear solvers for the elastic-plastic contact problem:: import tamaas as tm from tamaas.nonlinear_solvers import DFSANESolver # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 51, 51] # Order: [z, x, y] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Setup for plasticity - residual = tm.ModelFactory.createResidual(model, - sigma_y=0.1 * model.E, - hardening=0.01 * model.E) + material = tm.materials.IsotropicHardening(model, + sigma_y=0.1 * model.E, + hardening=0.01 * model.E) + residual = tm.Residual(model, material) epsolver = DFSANESolver(residual) # ... csolver = tm.PolonskyKeerRey(model, surface, 1e-12) epic = tm.EPICSolver(csolver, epsolver, 1e-7, relaxation=0.3) epic.solve(1e-3) # or use an accelerated scheme (which can sometimes not converge) epic.acceleratedSolve(1e-3) By default, :cpp:func:`tamaas::EPICSolver::solve` uses a relaxed fixed point. It can be tricky to make it converge: you need to decrease the relaxation parameter passed as argument of the constructor, but this also hinders the convergence rate. The function :cpp:func:`tamaas::EPICSolver::acceleratedSolve` does not require the tweaking of a relaxation parameter, so it can be faster if the latter does not have an optimal value. However, it is not guaranteed to converge. Finally, during the first iterations, the fixed point error will be large compared to the error of the non-linear solver. It can improve performance to have the tolerance of the non-linear solver (which is the most expensive step of the fixed point solve) decrease over the iterations. This can be achieved with the decorator class :py:class:`ToleranceManager `:: from tamaas.nonlinear_solvers import ToleranceManager, DFSANESolver @ToleranceManager(start=1e-2, end=1e-9, rate=1 / 3) class Solver(DFSANESolver): pass # ... epsolver = Solver(residual) # or epsolver = ToleranceManager(1e-2, 1e-9, 1/3)(DFSANESolver)(residual) Custom Contact Solvers ---------------------- The :cpp:class:`tamaas::ContactSolver` class can be derived in Python so that users can interface with Scipy's :py:func:`scipy.optimize.minimize` function or `PETSc `_'s solvers accessible in the Python `interface `_. See ``examples/scipy_penalty.py`` for an example on how to proceed. diff --git a/examples/plasticity.py b/examples/plasticity.py index 444d787..b1151bf 100644 --- a/examples/plasticity.py +++ b/examples/plasticity.py @@ -1,85 +1,90 @@ #!/usr/bin/env python3 # # Copyright (©) 2016-2023 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 . import numpy as np import tamaas as tm from tamaas.dumpers import H5Dumper as Dumper from tamaas.nonlinear_solvers import DFSANECXXSolver as Solver, ToleranceManager from tamaas.utils import publications - # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 51, 51] flat_domain = [1, 1] system_size = [0.5] + flat_domain # Creation of model -model = tm.ModelFactory.createModel(model_type, - system_size, - discretization) +model = tm.ModelFactory.createModel(model_type, system_size, discretization) model.E = 1. model.nu = 0.3 # Setup for plasticity -residual = tm.ModelFactory.createResidual(model, - sigma_y=0.1 * model.E, - hardening=0.01 * model.E) +material = tm.materials.IsotropicHardening(model, + sigma_y=0.1 * model.E, + hardening=0.01 * model.E) +residual = tm.Residual(model, material) + # Possibly change integration method # model.setIntegrationMethod(tm.integration_method.cutoff, 1e-12) # Setup non-linear solver with variable tolerance -epsolver = ToleranceManager(1e-5, 1e-9, 1/4)(Solver)(residual) +epsolver = ToleranceManager(1e-5, 1e-9, 1 / 4)(Solver)(residual) # Setup for contact -x = np.linspace(0, system_size[1], discretization[1], - endpoint=False, dtype=tm.dtype) -y = np.linspace(0, system_size[2], discretization[2], - endpoint=False, dtype=tm.dtype) +x = np.linspace(0, + system_size[1], + discretization[1], + endpoint=False, + dtype=tm.dtype) +y = np.linspace(0, + system_size[2], + discretization[2], + endpoint=False, + dtype=tm.dtype) xx, yy = np.meshgrid(x, y, indexing='ij') R = 0.2 surface = -((xx - flat_domain[0] / 2)**2 + (yy - flat_domain[1] / 2)**2) / (2 * R) # Scatter surface across MPI processes local_surface = tm.mpi.scatter(surface) csolver = tm.PolonskyKeerRey(model, local_surface, 1e-12, tm.PolonskyKeerRey.pressure, tm.PolonskyKeerRey.pressure) # EPIC setup epic = tm.EPICSolver(csolver, epsolver, 1e-7) # Dumper dumper_helper = Dumper('hertz', 'displacement', 'stress', 'plastic_strain') model.addDumper(dumper_helper) loads = np.linspace(0.001, 0.005, 3) for i, load in enumerate(loads): epic.acceleratedSolve(load) model.dump() tm.Logger().get(tm.LogLevel.info) \ << "---> Solved load step {}/{}".format(i+1, len(loads)) # Print list of relevant publications publications() diff --git a/examples/stresses.py b/examples/stresses.py index 1114ee2..f9bf888 100644 --- a/examples/stresses.py +++ b/examples/stresses.py @@ -1,123 +1,123 @@ #!/usr/bin/env python3 # # Copyright (©) 2016-2023 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 . import argparse import os import numpy as np import tamaas as tm from tamaas.dumpers import H5Dumper as Dumper from tamaas.dumpers._helper import hdf5toVTK from tamaas.utils import publications parser = argparse.ArgumentParser( description="Hertzian tractios applied on elastic half-space") parser.add_argument("radius", type=float, help="Radius of sphere") parser.add_argument("load", type=float, help="Applied normal force") parser.add_argument("name", help="Output file name") parser.add_argument("--plots", help='Show surface normal stress', action="store_true") args = parser.parse_args() # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [32, 127, 127] system_size = [0.25, 1., 1.] # Material contants E = 1. # Young's modulus nu = 0.3 # Poisson's ratio E_star = E/(1-nu**2) # Hertz modulus # Creation of model model = tm.ModelFactory.createModel(model_type, system_size, discretization) model.E = E model.nu = nu # Setup for integral operators -residual = tm.ModelFactory.createResidual(model, 0, 0) +tm.ModelFactory.registerVolumeOperators(model) # Coordinates x = np.linspace(0, system_size[1], discretization[1], endpoint=False) y = np.linspace(0, system_size[2], discretization[2], endpoint=False) x, y = np.meshgrid(x, y, indexing='ij') center = [0.5, 0.5] r = np.sqrt((x-center[0])**2 + (y-center[1])**2) # Span of local data local_size = model.boundary_shape local_offset = tm.mpi.local_offset(r.shape) local_slice = np.s_[local_offset:local_offset+local_size[0], :] # Sphere R = args.radius P = args.load # Contact area a = (3*P*R/(4*E_star))**(1/3) p_0 = 3 * P / (2 * np.pi * a**2) # Pressure definition traction = model.traction hertz_traction = np.zeros(discretization[1:]) hertz_traction[r < a] = p_0 * np.sqrt(1 - (r[r < a] / a)**2) traction[..., 2] = hertz_traction[local_slice] # Array references displacement = model.displacement stress = model['stress'] gradient = np.zeros_like(stress) # Getting integral operators boussinesq = model.operators["boussinesq"] boussinesq_gradient = model.operators["boussinesq_gradient"] # Applying operators boussinesq(traction, displacement) boussinesq_gradient(traction, gradient) model.operators["hooke"](gradient, stress) # Dumper dumper_helper = Dumper(args.name, 'stress') model.addDumper(dumper_helper) model.dump() # Converting HDF dump to VTK with tm.mpi.sequential(): if tm.mpi.rank() == 0: hdf5toVTK(os.path.join('hdf5', f'{args.name}_0000.h5'), args.name) if args.plots: import matplotlib.pyplot as plt # noqa fig, ax = plt.subplots(1, 2) fig.suptitle("Rank {}".format(tm.mpi.rank())) ax[0].set_title("Traction") ax[1].set_title("Normal Stress") ax[0].imshow(traction[..., 2]) ax[1].imshow(stress[0, ..., 2]) fig.tight_layout() plt.show() publications() diff --git a/python/wrap/compute.cpp b/python/wrap/compute.cpp index 8900c2e..21f7f1d 100644 --- a/python/wrap/compute.cpp +++ b/python/wrap/compute.cpp @@ -1,91 +1,92 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 "computes.hh" #include "wrap.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { namespace wrap { using namespace py::literals; void wrapCompute(py::module& mod) { auto compute_mod = mod.def_submodule("compute"); + compute_mod.doc() = "Module defining basic computations on fields."; compute_mod.def( "eigenvalues", [](model_type type, Grid& eigs, const Grid& field) { eigenvalues(type, eigs, field); }, "model_type"_a, "eigenvalues_out"_a, "field"_a, "Compute eigenvalues of a tensor field"); compute_mod.def( "von_mises", [](model_type type, Grid& vm, const Grid& field) { vonMises(type, vm, field); }, "model_type"_a, "von_mises"_a, "field"_a, "Compute the Von Mises invariant of a tensor field"); compute_mod.def( "deviatoric", [](model_type type, Grid& dev, const Grid& field) { deviatoric(type, dev, field); }, "model_type"_a, "deviatoric"_a, "field"_a, "Compute the deviatoric part of a tensor field"); compute_mod.def( "to_voigt", [](const Grid& field) { if (field.getNbComponents() != 9) TAMAAS_EXCEPTION("Wrong number of components to symmetrize"); Grid voigt(field.sizes(), 6); Loop::loop( [] CUDA_LAMBDA(MatrixProxy in, SymMatrixProxy out) { out.symmetrize(in); }, range>(field), range>(voigt)); return voigt; }, "Convert a 3D tensor field to voigt notation", py::return_value_policy::copy); compute_mod.def( "from_voigt", [](const Grid& field) { if (field.getNbComponents() != 6) TAMAAS_EXCEPTION("Wrong number of components to symmetrize"); Grid full(field.sizes(), 9); Loop::loop([] CUDA_LAMBDA( SymMatrixProxy in, MatrixProxy out) { out.fromSymmetric(in); }, range>(field), range>(full)); return full; }, "Convert a 3D tensor field to voigt notation", py::return_value_policy::copy); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/materials.cpp b/python/wrap/materials.cpp index d7bb2c6..fdbda13 100644 --- a/python/wrap/materials.cpp +++ b/python/wrap/materials.cpp @@ -1,81 +1,94 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 "materials/isotropic_hardening.hh" #include "materials/material.hh" #include "wrap.hh" #include #include namespace tamaas { namespace wrap { using namespace py::literals; class PyMaterial : public Material { using Field = typename Material::Field; public: using Material::Material; void computeStress(Field& stress, const Field& strain, const Field& strain_increment) override { // NOLINTNEXTLINE(readability-else-after-return) PYBIND11_OVERLOAD_PURE(void, Material, computeStress, stress, strain, strain_increment); } void computeEigenStress(Field& stress, const Field& strain, const Field& strain_increment) override { // NOLINTNEXTLINE(readability-else-after-return) PYBIND11_OVERLOAD_PURE(void, Material, computeEigenStress, stress, strain, strain_increment); } void update() override { // NOLINTNEXTLINE(readability-else-after-return) PYBIND11_OVERRIDE_PURE(void, Material, update); } }; void wrapMaterialInterface(py::module& mod) { py::class_>(mod, "Material") .def(py::init(), py::keep_alive<1, 2>()) - .def("computeStress", &Material::computeStress) - .def("computeEigenStress", &Material::computeEigenStress); + .def("computeStress", &Material::computeStress, "stress"_a, "strain"_a, + "strain_increment"_a, + "Compute the total stress from the total strain and total strain " + "increment") + .def("computeEigenStress", &Material::computeEigenStress, "stress"_a, + "strain"_a, "strain_increment"_a, + "Compute the stress from the eigen strain (e.g. plastic strain " + "increment)"); } void wrapIsotropicHardening(py::module& mod) { py::class_>( mod, "IsotropicHardening") .def(py::init(), "model"_a, "sigma_y"_a, - "hardening"_a, py::keep_alive<1, 2>()) + "hardening"_a, py::keep_alive<1, 2>(), + R"-(Create an isotropic linear hardening material. +:param model: the model on which to define the material +:param sigma_y: the (von Mises) yield stress +:param hardening: the hardening modulus)-") .def("computeInelasticDeformationIncrement", - &IsotropicHardening::computeInelasticDeformationIncrement); + &IsotropicHardening::computeInelasticDeformationIncrement, + "stress"_a, "strain"_a, "strain_increment"_a, + "Compute the plastic strain increment"); } void wrapMaterials(py::module& m) { auto mod = m.def_submodule("materials"); + mod.doc() = "Module defining classes that implement constitutive relations"; wrapMaterialInterface(mod); wrapIsotropicHardening(mod); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp index 05a52f5..206d5eb 100644 --- a/python/wrap/model.cpp +++ b/python/wrap/model.cpp @@ -1,524 +1,539 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 "elastic_functional.hh" #include "field_container.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 #include #include #include /* -------------------------------------------------------------------------- */ 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); } }; 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"); } } void wrapFieldContainer(py::module& mod) { py::class_(mod, "FieldContainer") .def(py::init<>()) .def("__getitem__", &FieldContainer::at, py::return_value_policy::reference_internal) .def( "__setitem__", [](FieldContainer& f, const std::string& key, numpy grid) { f[key] = instanciateFromNumpy(grid); }, py::keep_alive<1, 3>()) .def( "__setitem__", [](FieldContainer& f, const std::string& key, numpy grid) { f[key] = instanciateFromNumpy(grid); }, py::keep_alive<1, 3>()) .def( "__setitem__", [](FieldContainer& f, const std::string& key, numpy grid) { f[key] = instanciateFromNumpy(grid); }, py::keep_alive<1, 3>()) .def( "__contains__", [](const FieldContainer& m, const std::string& key) { const auto fields = m.fields(); return std::find(fields.begin(), fields.end(), key) != fields.end(); }, "Test field existence") .def( "__iter__", [](const FieldContainer& m) { const auto& fields = m.fields_map(); return py::make_key_iterator(fields.cbegin(), fields.cend()); }, py::keep_alive<0, 1>(), "Iterator on fields") /// Deprecated interface .def( "registerField", [](Model& m, const std::string& name, numpy field) { TAMAAS_DEPRECATE("registerField()", "the [] operator"); m[name] = instanciateFromNumpy(field); }, "field_name"_a, "field"_a, py::keep_alive<1, 3>()) .def( "getField", [](const Model& m, const std::string& name) -> decltype(m.at(name)) { TAMAAS_DEPRECATE("getField()", "the [] operator"); return m.at(name); }, "field_name"_a, py::return_value_policy::reference_internal) .def( "getFields", [](const Model& m) { TAMAAS_DEPRECATE("getFields()", "list(model)"); return m.fields(); }, "Return fields list"); } /// 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, "Compute functional value") .def("computeGradF", &functional::Functional::computeGradF, "Compute functional gradient"); py::class_( mod, "ElasticFunctionalPressure", func) .def(py::init&>()); py::class_(mod, "ElasticFunctionalGap", func) .def(py::init&>()); py::class_ adh(mod, "AdhesionFunctional", func); adh.def_property("parameters", &functional::AdhesionFunctional::getParameters, &functional::AdhesionFunctional::setParameters, "Parameters dictionary") .def("setParameters", [](functional::AdhesionFunctional& f, const std::map& m) { TAMAAS_DEPRECATE("setParameters()", "the parameters property"); f.setParameters(m); }); py::class_( mod, "ExponentialAdhesionFunctional", adh, "Potential of the form F = -γ·exp(-g/ρ)") .def(py::init&>(), "surface"_a); py::class_( mod, "MaugisAdhesionFunctional", adh, "Cohesive zone potential F = H(g - ρ)·γ/ρ") .def(py::init&>(), "surface"_a); py::class_( mod, "SquaredExponentialAdhesionFunctional", adh, "Potential of the form F = -γ·exp(-0.5·(g/ρ)²)") .def(py::init&>(), "surface"_a); } /// Wrap IntegralOperator void wrapIntegralOperator(py::module& mod) { py::class_(mod, "IntegralOperator") .def("apply", [](IntegralOperator& op, numpy input, numpy 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 input, numpy 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 X) -> GridBase { GridBaseNumpy x(X); auto y = op.matvec(x); return y; }, py::return_value_policy::move); py::enum_(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_(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 void wrapModelTypeTrait(py::module& mod) { using trait = model_type_traits; py::class_(mod, trait::repr) .def_property_readonly_static( "dimension", [](const py::object&) { return trait::dimension; }, "Dimension of computational domain") .def_property_readonly_static( "components", [](const py::object&) { return trait::components; }, "Number of components of vector fields") .def_property_readonly_static( "boundary_dimension", [](const py::object&) { return trait::boundary_dimension; }, "Dimension of boundary of computational domain") .def_property_readonly_static( "voigt", [](const py::object&) { return trait::voigt; }, "Number of components of symmetrical tensor fields") .def_property_readonly_static( "indices", [](const py::object&) { return trait::indices; }); } /// Wrap Models void wrapModelClass(py::module& mod) { py::enum_(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"); + trait_mod.doc() = "Defines type trait objects."; 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(py::init(py::overload_cast&, const std::vector&>( &ModelFactory::createModel)), 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_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("setIntegrationMethod", &Model::setIntegrationMethod) .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( "applyElasticity", [](Model& model, numpy stress, numpy 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( "__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& (Model::*)() const) & Model::getTraction, "Surface traction field") .def_property_readonly("displacement", (const GridBase& (Model::*)() const) & Model::getDisplacement, "Displacement field"); py::class_>( 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_(mod, "ModelFactory") .def_static( "createModel", py::overload_cast&, const std::vector&>( &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(&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. + .def_static( + "createResidual", + [](Model& model, Real sigma_y, Real hardening) { + TAMAAS_DEPRECATE("createResidual", "Residual constructor"); + return ModelFactory::createResidual(model, sigma_y, hardening); + }, + "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.") .def_static("setIntegrationMethod", &ModelFactory::setIntegrationMethod, "operator"_a, "method"_a, "cutoff"_a, "Set the integration method (linear or cutoff) for a volume " "integral operator"); } /// Wrap residual class void wrapResidual(py::module& mod) { // TODO adapt to n-dim py::class_(mod, "Residual") .def(py::init>(), - py::keep_alive<1, 2>()) + py::keep_alive<1, 2>(), "model"_a, "material"_a, + R"-(Create a residual object with a material. +Defines the following residual equation: + ɛ - ∇N[σ(ɛ)] - ∇M[t] = 0 + +Where σ(ɛ) is the *eigenstress* associated with the constitutive law. + +:param model: the model on which to define the residual +:param material: material object which defines a constitutive law +)-") .def("computeResidual", &Residual::computeResidual) .def("updateState", &Residual::updateState) .def("computeResidualDisplacement", &Residual::computeResidualDisplacement) .def(TAMAAS_DEPRECATE_ACCESSOR(getVector, Residual, "vector"), py::return_value_policy::reference_internal) .def( "getStress", [](const Residual& res) { TAMAAS_DEPRECATE("getStress", "model[\"stress\"]"); const auto& stress = res.getModel().field("stress"); Grid view(res.getModel().getDiscretization(), stress.getNbComponents(), stress.view()); return view; }, py::return_value_policy::reference_internal) .def("setIntegrationMethod", [](Residual& residual, integration_method method, Real cutoff) { TAMAAS_DEPRECATE("setIntegrationMethod", "model.setIntegrationMethod()"); residual.getModel().setIntegrationMethod(method, cutoff); }) .def_property_readonly("vector", &Residual::getVector) .def_property_readonly( "model", static_cast(&Residual::getModel)); } void wrapModel(py::module& mod) { wrapFieldContainer(mod); wrapBEEngine(mod); wrapModelClass(mod); wrapModelFactory(mod); wrapFunctionals(mod); wrapResidual(mod); wrapIntegralOperator(mod); } } // namespace wrap } // namespace tamaas diff --git a/python/wrap/mpi.cpp b/python/wrap/mpi.cpp index 18d8cdd..f45caa1 100644 --- a/python/wrap/mpi.cpp +++ b/python/wrap/mpi.cpp @@ -1,101 +1,102 @@ /* * SPDX-License-Indentifier: AGPL-3.0-or-later * * Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * Copyright (©) 2020-2023 Lucas Frérot * * 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 "mpi_interface.hh" #include "partitioner.hh" #include "wrap.hh" #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ namespace wrap { using namespace py::literals; /* -------------------------------------------------------------------------- */ void wrapMPI(py::module& mod) { auto mpi_mod = mod.def_submodule("mpi"); + mpi_mod.doc() = "Module defining convenience MPI routines."; py::class_(mpi_mod, "sequential") .def(py::init<>()) .def("__enter__", [](mpi::sequential& /*obj*/) { mpi::sequential::enter(); }) .def("__exit__", [](mpi::sequential& /*obj*/, py::object /*unused*/, py::object /*unused*/, py::object /*unused*/) { mpi::sequential::exit(); }); mpi_mod.def( "local_shape", [](std::vector global) { switch (global.size()) { case 1: return Partitioner<1>::local_size(global); case 2: return Partitioner<2>::local_size(global); default: TAMAAS_EXCEPTION("Please provide a 1D/2D shape"); } }, "global_shape"_a, "Gives the local size of a 1D/2D global shape"); mpi_mod.def( "global_shape", [](std::vector local) { switch (local.size()) { case 1: return Partitioner<1>::global_size(local); case 2: return Partitioner<2>::global_size(local); default: TAMAAS_EXCEPTION("Please provide a 1D/2D shape"); } }, "local_shape"_a, "Gives the global shape of a 1D/2D local shape"); mpi_mod.def( "local_offset", [](std::vector global) { switch (global.size()) { case 1: return Partitioner<1>::local_offset(global); case 2: return Partitioner<2>::local_offset(global); default: TAMAAS_EXCEPTION("Please provide a 1D/2D shape"); } }, "global_shape"_a, "Gives the local offset of a 1D/2D global shape"); mpi_mod.def("gather", &Partitioner<2>::gather, py::return_value_policy::move, "Gather 2D surfaces"); mpi_mod.def("scatter", &Partitioner<2>::scatter, py::return_value_policy::move, "Scatter 2D surfaces"); mpi_mod.def( "size", []() { return mpi::size(); }, "Returns the number of MPI processes"); mpi_mod.def( "rank", []() { return mpi::rank(); }, "Returns the rank of the local process"); } } // namespace wrap /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/tests/conftest.py b/tests/conftest.py index 62a92c6..f58f27b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,309 +1,310 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . from __future__ import division, print_function import subprocess import sys from functools import reduce from operator import mul as multiply import numpy as np import pytest from numpy.linalg import norm try: import tamaas as tm except ModuleNotFoundError as e: print(e, file=sys.stderr) print("Use 'scons test' to run tests", file=sys.stderr) sys.exit(1) type_list = [ tm.model_type.basic_1d, tm.model_type.basic_2d, tm.model_type.surface_1d, tm.model_type.surface_2d, tm.model_type.volume_1d, tm.model_type.volume_2d ] def get_libtamaas_file(): try: ldd_out = bytes(subprocess.check_output(['ldd', tm._tamaas.__file__])) for line in filter(lambda x: 'Tamaas' in x, ldd_out.decode().split('\n')): path = line.split(' => ')[1] return path.split(' ')[0] return "static link" except subprocess.CalledProcessError: return "N/A" def pytest_report_header(config): print('tamaas build type: {}\nmodule file: {}\nwrapper: {}\nlibTamaas: {}'. format(tm.TamaasInfo.build_type, tm.__file__, tm._tamaas.__file__, get_libtamaas_file())) def pytest_addoption(parser): parser.addoption('--log-debug', action='store_true', default=False) def mtidfn(val): return str(val).split('.')[-1] def pkridfn(val): return { tm.PolonskyKeerRey.pressure: "pkr_pressure", tm.PolonskyKeerRey.gap: "pkr_gap" }[val] def profile(func, domain, N, modes, amplitude): coords = [ np.linspace(0, d, n, endpoint=False, dtype=tm.dtype) for n, d in zip(N, domain) ] coords = np.meshgrid(*coords, indexing='ij') sines = [func(2 * np.pi * x * m) for x, m in zip(coords, modes)] return reduce(multiply, [amplitude] + sines) @pytest.fixture(scope='session', autouse=True) def log_level(pytestconfig): if pytestconfig.option.log_debug: tm.set_log_level(tm.LogLevel.debug) class HertzFixture: def __init__(self, n, load): self.domain_size = 1 self.n = n self.load = load self.curvature = 0.1 self.radius = 1. / self.curvature self.e_star = 1. self.a = (3 * load / (4 * self.curvature * self.e_star))**(1. / 3.) self.x = np.linspace(-self.domain_size / 2., self.domain_size / 2., self.n, dtype=tm.dtype) self.y = self.x.copy() self.x, self.y = np.meshgrid(self.x, self.y) self._computeSurface() self._computePressure() self._computeDisplacement() def _computeDisplacement(self): r = np.sqrt(self.x**2 + self.y**2) self.displacement = np.zeros_like(r) contact = r < self.a self.displacement[contact] = self.surface[contact] self.displacement[~contact] = \ (self.surface[~contact] + self.a / (np.pi * self.radius) * np.sqrt(r[~contact]**2 - self.a**2) + (r[~contact]**2 - 2 * self.a**2) / (np.pi * self.radius) * np.arccos(self.a / r[~contact])) def _computePressure(self): r = np.sqrt(self.x**2 + self.y**2) self.pressure = np.zeros_like(r) contact = np.where(r < self.a) self.pressure[contact] = \ 2 * self.e_star / (np.pi * self.radius) \ * np.sqrt(self.a**2 - r[contact]**2) def _computeSurface(self): self.surface = -1. / (2 * self.radius) * (self.x**2 + self.y**2) @pytest.fixture(scope="package") def hertz(): return HertzFixture(1024, 0.00001) @pytest.fixture(scope="package") def hertz_coarse(): return HertzFixture(512, 0.0001) @pytest.fixture(scope="package", params=[tm.PolonskyKeerRey.pressure], ids=pkridfn) def pkr(hertz, request): model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [hertz.domain_size, hertz.domain_size], [hertz.n, hertz.n]) model.E, model.nu = hertz.e_star, 0 solver = tm.PolonskyKeerRey(model, hertz.surface, 1e-12, request.param, request.param) solver.solve(hertz.load) return model, hertz class WestergaardFixture: def __init__(self, n, load): self.domain_size = 1. self.lamda = 1. self.delta = 0.1 self.e_star = 1. self.n = n self.p_star = np.pi * self.e_star * self.delta / self.lamda self.load = load * self.p_star self.a = self.lamda / np.pi \ * np.arcsin(np.sqrt(self.load / self.p_star)) self.x = np.linspace(-self.domain_size / 2., self.domain_size / 2., self.n, endpoint=False, dtype=tm.dtype) self._computeSurface() self._computePressure() self._computeDisplacement() def _computeSurface(self): self.surface = self.delta * np.cos(2 * np.pi * self.x / self.lamda) def _computePressure(self): self.pressure = np.zeros_like(self.surface) contact = np.where(np.abs(self.x) < self.a) self.pressure[contact] = 2 * self.load \ * (np.cos(np.pi * self.x[contact] / self.lamda) / np.sin(np.pi * self.a / self.lamda)**2) \ * np.sqrt(np.sin(np.pi * self.a / self.lamda)**2 - np.sin(np.pi * self.x[contact] / self.lamda)**2) def _computeDisplacement(self): psi = np.pi * np.abs(self.x) / self.lamda psi_a = np.pi * self.a / self.lamda with np.errstate(invalid='ignore'): # get some warnings out of the way self.displacement = ( np.cos(2 * psi) + 2 * np.sin(psi) * np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2) - 2 * np.sin(psi_a)**2 * np.log( (np.sin(psi) + np.sqrt(np.sin(psi)**2 - np.sin(psi_a)**2)) / np.sin(psi_a))) contact = np.where(np.abs(self.x) < self.a) self.displacement[contact] = np.cos(2 * psi[contact]) self.displacement *= self.load * self.lamda / (np.pi * self.e_star * np.sin(psi_a)**2) @pytest.fixture(scope="package") def westergaard(): return WestergaardFixture(19683, 0.1) class PatchWestergaard: def __init__(self, model_type): dim = tm.type_traits[model_type].dimension bdim = tm.type_traits[model_type].boundary_dimension domain = [1.] * dim size = [6] * dim self.modes = np.random.randint(1, 3, (bdim, )) self.model = tm.ModelFactory.createModel(model_type, domain, size) self.model.E = 3. self.model.nu = 0. self.pressure = profile(np.cos, self.model.boundary_system_size, self.model.boundary_shape, self.modes, 1) self.solution = profile( np.cos, self.model.boundary_system_size, self.model.boundary_shape, self.modes, 1 / (np.pi * self.model.E_star * norm(self.modes))) @pytest.fixture(scope="package", params=set(type_list) - {tm.model_type.volume_1d}, ids=mtidfn) def patch_westergaard(request): return PatchWestergaard(request.param) class UniformPlasticity: def __init__(self, model_type, domain, sizes): self.n = sizes self.domain = domain self.model = tm.ModelFactory.createModel(model_type, domain, sizes) self.E_h = 0.1 self.sigma_y = 0.01 - self.residual = tm.ModelFactory.createResidual(self.model, - sigma_y=self.sigma_y, - hardening=self.E_h) + self.mat = tm.materials.IsotropicHardening(self.model, + sigma_y=self.sigma_y, + hardening=self.E_h) + self.residual = tm.Residual(self.model, self.mat) self.model.E = 1. self.model.nu = 0. def solution(self, p): E, nu = self.model.E, self.model.nu E_h, sigma_y = self.E_h, self.sigma_y mu = E / (2 * (1 + nu)) strain = -1 / (mu + E_h) * (p * (3 * mu + E_h) / (2 * mu) - np.sign(p) * sigma_y) dep = (2 * mu * np.abs(strain) - sigma_y) / (3 * mu + E_h) plastic_strain = np.sign(strain) / 2 * dep * np.array([ -1, -1, 2, 0, 0, 0, ], dtype=tm.dtype) stress = 2 * mu * (np.array([0, 0, strain, 0, 0, 0], dtype=tm.dtype) - plastic_strain) return { "stress": stress, "plastic_strain": plastic_strain, "cumulated_plastic_strain": dep }, { "stress": mu, "plastic_strain": 1, "cumulated_plastic_strain": 1 } # Tests plane-strain situations with Ny = 1 @pytest.fixture(params=[1, 4], ids=mtidfn) def patch_isotropic_plasticity(request): return UniformPlasticity(tm.model_type.volume_2d, [1., 1., 1.], [4, 4, request.param]) diff --git a/tests/test_boussinesq_surface.py b/tests/test_boussinesq_surface.py index d29a231..f7165e3 100644 --- a/tests/test_boussinesq_surface.py +++ b/tests/test_boussinesq_surface.py @@ -1,79 +1,79 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . import numpy as np import tamaas as tm from numpy.linalg import norm def test_boussinesq_surface(): # Definition of modeled domain model_type = tm.model_type.volume_2d discretization = [2, 128, 128] system_size = [1., 1., 1.] - x, y = [np.linspace(0, size, n, endpoint=False, dtype=tm.dtype) - for size, n in zip(system_size[1:], discretization[1:])] + x, y = [ + np.linspace(0, size, n, endpoint=False, dtype=tm.dtype) + for size, n in zip(system_size[1:], discretization[1:]) + ] xx, yy = np.meshgrid(x, y, sparse=True) tractions = np.zeros(discretization[1:] + [3], dtype=tm.dtype) for i in range(3): - omega = np.random.randint(1, 20, (2,)) * 2 * np.pi - tractions[:, :, i] = np.cos(omega[0]*xx)*np.cos(omega[1]*yy) + omega = np.random.randint(1, 20, (2, )) * 2 * np.pi + tractions[:, :, i] = np.cos(omega[0] * xx) * np.cos(omega[1] * yy) # Material contants - E = 1. # Young's modulus - nu = 0.3 # Poisson's ratio + E = 1. # Young's modulus + nu = 0.3 # Poisson's ratio # Creation of model - model_vol = tm.ModelFactory.createModel(model_type, - system_size, + model_vol = tm.ModelFactory.createModel(model_type, system_size, discretization) - model_surf = tm.ModelFactory.createModel(tm.model_type.surface_2d, - [1, 1], + model_surf = tm.ModelFactory.createModel(tm.model_type.surface_2d, [1, 1], tractions.shape[:-1]) model_vol.E = model_surf.E = E model_vol.nu = model_surf.nu = nu # Setup for integral operators - tm.ModelFactory.createResidual(model_vol, 0, 0) + tm.ModelFactory.registerVolumeOperators(model_vol) # Pressure definition model_vol['traction'][...] = tractions model_surf['traction'][...] = tractions model_surf.solveNeumann() # Applying operator boussinesq = model_vol.operators["boussinesq"] boussinesq(model_vol['traction'], model_vol['displacement']) # # Dumper # dumper_helper = UVWDumper(residual, args.name) # model.addDumper(dumper_helper) # model.dump() # print("Done") for i in range(3): error = norm(model_vol.displacement[0, :, :, i] - model_surf.displacement[:, :, i]) \ / norm(tractions[:, :, i]) assert error < 1e-16 diff --git a/tests/test_epic.py b/tests/test_epic.py index ff4a199..cb3f87f 100644 --- a/tests/test_epic.py +++ b/tests/test_epic.py @@ -1,80 +1,78 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . from __future__ import division, print_function import pytest from conftest import pkridfn import tamaas as tm from numpy.linalg import norm from tamaas.nonlinear_solvers import DFSANESolver @pytest.fixture(scope="module", params=[tm.PolonskyKeerRey.pressure], ids=pkridfn) def model_setup(hertz_coarse, request): models, solvers = [], [] for mtype in [tm.model_type.basic_2d, tm.model_type.volume_2d]: dim = tm.type_traits[mtype].dimension n = [hertz_coarse.n] * dim d = [hertz_coarse.domain_size] * dim if dim == 3: n[0] = 2 d[0] = 0.001 model = tm.ModelFactory.createModel(mtype, d, n) model.E, model.nu = hertz_coarse.e_star, 0 solver = tm.PolonskyKeerRey(model, hertz_coarse.surface, 1e-12, - request.param, - request.param) + request.param, request.param) models.append(model) solvers.append(solver) - epsolver = DFSANESolver(tm.ModelFactory.createResidual( - model, - sigma_y=1e2 * model.E, - hardening=0 - )) + mat = tm.materials.IsotropicHardening(model, + sigma_y=1e2 * model.E, + hardening=0) + epsolver = DFSANESolver(tm.Residual(model, mat)) solvers[1] = tm.EPICSolver(solvers[1], epsolver, 1e-12) for solver in solvers: solver.solve(hertz_coarse.load) return models, hertz_coarse @pytest.mark.xfail(tm.TamaasInfo.backend == 'tbb', reason='TBB reductions are unstable?') def test_energy(model_setup): """Test the full elastic-plastic solve step in elasticity""" # We use computed values from basic PKR to test (model_elastic, model_plastic), hertz = model_setup error = norm( - (model_plastic.traction[..., 2] - model_elastic.traction) - * (model_plastic.displacement[0, ..., 2] - model_elastic.displacement)) + (model_plastic.traction[..., 2] - model_elastic.traction) * + (model_plastic.displacement[0, ..., 2] - model_elastic.displacement)) error /= norm(hertz.pressure * hertz.displacement) assert error < 1e-16 diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index 10cbf54..4241b8f 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,230 +1,225 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . import tamaas as tm import numpy as np import pytest from numpy.linalg import norm from register_integral_operators import register_kelvin_force, \ set_integration_method @pytest.fixture(params=[tm.integration_method.linear, tm.integration_method.cutoff], ids=['linear', 'cutoff']) def integral_fixture(request): return request.param def test_kelvin_volume_force(integral_fixture): N = 65 E = 1. nu = 0.3 mu = E / (2*(1+nu)) domain = np.array([1.] * 3) omega = 2 * np.pi * np.array([1, 1]) / domain[:2] omega_ = norm(omega) discretization = [N] * 3 model = tm.ModelFactory.createModel(tm.model_type.volume_2d, domain, discretization) model.E = E model.nu = nu register_kelvin_force(model) kelvin = model.operators['kelvin_force'] set_integration_method(kelvin, integral_fixture, 1e-12) coords = [np.linspace(0, domain[i], discretization[i], endpoint=False, dtype=tm.dtype) for i in range(2)]\ + [np.linspace(0, domain[2], discretization[2])] x, y = np.meshgrid(*coords[:2], indexing='ij') displacement = model['displacement'] source = np.zeros_like(displacement) # The integral of forces should stay constant source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y) * (N-1) kelvin(source, displacement) z = coords[2] - 0.5 z, x, y = np.meshgrid(z, *coords[:2], indexing='ij') solution = np.zeros_like(source) solution[:, :, :, 0] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * omega[0]*z*np.cos(omega[0]*x)*np.sin(omega[1]*y) solution[:, :, :, 1] = -np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * omega[1]*z*np.sin(omega[0]*x)*np.cos(omega[1]*y) solution[:, :, :, 2] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) \ * (3-4*nu + omega_*np.abs(z))*np.sin(omega[0]*x)*np.sin(omega[1]*y) error = norm(displacement - solution) / norm(solution) assert error < 5e-2 def test_mindlin(integral_fixture): # Definition of modeled domain # tm.set_log_level(tm.LogLevel.debug) model_type = tm.model_type.volume_2d discretization = [126, 128, 128] system_size = [1., 3., 3.] integration_method = integral_fixture print(integration_method) # Material contants E = 1. # Young's modulus nu = 0.3 # Poisson's ratio - h = 0.01 * E # Hardening modulus - sigma_y = 0.3 * E # Initial yield stress - # Creation of model - model = tm.ModelFactory.createModel(model_type, - system_size, - discretization) + model = tm.Model(model_type, system_size, discretization) model.E = E model.nu = nu mu = E / (2*(1+nu)) lamda = E * nu / ((1+nu) * (1-2*nu)) # Setup for integral operators - residual = tm.ModelFactory.createResidual(model, sigma_y, h) + tm.ModelFactory.registerVolumeOperators(model) model.setIntegrationMethod(integration_method, 1e-12) # Coordinates x = np.linspace(0, system_size[1], discretization[1], endpoint=False, dtype=tm.dtype) y = np.linspace(0, system_size[2], discretization[2], endpoint=False, dtype=tm.dtype) z = np.linspace(0, system_size[0], discretization[0], endpoint=True, dtype=tm.dtype) z, x, y = np.meshgrid(z, x, y, indexing='ij') # Inclusion definition a, c = 0.1, 0.2 center = [system_size[1] / 2, system_size[2] / 2, c] r = np.sqrt((x-center[0])**2 + (y-center[1])**2 + (z-center[2])**2) ball = r < a # Eigenstrain definition alpha = 1 # constant isotropic strain beta = (3 * lamda + 2 * mu) * alpha * np.eye(3) eigenstress = np.zeros(discretization + [3, 3], dtype=tm.dtype) eigenstress[ball, ...] = beta eigenstress = tm.compute.to_voigt(eigenstress.reshape(discretization + [9])) # Array references - gradient = residual.vector - stress = model["stress"] + stress = np.zeros(discretization + [6], dtype=tm.dtype) + gradient = np.zeros_like(stress) # Applying operator # mindlin = model.operators["mindlin"] mindlin_gradient = model.operators["mindlin_gradient"] # Not testing displacements yet # mindlin(eigenstress, model['displacement']) # Applying gradient mindlin_gradient(eigenstress, gradient) model.operators['hooke'](gradient, stress) stress -= eigenstress # Normalizing stess as in Mura (1976) T, alpha = 1., 1. beta = alpha * T * (1+nu) / (1-nu) stress *= 1. / (2 * mu * beta) n = discretization[1] // 2 vertical_stress = stress[:, n, n, :] z_all = np.linspace(0, 1, discretization[0], dtype=tm.dtype) sigma_z = np.zeros_like(z_all) sigma_t = np.zeros_like(z_all) inclusion = np.abs(z_all - c) < a # Computing stresses for exterior points z = z_all[~inclusion] R_1 = np.abs(z - c) R_2 = np.abs(z + c) sigma_z[~inclusion] = 2*mu*beta*a**3/3 * ( 1 / R_1**3 - 1 / R_2**3 - 18*z*(z+c) / R_2**5 + 3*(z+c)**2 / R_2**5 - 3*(z-c)**2 / R_1**5 + 30*z*(z+c)**3 / R_2**7 ) sigma_t[~inclusion] = 2*mu*beta*a**3/3 * ( 1 / R_1**3 + (3-8*nu) / R_2**3 - 6*z*(z+c) / R_2**5 + 12*nu*(z+c)**2 / R_2**5 ) # Computing stresses for interior points z = z_all[inclusion] R_1 = np.abs(z - c) R_2 = np.abs(z + c) sigma_z[inclusion] = 2*mu*beta*a**3/3 * ( - 2 / a**3 - 1 / R_2**3 - 18*z*(z+c) / R_2**5 + 3*(z+c)**2 / R_2**5 + 30*z*(z+c)**3 / R_2**7 ) sigma_t[inclusion] = 2*mu*beta*a**3/3 * ( - 2 / a**3 + (3-8*nu) / R_2**3 - 6*z*(z+c) / R_2**5 + 12*nu*(z+c)**2 / R_2**5 ) # This test can be used to debug if False: import matplotlib.pyplot as plt plt.plot(z_all, vertical_stress[:, 2]) plt.plot(z_all, sigma_z / (2 * mu * beta)) plt.figure() plt.plot(z_all, vertical_stress[:, 0]) plt.plot(z_all, sigma_t / (2 * mu * beta)) plt.show() z_error = norm(vertical_stress[:, 2] - sigma_z / (2 * mu * beta)) \ / discretization[0] t_error = norm(vertical_stress[:, 0] - sigma_t / (2 * mu * beta)) \ / discretization[0] assert z_error < 1e-3 and t_error < 1e-3 diff --git a/tests/test_restart.py b/tests/test_restart.py index 753e85e..cce5e58 100644 --- a/tests/test_restart.py +++ b/tests/test_restart.py @@ -1,54 +1,55 @@ # -*- coding: utf-8 -*- # # Copyright (©) 2016-20 EPFL (École Polytechnique Fédérale de Lausanne), # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # Copyright (©) 2020-2023 Lucas Frérot # # 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 . import pytest import tamaas as tm from numpy.testing import assert_allclose from tamaas.dumpers import NumpyDumper from tamaas.nonlinear_solvers import DFSANECXXSolver as CXXSolver def test_restart_plasticity(patch_isotropic_plasticity, tmp_path): model = patch_isotropic_plasticity.model residual = patch_isotropic_plasticity.residual applied_pressure = 0.1 solver = CXXSolver(residual) solver.tolerance = 1e-15 pressure = model['traction'][..., 2] pressure[:] = applied_pressure solver.solve() solver.updateState() path = tmp_path / 'restart' NumpyDumper(path, all_fields=True, mkdir=False) << model new_model = NumpyDumper.read( path.with_stem("restart_0000").with_suffix('.npz') ) - new_res = tm.ModelFactory.createResidual(model, sigma_y=1e-2, hardening=0.1) + mat = tm.materials.IsotropicHardening(model, sigma_y=1e-2, hardening=0.1) + new_res = tm.Residual(model, mat) new_solver = CXXSolver(residual) for key in new_model: print(key) assert_allclose((model[key] - new_model[key]), 0, atol=3e-15)