pysdic.assemble_property_projection#

assemble_property_projection(property_array, shape_function_matrix, point_weights=None, *, return_unaffected=False)[source]#

Project a property defined at integration points within elements back to the nodes of a mesh using a precomputed shape function matrix.

projected_properties[i] contains the projected property at the node \(i\).

Note

  • Inputs property_array, shape_function_matrix and point_weights will be converted to numpy.float64.

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.

Parameters:

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

  • shape_function_matrix (Union[ArrayLike, scipy.sparse.csr_matrix]) – A precomputed shape functions matrix of shape \((N_{p}, N_{v})\) where \(N_{v}\) is the number of nodes in the mesh.

  • point_weights (Optional[ArrayLike], optional) – An array of shape \((N_{p},)\) containing the weights associated with each integration point. If not provided, all weights will be assumed to be equal to one.

  • return_unaffected (bool, optional) – If set to True, the function will return a tuple containing the projected properties and a boolean mask indicating which nodes were unaffected by any integration point. Default is False.

Returns:

  • projected_properties (numpy.ndarray) – An array of shape \((N_{v}, P)\) containing the projected property values at the nodes of the mesh.

  • unaffected_mask (numpy.ndarray, optional) – A boolean array of shape \((N_{v},)\) indicating which nodes were unaffected by any integration point. This is only returned if return_unaffected is set to True. True indicates the node was unaffected, while False indicates it was affected. So using projected_properties[unaffected_mask, :] = default will set the unaffected nodes to default.

Return type:

ndarray | Tuple[ndarray, 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.

The evaluation of a property at the given integration points is represented by an array of shape \((N_{p}, P)\) and is performed as:

\[P_{points} = N_f \cdot P_{nodes} \Leftrightarrow P(\xi, \eta, \zeta, ...) = \sum_{i=1}^{N_{vpe}} N_i(\xi, \eta, \zeta, ...) P_i\]

where \(N_f\) is the shape functions matrix assembled for all \(N_{p}\) integration points, \(P_{points}\) is the property array at the integration points, and \(P_{nodes}\) is the property array at the mesh nodes.

To project back to the nodes, we solve the following system using the pseudo-inverse of the shape functions matrix:

\[P_{nodes} = (N_f^T W N_f)^{-1} N_f^T W P_{points}\]

where \(W\) is a diagonal matrix of weights associated with each integration point.

Demonstration#

Lets consider the following notations:

  • \(N_v\) the number of nodes in the mesh.

  • \(N_{p}\) the number of integration points.

  • \(P\) the number of property components (e.g., 1 for scalar properties, 3 for vector properties).

  • \(P_{points}\) the property array at the integration points with shape \((N_{p}, P)\).

  • \(P_{nodes}\) the property array at the mesh nodes with shape \((N_v, P)\).

  • \(N_f\) the shape functions matrix with shape \((N_{p}, N_v)\).

  • \(W\) the diagonal matrix of weights with shape \((N_{p}, N_{p})\).

We know that:

\[P_{points} = N_f \cdot P_{nodes}\]

To project back to the vertices of the mesh, we want tot minimize the following weighted least squares problem:

\[\hat{P_{nodes}} = \min_{P_{nodes}} ( \frac{1}{2}\sum_{i=1}^{N_p} w_i \| P_{points,i} - N_{f,i} \cdot P_{nodes} \|^2 )\]
\[\hat{P_{nodes}} = \min_{P_{nodes}} ( \frac{1}{2} (P_{points} - N_f \cdot P_{nodes})^T W (P_{points} - N_f \cdot P_{nodes}) )\]

The gradient of the objective function is:

\[\nabla_{P_{nodes}} \hat{P_{nodes}} = - N_f^T W (P_{points} - N_f \cdot P_{nodes})\]

To find the optimal solution, we set the gradient to zero and solve for \(P_{nodes}\):

\[\nabla_{P_{nodes}} \hat{P_{nodes}} = 0 \Rightarrow N_f^T W P_{points} = N_f^T W N_f \cdot P_{nodes}\]

So if we denote \(A = N_f^T W N_f\) invertible, we have:

\[\hat{P_{nodes}} = A^{-1} N_f^T W P_{points}\]

See also

pysdic.compute_shape_functions

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

pysdic.compute_property_projection

To project properties defined at integration points back to the nodes of the mesh using mesh and integration point information.

pysdic.assemble_property_interpolation

To interpolate properties defined at the nodes of the mesh to integration points within elements.

Examples

Lets construct a simple 2D mesh and project a scalar property defined at given integration points back to the nodes for triangular elements.

 1import numpy
 2from pysdic import (
 3    compute_shape_function_matrix,
 4    assemble_property_projection,
 5    compute_property_interpolation
 6)
 7
 8vertices_coordinates = numpy.array(
 9    [[0.0, 0.0],
10     [1.0, 0.0],
11     [1.0, 1.0],
12     [0.0, 1.0]]
13)
14
15connectivity = numpy.array(
16    [[0, 1, 2],
17     [0, 2, 3]]
18)
19
20N_e = connectivity.shape[0]
21
22property_array = numpy.array([10.0, 20.0, 30.0,  40.0])  # Scalar property
23
24# Define natural coordinates for interpolation (e.g., 4 points per element)
25natural_coordinates = numpy.array(
26    [[0.3, 0.3],
27     [0.2, 0.5],
28     [0.5, 0.2],
29     [0.1, 0.1]]
30)
31N_p = natural_coordinates.shape[0]
32natural_coordinates = numpy.vstack([natural_coordinates] * N_e)
33element_indices = numpy.repeat(numpy.arange(N_e), N_p)
34N_p = natural_coordinates.shape[0]
35
36interpolate_property = compute_property_interpolation(
37    property_array=property_array,
38    connectivity=connectivity,
39    element_type='triangle_3',
40    natural_coordinates=natural_coordinates,
41    element_indices=element_indices
42)
43
44# Now project back to the nodes
45shape_function_matrix = compute_shape_function_matrix(
46    connectivity=connectivity,
47    element_type='triangle_3',
48    natural_coordinates=natural_coordinates,
49    element_indices=element_indices,
50    n_vertices=vertices_coordinates.shape[0],
51    sparse=False
52)
53
54projected_props = assemble_property_projection(
55    property_array=property_array,
56    shape_function_matrix=shape_function_matrix
57)
58
59print(f"projected properties (shape={projected_props.shape}):")
60print(projected_props)

projected properties (shape=(4, 1)):
[[10.]
 [20.]
 [30.]
 [40.]]