Python Bindings¶
Fields¶
The µGrid library handles field quantities, i.e. scalar, vectors or tensors, that vary in space. It supports only a uniform discretization of space. In the language of µGrid, we call the discrete coordinates pixels. Each pixel is associated with a physical position in space. µGrid supports fields in two and three dimensional Cartesian grids. Note that a common name for a pixel in three dimensions is a voxel, but we refer to it as a pixel throughout this documentation.
Each pixel can be subdivided into logical elements, for example into a number of support points for numerical quadrature. These subdivisions are called sub-points. Note that while µGrid understands the physical location of each pixel, the sub-points are logical subdivisions and not associated with any position within the pixel. Each sub-point carries the field quantity, which can be a scalar, vector or any type of tensor. The field quantity is called the component.
Multidimensional arrays¶
Each µGrid field has a representation in a multidimensional array. This representation is used when accessing a field with Python via numpy or when writing the field to a file. The default representation has the shape
(components, sub-points, pixels)
As a concrete example, a second-rank tensor (for example the deformation gradient) living on two quadrature points in three dimensions with a spatial discretization of 11 x 12 grid points would have the following shape:
(3, 3, 2, 11, 12)
Note that the components are omitted if there is only a single component (i.e. a scalar value), but the sub-points are always included even if there is only a single sub-point.
Throughout the code, we call the field quantities entries. For example, nb_pixels refers to the number of pixels while nb_sub_pts refers to the number of sub-points per pixel. nb_entries is then the product of the two.
The above field has another multidimensional array representation that folds the sub-point into the last dimension of the components. For the above example, this means a multidimensional array of shape:
(3, 6, 11, 12)
Depending on the numerical use case, it can be useful to work with either representation.
Field collections¶
In µGrid, fields are grouped into field collections. A field collection knows about the spatial discretization of the fields, but each field can differ in number of sub-points and components.
There are two kinds of field collections:
A global collection that groups fields that have values at all pixels.
A local collection that groups fields that have values only at a subset of pixels. An example would be the elastic constants of a material that only exists at certain locations in the domain.
The following example shows how to initialize a field collection and create scalar fields:
import numpy as np
from muGrid import GlobalFieldCollection
# Two dimensional grid
nb_grid_pts = (11, 12)
fc = GlobalFieldCollection(nb_grid_pts)
# Get a scalar (real-valued) field
# This call creates a new field if it does not exist.
field = fc.real_field('my-real-valued-field')
# Check the registered fields
assert fc.field_names == ['my-real-valued-field']
# Fill the field with random numbers
field.p = np.random.rand(*nb_grid_pts)
# Set a specific value
field.p[5, 6] = 42
# Get the field, but now into a different variable
field2 = fc.real_field('my-real-valued-field')
# Check that the value is the same
assert field2.p[5, 6] == 42
The first argument to the constructor of FieldCollection is the spatial dimension. Fields are then registered with the field accessor methods of the field collection, e.g. real_field in the example above. Fields are named, and the name needs to be unique. Accessing a field of the same name yield the same field object.
Note the FieldCollection additionally has register methods, e.g. register_real_field. The different to real_field method is that the explicit registration of the field fails if it already exists, while the accessor method real_field registers it if it does not exist but returns the existing field if it does. You can also query a field with get_field, which will raise an exception if the field does not exist.
Components¶
In the above example, we registered a scalar field, which has one component. Vector or tensor-valued field can be defined by specifying either simply a number of components or the shape of the components explicitly. The following example shows how to create a tensor-valued field that contains 2 x 2 matrices:
import numpy as np
from muGrid import GlobalFieldCollection
# Two dimensional grid
nb_grid_pts = (11, 12)
fc = GlobalFieldCollection(nb_grid_pts)
# Get a tensor-field (for example to represent the strain)
strain = fc.real_field('strain', (2, 2))
# Fill the field with random numbers
strain.p = np.random.rand(*((2, 2) + nb_grid_pts))
Sub-points¶
A pixel can be subdivided into multiple sub-points, each of which holds a value (scalar, vector or tensor) of the field. Examples for sub-points are elements or quadrature points in a the finite element method.
Sub-points are named, e.g. common names would be element for subdivision into elements or quad for quadrature points. The name of the subdivision must be specified when the field is created.
The following example initializes a three-dimension grid with a sub-division of each pixel into five elements:
import numpy as np
from muGrid import GlobalFieldCollection
# Two dimensional grid
nb_grid_pts = (11, 12, 13)
fc = GlobalFieldCollection(nb_grid_pts, sub_pts={'element': 5})
# Get a tensor-field (for example to represent the strain)
strain = fc.real_field(
'strain', # name of the field
(3, 3), # shape of components
'element' # sub-point type
)
# Fill the field with random numbers
strain.s = np.random.rand(*((3, 3, 5) + nb_grid_pts))
# Note that the covenience accessor `p` yields the same field but a different shape
strain.p = np.random.rand(*((3, 3 * 5) + nb_grid_pts))
The example demonstrates two ways of accessing the field. The convenience accessor p (that we also used in the above examples) and the new accessor s. Both yield a numpy array that is a view on the underlying data, but with different shapes. The s accessor has the shape (components, sub-points, pixels) and exposes the sub-points explicitly. The p accessor folds the sub-points into the last dimension of the components.
numpy views¶
The multidimensional array representation of a field is accessible via the array method.
a = displacement_field.array(muGrid.IterUnit.SubPt)
yields the multidimensional array with the explicit sub-point dimension. The pixel-representation can be obtained by
a = displacement_field.array(muGrid.IterUnit.Pixel)
Because those operations are used to frequently, there are shortcuts already introduced in the examples above:
displacement_field.s # sub-point representation
displacement_field.p # pixel representation
The entries of the field occur as the first indices in the multidimensional because a numerical code is typically vectorized of the spatial domain, i.e. we carry out the same operation on each pixel but not on each component. This means for the above displacement field, we can simply get the components of the field from
ux, uy, uz = displacement_field.s
Each of the variable ux, uy, and uz is a three-dimensional array with shape (2, 11, 12). numpy’s broadcasting rules make it simple to vectorize over pixels, for example normalizing the displacement field could look like:
displacement_field.s /= np.sqrt(ux**2 + uy**2 + uz**2)
Note that the default storage order of the field is column-major, which means the field components are stored next to each other in memory.
I/O¶
Fields can be written to disk in the NetCDF format. µGrid uses Unidata NetCDF when build with just serial capabilities and PnetCDF when build with MPI enabled.
I/O is handled by the FileIONetCDF class. The following example shows how to write all fields from a field collection to disk:
import numpy as np
from muGrid import GlobalFieldCollection, FileIONetCDF, OpenMode
# Two dimensional grid
nb_grid_pts = (11, 12, 13)
fc = GlobalFieldCollection(nb_grid_pts, sub_pts={'element': 5})
# Get a tensor-field (for example to represent the strain)
strain = fc.real_field(
'strain', # name of the field
(3, 3), # shape of components
'element' # sub-point type
)
# Fill the field with random numbers
strain.s = np.random.rand(*((3, 3, 5) + nb_grid_pts))
# Initialize a file I/O object
file = FileIONetCDF('example.nc', OpenMode.Overwrite)
file.register_field_collection(fc) # Register the field collection with the file
file.append_frame().write() # Write all fields of the field collection to the file
The file has the following structure (output of ncdump -h):
netcdf example {
dimensions:
frame = UNLIMITED ; // (1 currently)
tensor_dim__strain-0 = 3 ;
tensor_dim__strain-1 = 3 ;
subpt__element-5 = 5 ;
nx = 11 ;
ny = 12 ;
nz = 13 ;
variables:
double strain(frame, tensor_dim__strain-0, tensor_dim__strain-1, subpt__element-5, nx, ny, nz) ;
strain:unit = "no unit provided" ;
// global attributes:
:creation_date = "01-06-2024 (d-m-Y)" ;
:creation_time = "23:02:06 (H:M:S)" ;
:last_modified_date = "01-06-2024 (d-m-Y)" ;
:last_modified_time = "23:02:06 (H:M:S)" ;
:muGrid_version_info = "µGrid version: 0.90.1+40-g5bc73b30\n",
"" ;
:muGrid_git_hash = "5bc73b305881ef837ce6598568d41eb3d0307c41" ;
:muGrid_description = "0.90.1+40-g5bc73b30" ;
:muGrid_git_branch_is_dirty = "false" ;
}