µGrid

As the lowest-level component of µSpectre, µGrid defines all the commonly used type aliases and data structures used througout the project. The most common aliases are described below, but it is worth having a look at the file grid_common.hh for the details.

Common Type Aliases

All mathematical calculations should use the types.

Warning

doxygengroup: Cannot find group “Scalars” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

While it is possible to use other types in principle, these are the ones for which all datastructures are tested and which are known to work. Also, other µSpectre developpers will expect and understand these types.

Dimensions are counted using the signed integer type muGrid::Dim_t. This is necessary because Eigen uses -1 to signify a dynamic number of dimensions.

The types muGrid::Rcoord_t and muGrid::Ccoord_t are used to represent real-valued coordinates and integer-valued coordinates (i.e., pixel- or cell-coordinates).

group Coordinates

Typedefs

using Ccoord_t = std::array<Dim_t, Dim>

Ccoord_t are cell coordinates, i.e. integer coordinates.

using Rcoord_t = std::array<Real, Dim>

Real space coordinates.

using DynCcoord_t = DynCcoord<threeD>

usually, we should not need omre than three dimensions

using DynRcoord_t = DynCcoord<threeD, Real>

usually, we should not need omre than three dimensions

Functions

template<typename T, size_t Dim>
Eigen::Map<Eigen::Matrix<T, Dim, 1>> eigen(std::array<T, Dim> &coord)

return a Eigen representation of the data stored in a std::array (e.g., for doing vector operations on a coordinate)

template<typename T, size_t Dim>
Eigen::Map<const Eigen::Matrix<T, Dim, 1>> eigen(const std::array<T, Dim> &coord)

return a constant Eigen representation of the data stored in a std::array (e.g., for doing vector operations on a coordinate)

template<typename T, size_t MaxDim>
Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>> eigen(DynCcoord<MaxDim, T> &coord)

return a Eigen representation of the data stored in a std::array (e.g., for doing vector operations on a coordinate)

template<typename T, size_t MaxDim>
Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>> eigen(const DynCcoord<MaxDim, T> &coord)

return a const Eigen representation of the data stored in a std::array (e.g., for doing vector operations on a coordinate)

template<size_t MaxDim, typename T = Dim_t>
class DynCcoord

#include <grid_common.hh>

Class to represent integer (cell-) coordinates or real-valued coordinates. This class can dynamically accept any spatial-dimension between 1 and MaxDim, and DynCcoord references can be cast to muGrid::Ccoord_t & or muGrid::Rcoord_t & references. These are used when templating with the spatial dimension of the problem is undesireable/impossible.

Public Types

using iterator = typename std::array<T, MaxDim>::iterator

iterator type

using const_iterator = typename std::array<T, MaxDim>::const_iterator

constant iterator type

Public Functions

inline DynCcoord()

default constructor

inline DynCcoord(std::initializer_list<T> init_list)

constructor from an initialiser list for compound initialisation.

Parameters

init_list – The length of the initialiser list becomes the spatial dimension of the coordinate, therefore the list must have a length between 1 and MaxDim

inline explicit DynCcoord(Dim_t dim)

Constructor only setting the dimension. WARNING: This constructor needs regular (round) braces ‘()’, using curly braces ‘{}’ results in the initialiser list constructor being called and creating a DynCcoord with spatial dimension 1

Parameters

dim – spatial dimension. Needs to be between 1 and MaxDim

template<size_t Dim>
inline explicit DynCcoord(const std::array<T, Dim> &ccoord)

Constructor from a statically sized coord.

DynCcoord(const DynCcoord &other) = default

Copy constructor.

DynCcoord(DynCcoord &&other) = default

Move constructor.

~DynCcoord() = default

nonvirtual Destructor

template<size_t Dim>
inline DynCcoord &operator=(const std::array<T, Dim> &ccoord)

Assign arrays.

DynCcoord &operator=(const DynCcoord &other) = default

Copy assignment operator.

DynCcoord &operator=(DynCcoord &&other) = default

Move assignment operator.

template<size_t Dim2>
inline bool operator==(const std::array<T, Dim2> &other) const

comparison operator

inline bool operator==(const DynCcoord &other) const

comparison operator

template<typename T2>
inline DynCcoord<MaxDim, decltype(T{} / T2{})> operator/(const DynCcoord<MaxDim, T2> &other) const

element-wise division

inline T &operator[](const size_t &index)

access operator

inline const T &operator[](const size_t &index) const

access operator

template<size_t Dim>
inline operator std::array<T, Dim>() const

conversion operator

template<Dim_t Dim>
inline std::array<T, Dim> &get()

cast to a reference to a statically sized array

template<Dim_t Dim>
inline const std::array<T, Dim> &get() const

cast to a const reference to a statically sized array

inline const Dim_t &get_dim() const

return the spatial dimension of this coordinate

inline iterator begin()

iterator to the first entry for iterating over only the valid entries

inline iterator end()

iterator past the dim-th entry for iterating over only the valid entries

inline const_iterator begin() const

const iterator to the first entry for iterating over only the valid entries

inline const_iterator end() const

const iterator past the dim-th entry for iterating over only the valid entries

inline T *data()

return the underlying data pointer

inline const T *data() const

return the underlying data pointer

inline T &back()

return a reference to the last valid entry

inline const T &back() const

return a const reference to the last valid entry

These types are also used to define nb_grid_pts or spatial lengths for computational domains.

Field Data Types

The most important part of µGrid to understand is how it handles field data and access to it. By field we mean the discretisation of a mathematical field on the grid points, i.e., numerical data associated with all pixels/voxels of an FFT grid or a subset thereof. The numerical data can be scalar, vectorial, matricial, tensorial or a generic array of integer, real or complex values per pixel.

Fields that are defined on every pixel/voxel of a grid are called global fields while fields defined on a subset of pixels/voxels local fields. As an example, the strain field is a global field for any calculation (it exists in the entire domain), while for instance the field of an internal (state) variable of a material in a composite is only defined for the pixels that belong to that material.

There are several ways in which we interact with fields, and the same field might be interacted with in different ways by different parts of a problem. Let’s take the (global) strain field in a three-dimensional finite strain problem with 255 × 255 × 255 voxels as an example: the solver treats it as a long vector (of length 3² · 255³), the FFT sees it as a four-dimensional array of shape 255 × 255 × 255 × 3², and from the constitutive laws’ perspective, it is just a sequence of second-rank tensors (i.e., shape 255³ × 3 × 3).

Basic µGrid Field Concepts

In order to reconcile these different interpretations without copying data around, µGrid splits the concept of a field into three components:

• storage

  This refers managing the actual memory in which field data is held. For this, the storage abstraction needs to know the scalar type of data (Int, Real, Complex, e.t.c.), the number of pixels/voxels for which the field is defined, and the number of scalar components per pixel/voxel (e.g., 9 for a second-rank asymmetric tensor in a three-dimensional problem).

  µGrid’s abstraction for field data storage is the field represented by a child class of FieldBase<FieldCollection_t>, see fields.

• representation

  Meaning how to interpret the data at a given pixel/voxel (i.e., is it a vector, a matrix, …). This will also determine which mathematical operations can be performed on per-pixel/voxel data. The representation allows also to iterate over a field pixel/voxel by pixel/voxel.

  µGrid’s abstraction for field representations is the field map represented by a child class of FieldMap<FieldCollection_t, Scalar_t, NbComponents[, IsConst]>, see field_map.

• per-pixel/voxel access/iteration

  Given a pixel/voxel coordinate or index, the position of the associated pixel/voxel data is a function of the type of field (global or local). Since the determination procedure is identical for every field defined on the same domain, this ability (and the associated overhead) can be centralised into a manager of field collections.

  µGrid’s abstraction for field access and management is the field collection represented by the two classes LocalFieldCollection<Dim> and GlobalFieldCollection<Dim>, see field_collection.

Fields

Fields are where the data is stored, so they are mainly distinguished by the scalar type they store (Int, Real or Complex), and the number of components (statically fixed size, or dynamic size).

The most commonly used fields are the statically sized ones, TensorField, MatrixField, and the ScalarField (which is really just a 1×1 matrix field).

Less commonly, we use the dynamically sized TypedField, but more on this later.

Fields have a protected constructor, which means that you cannot directly build a field object, instead you have to go through the factory function make_field<Field_t>(name, collection) (or make_statefield<Field_t>(name, collection) if you’re building a statefield, see state_field) to create them and you only receive a reference to the built field. The field itself is stored in a std::unique_ptr which is registered in and managed by a field collection. This mechanism is meant to ensure that fields are not copied around or free’d so that field maps always remain valid and unambiguously linked to a field.

Fields give access to their bulk memory in form of an Eigen:Map. This is useful for instance for accessing the global strain, stress, and tangent moduli fields in the solver.

Example: Using fields as global arrays:

The following is a code example from the standard Cell:

template <Dim_t DimS, Dim_t DimM>
auto CellBase<DimS, DimM>::get_strain_vector() -> Vector_ref {
  return this->get_strain().eigenvec();
}

The return value of :cpp:function:`Cell::get_strain_vector()<muSpectre::CellBase::get_strain_vector>` is an Eigen::Map onto a matrix of shape 1×N.

If you wish to handle field data on a per-pixel/voxel basis, the mechanism for that is the field map and is described in field_map.

Field Maps

Field maps are light-weight resource handles (meaning they can be created and destroyed cheaply) that are iterable and provide direct per-pixel/voxel access to the data stored in the mapped field.

The choice of field map defines the type of reference you obtain when dereferencing an iterator or using the direct random acccess operator [].

Typically used field maps include:

• ScalarFieldMap,

• ArrayFieldMap,

• MatrixFieldMap, and the

• T4MatrixFieldMap.

All of these are fixed size (meaning their size is set at compile time) and therefore support fast linear algebra on the iterates. There is also a dynamically sized field map type, the TypedFieldMap which is useful for debugging and python bindings. It supports all the features of the fixed-size maps, but linear algebra on the iterates will be slow because it cannot be effectively vectorised.

Example 1: Iterating over fields and do math on the iterates:

The following is a code example from the tests of the finite-strain projection operator defined by T.W.J. de Geus, J. Vondřejc, J. Zeman, R.H.J. Peerlings, M.G.D. Geers.

for (auto && tup :
     akantu::zip(fields.get_pixels().template get_dimensioned_pixels<dim>(),
                 grad, var)) {
  auto & ccoord = std::get<0>(tup);  // iterate from fields
  auto & g = std::get<1>(tup);       // iterate from grad
  auto & v = std::get<2>(tup);       // iterate from var

  // use iterate in arbitrary expressions
  Vector vec = muGrid::CcoordOps::get_vector(
      ccoord, (fix::projector.get_domain_lengths() /
               fix::projector.get_nb_domain_grid_pts())
                  .template get<dim>());
  // do efficient linear algebra on iterates
  g.row(0) = k.transpose() *
             cos(k.dot(vec));  // This is a plane wave with wave vector k in
                               // real space. A valid gradient field.
  v.row(0) = g.row(0);
}

Field Collections

Field collections come in two flavours; LocalFieldCollection<Dim> and GlobalFieldCollection<Dim> and are templated by the spatial dimension of the problem. They adhere to the interface defined by their common base class, FieldCollectionBase. Both types are iterable (the iterates are the coordinates of the pixels/voxels for which the fields of the collection are defiened.

Global field collections need to be given the problem nb_grid_pts (i.e. the size of the grid) at initialisation, while local collections need to be filled with pixels/voxels through repeated calls to add_pixel(pixel). At initialisation, they derive their size from the number of pixels that have been added.

Fields (State Fields) are identified by their unique name (prefix) and can be retrieved in multiple ways:

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::operator[]” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::at” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_typed_field” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_statefield” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_typed_statefield” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

• per-pixel/voxel access/iteration

  Given a pixel/voxel coordinate or index, the position of the associated pixel/voxel data is a function of the type of field (global or local). Since the determination procedure is identical for every field defined on the same domain, this ability (and the associated overhead) can be centralised into a manager of field collections.

  µGrid’s abstraction for field access and management is the field collection represented by the two classes LocalFieldCollection<Dim> and LocalFieldCollection<Dim>, see field_collection.

Fields

Fields are where the data is stored, so they are mainly distinguished by the scalar type they store (Int, Real or Complex), and the number of components (statically fixed size, or dynamic size).

The most commonly used fields are the statically sized ones, TensorField, MatrixField, and the ScalarField (which is really just a 1×1 matrix field).

Less commonly, we use the dynamically sized TypedField, but more on this later.

Fields have a protected constructor, which means that you cannot directly build a field object, instead you have to go through the factory function make_field<Field_t>(name, collection) (or make_statefield<Field_t>(name, collection) if you’re building a statefield, see state_field) to create them and you only receive a reference to the built field. The field itself is stored in a std::unique_ptr which is registered in and managed by a field collection. This mechanism is meant to ensure that fields are not copied around or free’d so that field maps always remain valid and unambiguously linked to a field.

Fields give access to their bulk memory in form of an Eigen:Map. This is useful for instance for accessing the global strain, stress, and tangent moduli fields in the solver.

If you wish to handle field data on a per-pixel/voxel basis, the mechanism for that is the field map and is described in field_map.

Example: Using fields as global arrays:

The following is a code example from the standard Cell:

template <Dim_t DimS, Dim_t DimM>
auto CellBase<DimS, DimM>::get_strain_vector() -> Vector_ref {
  return this->get_strain().eigenvec();
}

The return value of :cpp:function:`Cell::get_strain_vector()<muSpectre::CellBase::get_strain_vector>` is an Eigen::Map onto a matrix of shape 1×N.

If you wish to handle field data on a per-pixel/voxel basis, the mechanism for that is the field map and is described in field_map.

Field Maps

Field maps are light-weight resource handles (meaning they can be created and destroyed cheaply) that are iterable and provide direct per-pixel/voxel access to the data stored in the mapped field.

The choice of field map defines the type of reference you obtain when dereferencing an iterator or using the direct random acccess operator [].

Typically used field maps include:

• ScalarFieldMap,

• ArrayFieldMap,

• MatrixFieldMap, and the

• T4MatrixFieldMap.

All of these are fixed size (meaning their size is set at compile time) and therefore support fast linear algebra on the iterates. There is also a dynamically sized field map type, the TypedFieldMap which is useful for debugging and python bindings. It supports all the features of the fixed-size maps, but linear algebra on the iterates will be slow because it cannot be effectively vectorised.

Example 1: Iterating over fields and do math on the iterates:

The following is a code example from the tests of the finite-strain projection operator defined by T.W.J. de Geus, J. Vondřejc, J. Zeman, R.H.J. Peerlings, M.G.D. Geers.

for (auto && tup :
     akantu::zip(fields.get_pixels().template get_dimensioned_pixels<dim>(),
                 grad, var)) {
  auto & ccoord = std::get<0>(tup);  // iterate from fields
  auto & g = std::get<1>(tup);       // iterate from grad
  auto & v = std::get<2>(tup);       // iterate from var

  // use iterate in arbitrary expressions
  Vector vec = muGrid::CcoordOps::get_vector(
      ccoord, (fix::projector.get_domain_lengths() /
               fix::projector.get_nb_domain_grid_pts())
                  .template get<dim>());
  // do efficient linear algebra on iterates
  g.row(0) = k.transpose() *
             cos(k.dot(vec));  // This is a plane wave with wave vector k in
                               // real space. A valid gradient field.
  v.row(0) = g.row(0);
}

Field Collections

Field collections come in two flavours; LocalFieldCollection<Dim> and GlobalFieldCollection<Dim> and are templated by the spatial dimension of the problem. They adhere to the interface defined by their common base class, FieldCollectionBase. Both types are iterable (the iterates are the coordinates of the pixels/voxels for which the fields of the collection are defiened.

Global field collections need to be given the problem nb_grid_pts (i.e. the size of the grid) at initialisation, while local collections need to be filled with pixels/voxels through repeated calls to add_pixel(pixel). At initialisation, they derive their size from the number of pixels that have been added.

Fields (State Fields) are identified by their unique name (prefix) and can be retrieved in multiple ways:

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::operator[]” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::at” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_typed_field” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_statefield” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

Warning

doxygenfunction: Cannot find function “muGrid::FieldCollectionBase::get_typed_statefield” in doxygen xml output for project “µSpectre” from directory: ./doxygenxml

State or History Variables

Some fields hold state or history variables, i.e., such fields have a current value and one or more old values. This is particularly common for internal variables of inelastic materials (e.g., damage variables, plastic flow, e.t.c.). The straight-forward way of handling this situation is to define a current field, and one or more fields of the same type to hold old values. This approach has the disadvantages that it leads to a multitude of variables to keep track of, and that the values need to by cycled between the fields using a copy; this approach is both inefficient and error-prone.

µGrid addresses this situation with the state field abstraction. A state

field is an encapsulated container of fields in a single variable. It allows to access the current field values globally, and gives read-only access to old field values globally. Iterative per-pixel access is handled through state field maps which, similarly to the field map, allow to iterate though all pixels/voxels on which the field is defined, and the iterates give access to the current value at the pixel/voxel or read-only access to the old values.

Mapped Fields

Some fields are only ever going to be used by one entity (e.g., internal variables of a material). For these fields, the flexibility of the field/field collection/field map paradigm can be a burden. Mapped fields are an encapsulation of a field and a corresponding map into a single object, drastically reducing boilerplate code.