pysdic.assemble_property_derivative#

assemble_property_derivative(property_array, connectivity, element_indices, shape_function_derivatives, jacobian_matrix, *, skip_m1=True, default=nan)[source]#

Assemble the derivative of a property defined at the nodes of a mesh to given integration points within elements with respect to global coordinates \((x,y,z,...)\) from precomputed shape function derivatives and Jacobian matrices.

derivated_properties[i] contains the derivated property values at the integration point \(i\).

Note

  • Inputs property_array, shape_function_derivatives and jacobian_matrix will be converted to numpy.float64.

  • Inputs connectivity and element_indices will be converted to numpy.int64.

Warning

No tests are performed to check if the inputs array are consistent (e.g., if the connectivity contains invalid vertex indices, …). Only shapes are validated. The behavior of the function is undefined in this case.

Important

When using -1 in element_indices for invalid elements, ensure to set skip_m1 to True to avoid indexing errors.

Parameters:

  • property_array (ArrayLike) – An array of shape \((N_{v}, P)\) containing the property values defined at the nodes of the mesh. If 1D-array is provided, it will be treated as a single-component property of shape \((N_{v}, 1)\).

  • connectivity (ArrayLike) – An array of shape \((N_{e}, N_{vpe})\) 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) – An array of shape \((N_{p},)\) containing the indices of each element corresponding to the \(N_{p}\) integration points.

  • shape_function_derivatives (ArrayLike) – An array of shape \((N_{p}, N_{vpe}, K)\) containing the derivatives of the shape functions evaluated at \(N_{p}\) points for the \(N_{vpe}\) nodes of the element.

  • jacobian_matrix (ArrayLike) – An array of shape \((N_{p}, E, K)\) containing the Jacobian matrices for each of the \(N_{p}\) points. Each Jacobian matrix has dimensions \((E \times K)\).

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

  • default (Real, optional) – The default value to assign to derivated properties when the Jacobian matrix is singular or not invertible. Default is numpy.nan.

Returns:

derivated_properties – An array of shape \((N_{p}, P, E)\) containing the derivated property values at each of the \(N_{p}\) integration points with respect to the global coordinates \((x, y, z, ...)\).

Return type:

numpy.ndarray

Notes

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.

For a given set of \(N_p\) integration points located within elements, the derivated property array has shape \((N_p, P, K)\) where \(P\) is the number of property components (e.g., 1 for scalar properties, 3 for vector properties), and \(K\) is the number of local coordinates.

The derivative of the property at each integration point is computed as:

If the dimension of the elements \(K\) is equal to the space dimension \(E\):

\[\nabla_X P(\xi, \eta, \zeta, ...) = J^{-T} \cdot \sum_{i=1}^{N_{vpe}} \nabla_{\Xi} N_i(\xi, \eta, \zeta, ...) P_i\]

If the dimension of the elements \(K\) is less than the space dimension \(E\):

\[\nabla_X P(\xi, \eta, \zeta, ...) = J (J^T J)^{-1} \cdot \sum_{i=1}^{N_{vpe}} \nabla_{\Xi} N_i(\xi, \eta, \zeta, ...) P_i\]

where \(N_i\) are the shape functions associated with each node, \(P_i\) are the property values at the nodes of the element, and \(J\) is the Jacobian matrix of the transformation from natural to global coordinates.

Demonstration#

The property field is defined at the nodes of the mesh and can be interpolate to any point within the elements using shape functions.

\[P(\Xi) = \sum_{i=1}^{N_{vpe}} N_i(\Xi) P_i\]

The spatial derivative of the property field with respect to global coordinates is defined as:

\[\nabla_X P(\Xi) = \sum_{i=1}^{N_{vpe}} \nabla_{X} N_i(\Xi) P_i\]

The derivative of the shape functions with respect to global coordinates is obtained via the chain rule:

\[\nabla_{X} N = J_{\Xi/X}^T \nabla_{\Xi} N\]

where \(J_{\Xi/X}\) is the Jacobian matrix of the transformation from global coordinates to natural coordinates, which is the inverse of the Jacobian matrix \(J = J_{X/\Xi}\) of the transformation from natural coordinates to global coordinates.

For \(K = E\), the Jacobian matrix is square and invertible :

\[J_{\Xi/X} = J^{-1} \rightarrow \nabla_{X} N = J^{-T} \nabla_{\Xi} N\]

And we obtain the formula as mentioned above:

\[\nabla_X P(\Xi) = J^{-T} \cdot \sum_{i=1}^{N_{vpe}} \nabla_{\Xi} N_i(\Xi) P_i\]

For \(K < E\), the Jacobian matrix is not square and thus not invertible. The pseudo-inverse is :

\[J_{\Xi/X} = (J^T J)^{-1} J^T \rightarrow \nabla_{X} N = J (J^T J)^{-1} \nabla_{\Xi} N\]

And we obtain the formula as mentioned above:

\[\nabla_X P(\Xi) = J (J^T J)^{-1} \cdot \sum_{i=1}^{N_{vpe}} \nabla_{\Xi} N_i(\Xi) P_i\]

See also

pysdic.compute_shape_functions

To compute the shape functions and their derivatives at given natural coordinates within elements.

pysdic.assemble_jacobian_matrix

To assemble the Jacobian matrices from shape function derivatives and remapped coordinates.

pysdic.compute_property_derivative

To compute the property derivatives directly from mesh and integration point information.

Examples

Lets construct a simple 2D mesh and compute the derivative of a scalar property at given integration points for triangular elements.

 1import numpy
 2from pysdic import (
 3    remap_vertices_coordinates,
 4    assemble_jacobian_matrix,
 5    assemble_property_derivative,
 6    compute_triangle_3_shape_functions
 7)
 8
 9vertices_coordinates = numpy.array(
10    [[0.0, 0.0],
11     [1.0, 0.0],
12     [1.0, 1.0],
13     [0.0, 1.0]]
14)
15
16connectivity = numpy.array(
17    [[0, 1, 2],
18     [0, 2, 3]]
19)
20
21natural_coordinates = numpy.array(
22    [[0.2, 0.3],
23     [0.6, 0.2]]
24)
25
26element_indices = numpy.array([0, 1])
27
28property_array = numpy.array([10.0, 20.0, 30.0, 40.0])  # Scalar property
29
30remapped_coords = remap_vertices_coordinates(
31    vertices_coordinates=vertices_coordinates,
32    connectivity=connectivity,
33    element_indices=element_indices
34)
35
36_, shape_func_derivs = compute_triangle_3_shape_functions(
37    natural_coordinates=natural_coordinates,
38    return_derivatives=True
39)
40
41jacobians = assemble_jacobian_matrix(
42    shape_function_derivatives=shape_func_derivs,
43    remapped_coordinates=remapped_coords
44)
45
46derivated_properties = assemble_property_derivative(
47    property_array=property_array,
48    connectivity=connectivity,
49    element_indices=element_indices,
50    shape_function_derivatives=shape_func_derivs,
51    jacobian_matrix=jacobians
52)
53
54print(f"derivated properties (shape={derivated_properties.shape}):")
55print(derivated_properties)

derivated properties (shape=(2, 1, 2)):
[[[ 10.  10.]]

 [[-10.  30.]]]