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. .. _units: Units in Tamaas --------------- All quantities in Tamaas are unitless. To put real units onto them, one can use the following equation, which must have a consistent set of units: .. math:: \frac{u}{L} = \frac{\pi}{k} \frac{p}{\frac{E}{1-\nu^2}}. It relates the surface displacement :math:`u` to the pressure :math:`p`, with the Young's modulus :math:`E`, the Poisson coefficient :math:`\nu` and the *horizontal* physical size of the system :math:`L` (i.e. the horizontal period of the system). The wavenumber :math:`k` is adimensional. This means that :math:`L` is the natural unit for lengths (which can be specified when creating the ``Model`` object), and :math:`E^* := E / (1 - \nu^2)` is the natural unit for pressures (which can be accessed with ``model.E_star``). All other units (areas, forces, energies, etc.) are derived from these two units. In general, it is a good idea to work with a unit set that gives numerical values of the same order of magnitude to minimize floating point errors in algebraic operations. While this is not always possible, because in elastic contact we have :math:`u \sim h_\mathrm{rms}` and :math:`p \sim h'_\mathrm{rms} E^*`, which are dependent, one should avoid using meters if expected lengths are nanometers for example. 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 ` requires the model type, the physical size of the system, and the number of points in each direction:: physical_size = [1., 1.] discretization = [512, 512] model = tm.Model(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.Model(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 constructor to :py:class:`Model ` expects the *global* system sizes and discretization of the model. .. tip:: Deep copies of model objects can be done with the :py:func:`copy.deepcopy` function:: import copy model_copy = copy.deepcopy(model) Note that it only copies fields, not operators or dumpers. Model properties ~~~~~~~~~~~~~~~~ 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 .. tip:: The Greenwood--Tripp equivalence, sometimes called Johnson equivalence, allows one to model a normal, frictionless contact between two homogeneous linear, isotropic, elastic materials with rough surfaces :math:`h_\alpha` and :math:`h_\beta` as the contact of a rigid-rough surface :math:`h_\beta - h_\alpha` and a flat elastic solid with an equivalent contact modulus :math:`E^*`: .. math:: \frac{1}{E^*} = \frac{1 - \nu_\alpha^2}{E_\alpha^2} + \frac{1 - \nu_\beta^2}{E_\beta^2} = \frac{1}{E^*_\alpha} + \frac{1}{E^*_\beta} Since a model only defines a single set of elastic properties :math:`(E, \nu)`, one can model contact between two different elastic materials by computing :math:`E^*` *a priori*, and setting :math:`(E, \nu) = (E^*, 0)`. The total displacement :math:`u` computed by, for example, a contact solver can then be redistributed in proportion according to: .. math:: \frac{1}{E^*}u = \frac{1}{E_\alpha^*}u_\alpha + \frac{1}{E_\beta^*}u_\beta Note that is equivalence breaks down with friction, plasticity, heterogeneities, anisotropy, etc. Fields ~~~~~~ Fields can be easlily accessed with the ``[]`` operator, similar to Python's dictionaries:: surface_traction = model['traction'] # surface_traction is a numpy array To know what fields are available, you can call the :py:class:`list` function on a model (``list(model)``). You can add new fields to a model object with the ``[]`` operator:``model['new_field'] = new_field``, which is convenient for dumping fields that are computed outside of Tamaas. .. note:: The fields ``traction`` and ``displacement`` are always registered in models, and are accessible via :py:attr:`model.traction ` and :py:attr:`model.displacement `. A model can also be used to compute stresses from a strain field:: import numpy strain = numpy.zeros(model.shape + [6]) # Mandel--Voigt notation stress = numpy.zeros_like(strain) model.applyElasticity(stress, strain) .. tip:: ``print(model)`` gives a lot of information about the model: the model type, shape, registered fields, and more! Manipulating fields ................... Fields are stored in C-contiguous order. The last dimension of any field is always the vector/tensor components, if the number of components is greater than one. The other dimensions are the space dimensions, in ``z, x, (y)`` order for volumetric quantities and ``x, (y)`` for boundary quantities. .. tip:: To select normal component of the surface displacement of a ``volume_2d`` model, one can write:: model.displacement[0, ..., 2] Symmetric tensors are stored in Mandel--Voigt notation as vectors in the form: .. math:: \underline \sigma = (\sigma_{11}, \sigma_{22}, \sigma_{33}, \sqrt{2}\sigma_{12}, \sqrt{2}\sigma_{13}, \sqrt{2}\sigma_{23}). Integral operators ------------------ Integral operators are a central part of Tamaas: they are carefully designed for performance in periodic system. When a :py:class:`Model ` object is used with contact solvers or with a residual object (for plasticty), the objects using the model register integral operators with the model, so the user typically does not have to worry about creating integral operators. Integral operators are accessed through the :py:attr:`operators ` property of a model object. The ``[]`` operator gives access to the operators, and ``list(model.operators)`` gives the list of registered operators:: # Accessing operator elasticity = model.operators['hooke'] # Applying operator elasticity(strain, stress) # Print all registered operators print(list(model.operators)) .. note:: At model creation, these operators are automatically registered: - ``hooke``: Hooke's elasticity law - ``von_mises``: computes Von Mises stresses - ``deviatoric``: computes the deviatoric part of a stress tensor - ``eigenvalues``: computes the eigenvalues of a symetric tensor field :cpp:class:`Westergaard ` operators are automatically registered when :py:meth:`solveNeumann ` or :py:meth:`solveDirichlet ` are called. Additional operators ~~~~~~~~~~~~~~~~~~~~ Tamaas defines operators beyond the default ones listed above, but does not instanciate them by default becuse they require additional storage. Most are automatically instanciated as needed (e.g. with :cpp:class:`Residual `), but the :py:class:`ModelFactory ` class can be used to register these additional operators. .. tip:: Some operators define their own fields, just like a model. For example, the Fourrier coefficients are a field of the Westergaard operator. :cpp:class:`HookeField ` ............................................ The :cpp:class:`HookeField ` operator defines a linear elastic constitutive law with space-dependent property, which is assigned *via* the ``mu`` field belonging to the operator:: # Register the operator with a name tm.ModelFactory.registerHookeField(model, 'heterogeneous') # Function defining the shear modulus mu = lambda z, x, y: ... # Assign the shear modulus values model.operators['heterogeneous']['mu'][:] = mu(z, x, y) Non-periodic surface operator ............................. The Westergaard operator computes equilibrium surface displacement from normal pressure for a periodic system, but a non-periodic equilibrium computation is possible with the DC-FFT method. The operator is registered with :py:meth:`ModelFactory.registerNonPeriodic `. Volume operators ................ When the model type is ``volume_2d``, two volume-based half-space operators are available (with their gradients): the Boussinesq and the Mindlin operators. The Boussinesq operator acts on a traction field defined on the boundary (i.e. the half-space surface) and constructs a volumetric displacement field in equilibrium with this traction field :math:`\mathbf t`, with a divergence-free stress. The operator is often referred-to with the symbol :math:`\mathcal{M}`, and satisfies: .. math:: \mathrm{div}(\mathcal C : \varepsilon[\mathcal{M}[\mathbf t]]) = \mathbf 0,\quad\text{with}\quad \left.(\mathcal C:\varepsilon[\mathcal M[\mathbf t]])\cdot \mathbf n\right|_\text{surface} = \mathbf t The Mindlin operator acts on a volumetric eigenstress field and constructs a displacement field in equilibrium with the eigenstress and with a free surface (i.e. zero traction on the boundary). The operator is often referred to with the symbol :math:`\mathcal{N}`, and satisfies: .. math:: \mathrm{div}(\mathcal C : \varepsilon[\mathcal{N}[\tau]]) = -\mathrm{div}(\tau),\quad\text{with}\quad \left.(\mathcal C:\varepsilon[\mathcal N[\tau]])\cdot \mathbf n\right|_\text{surface} = \mathbf 0 The gradient of these operators computes the strain instead of displacement. These fields are instantiated with :py:meth:`ModelFactory.registerVolumeOperators `. They have the following names: - ``boussinesq``, ``boussinesq_gradient`` - ``mindlin``, ``mindlin_gradient`` Scipy Sparse interoperability ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Scipy defines an interface for matrix-vector products with the class :py:class:`scipy.sparse.linalg.LinearOperator`. This interface is implemented by several operators in Tamaas, and can be used easily with the function :py:func:`scipy.sparse.linalg.aslinearoperator`:: from scipy.sparse.linalg import aslinearoperator ɛ = np.ravel(model['strain']) # make one-dimensional C = aslinearoperator(model.operators['hooke']) σ = C * ɛ This makes these linear operators easily composable (with the operators ``+, -, *``) and interoperable with Scipy's sparse linear solvers. 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. Saves separate files for each dump of a model. :py:class:`NetCDFDumper ` Dumps a model to a compressed `NetCDF `_ file. Requires the `netCDF4 `_ package. Saves sequential dumps of a model into a single file, with the ``frame`` dimension containing the model dumps. 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 trigger all 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 .. important:: Currently, only :py:class:`H5Dumper ` and :py:class:`NetCDFDumper ` support parallel output with MPI.