Build Displacement Operator#

build_displacement_operator(jacobian_dxyz, shape_function, element_connectivity, element_indices, vertices_number=None, *, sparse=False, skip_m1=True, default=0.0)[source]#

Build the displacement operator that maps nodal displacements to image displacements at integration points.

In a space of dimension \(E\) with a mesh constituted of \(N_{e}\) elements and \(N_{v}\) nodes, The mesh is composed of \(K\)-dimensional elements (with \(K \leq E\)) defined by \(N_{vpe}\) nodes for each element.

Lets consider \(M\) integration points located within elements. We have a jacobian array of shape \((M, P, E)\) that describes how displacements in the \(E\)-dimensional space affect each property at the integration points, where \(P\) (e.g., gray levels). We want to build the displacement operator that maps nodal displacements from property displacements at these integration points with shape \((M \times P, N_{v} \times E)\).

The objectif is to build the jacobian operator \(\mathbf{J} \in \mathbb{R}^{M \times (N_{v} \cdot E)}\) such that:

\[J dU = R\]

where \(dU \in \mathbb{R}^{N_{v} \cdot E}\) is the vector of nodal displacements and \(R \in \mathbb{R}^{M \times P}\) is the residual vector.

Note

For (\(P \neq 1\)), the lines are groupes by property first (i.e., all integration points for property 1, followed by all integration points for property 2, etc.). Thus the residual \(R\) associated with the output must be constructed accordingly.

\[\begin{split}R = \begin{bmatrix} R_{1}^{1} \\ R_{2}^{1} \\ \vdots \\ R_{M}^{1} \\ R_{1}^{2} \\ R_{2}^{2} \\ \vdots \\ R_{M}^{2} \\ \vdots \\ R_{1}^{P} \\ R_{2}^{P} \\ \vdots \\ R_{M}^{P} \end{bmatrix}\end{split}\]

where \(R_{i}^{j}\) represents the residual at integration point \(i\) for property \(j\).

residual = numpy.random.rand(M, P)
residual_vector = residual.flatten(order='F')  # Column-major flattening

Note

The columns of the operator correspond to the nodal displacements organized by direction first (i.e., all nodes for direction 1, followed by all nodes for direction 2, etc.). Thus the displacement vector \(dU\) must be constructed accordingly.

\[\begin{split}dU = \begin{bmatrix} dU_{1}^{x} \\ dU_{2}^{x} \\ \vdots \\ dU_{N_{v}}^{x} \\ dU_{1}^{y} \\ dU_{2}^{y} \\ \vdots \\ dU_{N_{v}}^{y} \\ \vdots \\ dU_{1}^{E} \\ dU_{2}^{E} \\ \vdots \\ dU_{N_{v}}^{E} \end{bmatrix}\end{split}\]

where \(dU_{i}^{j}\) represents the displacement of node \(i\) in direction \(j\).

displacement_vector = numpy.random.rand(N_v * E) # resulting vector of size N_v * E
displacement = displacement_vector.reshape((N_v, E), order='F')  # Column-major reshaping

Note

The inputs jacobian_dxyz, and shape_function will be converted to numpy.float64 The inputs element_connectivity and element_indices will be converted to numpy.int64

Warning

No tests are performed to check the consistency of the input data (e.g., whether the element indices are valid, …). Only the shapes and types of the input arrays are checked. The behavior of the function is undefined if the input data is not consistent.

Parameters:

  • jacobian_dxyz (ArrayLike of shape \((M, E)\) or \((M, P, E)\)) – Array containing the jacobian of the image properties with respect to displacements in the \(E\)-dimensional space at \(M\) integration points.

  • shape_function (ArrayLike of shape \((M, N_{vpe})\)) – Array containing the shape function values evaluated at \(M\) points for the \(N_{vpe}\) nodes of the element.

  • element_connectivity (ArrayLike of shape \((N_{e}, N_{vpe})\)) – Array defining the connectivity of the elements in the mesh, where each row contains the indices of the nodes that form an element.

  • element_indices (ArrayLike of shape \((M,)\)) – An array containing the indices of each element corresponding to the \(M\) integration points.

  • vertices_number (Optional[Integral], optional) – Total number of vertices \(N_{v}\) in the mesh. If not provided, it will be inferred from the element_connectivity.

  • sparse (bool, optional) – If set to True, the function will return a sparse matrix representation of the displacement operator. Default is False.

  • skip_m1 (bool, optional) – If set to True, any element index of -1 in element_indices will result in the corresponding rows in the displacement operator being set to default. Default is True.

  • default (Real, optional) – The default value to use for entries corresponding to integration points with element index -1 when skip_m1 is True. Default is 0.0.

Returns:

displacement_operator – The displacement operator mapping nodal displacements from property displacements at integration points.

Return type:

Union[numpy.ndarray, scipy.sparse.csr_matrix] of shape \((M \times P, N_{v} \times E)\)

Raises:

ValueError – If the input arrays do not have the expected shapes or types, or if the parameters are not valid.