# SPDX-License-Identifier: MIT
# Copyright (c) 2011–2026 Joris J.C. Remmers
import numpy as np
from numpy import zeros, ones, ix_, append, repeat, array
from scipy.sparse import coo_matrix
from typing import Any, Tuple, Iterable, Optional
from numpy.typing import NDArray
from pyfem.util.dataStructures import Properties
from pyfem.util.dataStructures import elementData
#-------------------------------------------------------------------------------
# Matrix builder
#-------------------------------------------------------------------------------
[docs]
class MatrixBuilder:
"""
Incrementally build sparse matrix data in COO format for finite element assembly.
Also stores a global vector B and a scalar c for additional model data.
Attributes:
nDofs (int): Number of global degrees of freedom.
val (np.ndarray): Values of the matrix entries.
row (np.ndarray): Row indices for COO format.
col (np.ndarray): Column indices for COO format.
B (np.ndarray): Global vector (initialized to zeros).
c (float): Scalar value (initialized to 0.0).
"""
def __init__(self, nDofs: int) -> None:
"""
Initialize a new builder for a system with nDofs.
Args:
nDofs (int): Number of global degrees of freedom.
"""
self.nDofs: int = nDofs
self.clear()
[docs]
def clear(self) -> None:
"""
Reset stored row/col/value arrays, B, and c.
"""
self.val: NDArray[np.floating] = array([], dtype=float)
self.row: NDArray[np.integer] = array([], dtype=int)
self.col: NDArray[np.integer] = array([], dtype=int)
self.B: NDArray[np.floating] = zeros(self.nDofs, dtype=float)
self.c: float = 0.0
[docs]
def append(self, a: NDArray[np.floating], dofs: NDArray[np.integer]) -> None:
"""
Append element matrix `a` using associated dof indices.
Args:
a (np.ndarray): Element matrix (square, shape [n, n]).
dofs (np.ndarray): DOF indices for the element (length n).
"""
n = len(dofs)
self.row = append(self.row, repeat(dofs, n))
self.col = append(self.col, np.tile(dofs, n))
self.val = append(self.val, a.reshape(n * n))
[docs]
def getMatrix(self) -> coo_matrix:
"""
Return the assembled COO sparse matrix.
Returns:
scipy.sparse.coo_matrix: Assembled sparse matrix.
"""
return coo_matrix((self.val, (self.row, self.col)), shape=(self.nDofs, self.nDofs))
#-------------------------------------------------------------------------------
# Prepare
#-------------------------------------------------------------------------------
[docs]
def prepare(props: Properties, globdat: Any) -> None:
"""
Commit element states by calling the element 'commit' method.
This function is called after a successful time step or load step to
finalize and store the current element states (e.g., history variables,
plastic strains, damage parameters).
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
"""
globdat.models.takeAction( "prepare" , props , globdat )
return None
#-------------------------------------------------------------------------------
# Assemble Internal force
#-------------------------------------------------------------------------------
[docs]
def assembleInternalForce(props: Properties, globdat: Any) -> NDArray[np.floating]:
"""
Assemble and return the global internal force vector.
Computes the internal force vector by calling the 'getInternalForce'
method on all elements and assembling their contributions.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Returns:
np.ndarray: The assembled internal force vector.
"""
globdat.mbuilder = MatrixBuilder(len(globdat.dofs))
globdat.resetNodalOutput()
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "getInternalForce"):
element.getInternalForce(elemdat)
globdat.mbuilder.B[elemdat.el_dofs] += elemdat.fint
globdat.models.takeAction( "getInternalForce" , props , globdat )
return globdat.mbuilder.B
#-------------------------------------------------------------------------------
# Assemble Internal force
#-------------------------------------------------------------------------------
[docs]
def assembleExternalForce(props: Properties, globdat: Any) -> NDArray[np.floating]:
"""
Assemble and return the global external force vector.
Computes the external force vector by calling the 'getExternalForce'
method on all elements and assembling their contributions. The external
force returned includes contributions assembled from elements plus any
scaled forcing term stored on ``globdat``.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Returns:
np.ndarray: The assembled external force vector, including
the scaled load factor contribution (globdat.fhat * globdat.solverStatus.lam).
"""
globdat.mbuilder = MatrixBuilder(len(globdat.dofs))
globdat.resetNodalOutput()
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "getExternalForce"):
element.getExternalForce(elemdat)
globdat.mbuilder.B[elemdat.el_dofs] += elemdat.fint
globdat.models.takeAction( "getExternalForce" , props , globdat )
return globdat.mbuilder.B + globdat.fhat * globdat.solverStatus.lam
#-------------------------------------------------------------------------------
# Assemble Dissipation
#-------------------------------------------------------------------------------
[docs]
def assembleDissipation(props: Properties, globdat: Any) -> Tuple[NDArray[np.floating], float]:
"""
Assemble and return dissipation contributions.
Computes dissipation by calling the 'getDissipation' method on all elements.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Returns:
tuple[np.ndarray, float]:
- dissipation_vector: Assembled dissipation force vector
- accumulated_dissipation: Total scalar dissipation from all elements
"""
globdat.mbuilder = MatrixBuilder(len(globdat.dofs))
globdat.resetNodalOutput()
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "getDissipation"):
element.getDissipation(elemdat)
globdat.mbuilder.B[elemdat.el_dofs] += elemdat.fint
globdat.mbuilder.c += elemdat.diss
globdat.models.takeAction( "getDissipation" , props , globdat )
return globdat.mbuilder.B, globdat.mbuilder.c
#-------------------------------------------------------------------------------
# Assemble Tangent stiffness
#-------------------------------------------------------------------------------
[docs]
def assembleTangentStiffness(props: Properties, globdat: Any) -> Tuple[coo_matrix, NDArray[np.floating]]:
"""
Assemble and return the global tangent stiffness matrix and residual.
Computes the tangent stiffness matrix by calling the 'getTangentStiffness'
method on all elements and assembling their contributions into a sparse matrix.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Returns:
tuple[scipy.sparse.coo_matrix, np.ndarray]:
- stiff_matrix: Global tangent stiffness matrix in COO sparse format
- residual_vector: Assembled internal force residual vector
"""
globdat.mbuilder = MatrixBuilder(len(globdat.dofs))
globdat.resetNodalOutput()
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "getTangentStiffness"):
element.getTangentStiffness(elemdat)
globdat.mbuilder.append(elemdat.stiff, elemdat.el_dofs)
globdat.mbuilder.B[elemdat.el_dofs] += elemdat.fint
globdat.models.takeAction( "getTangentStiffness" , props , globdat )
return globdat.mbuilder.getMatrix(), globdat.mbuilder.B
#-------------------------------------------------------------------------------
# Assemble Mass Matrix
#-------------------------------------------------------------------------------
[docs]
def assembleMassMatrix(props: Properties, globdat: Any) -> Tuple[coo_matrix, NDArray[np.floating]]:
"""
Assemble and return the global mass matrix and lumped mass vector.
Computes the mass matrix by calling the 'getMassMatrix' method on all
elements and assembling their contributions into a sparse matrix.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Returns:
tuple[scipy.sparse.coo_matrix, np.ndarray]:
- mass_matrix: Global mass matrix in COO sparse format
- lumped_mass_vector: Assembled lumped mass vector (diagonal approximation)
"""
globdat.mbuilder = MatrixBuilder(len(globdat.dofs))
globdat.resetNodalOutput()
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "getMassMatrix"):
element.getMassMatrix(elemdat)
globdat.mbuilder.append(elemdat.mass, elemdat.el_dofs)
globdat.mbuilder.B[elemdat.el_dofs] += elemdat.lumped
globdat.models.takeAction( "getMassMatrix" , props , globdat )
return globdat.mbuilder.getMatrix(), globdat.mbuilder.B
#-------------------------------------------------------------------------------
# Commit
#-------------------------------------------------------------------------------
[docs]
def commit(props: Properties, globdat: Any) -> None:
"""
Commit element states by calling the element 'commit' method.
This function is called after a successful time step or load step to
finalize and store the current element states (e.g., history variables,
plastic strains, damage parameters).
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
"""
for elementGroup in globdat.elements.iterGroupNames():
el_props = getattr(props, elementGroup)
for iElm, element in enumerate(globdat.elements.iterElementGroup(elementGroup)):
elemdat = getElementData(iElm, element, el_props, globdat)
if hasattr(element, "mat"):
element.mat.reset()
if hasattr(element, "commit"):
element.commit(elemdat)
globdat.models.takeAction( "commit" , props , globdat )
return None
#-------------------------------------------------------------------------------
# getAllConstraints
#-------------------------------------------------------------------------------
[docs]
def getAllConstraints(props: Properties, globdat: Any) -> None:
"""
Invoke 'getConstraints' on all elements to collect constraint data.
This function iterates over all element groups and elements, calling their
'getConstraints' method if available. This is typically used for multi-point
constraints, contact constraints, or other element-level constraint definitions.
Args:
props (Properties): Global properties container.
globdat (Any): Global data/state object.
Note:
The current implementation creates a minimal elemdat structure for each
element. This mirrors the original behavior and intentionally does not
change the logic.
"""
# Loop over all element groups
for elementGroup in globdat.elements.iterGroupNames():
# Get the properties corresponding to the element group
el_props = getattr(props, elementGroup)
# Pre-create a placeholder element data container so that
# calls to element.getConstraints can always receive an object
elemdat = elementData(np.array([]), np.array([]))
# Loop over all elements in the element group
for element in globdat.elements.iterElementGroup(elementGroup):
# Get the element node indices
el_nodes = element.getNodes()
# Populate element data with nodes and properties
elemdat.nodes = el_nodes
elemdat.props = el_props
# Call the getConstraints method if it exists on the element
getattr(element, "getConstraints", None)(elemdat)
#-------------------------------------------------------------------------------
# getElementData
#-------------------------------------------------------------------------------
[docs]
def getElementData(iElm: int, element: Any, el_props: Properties, globdat: Any) -> elementData:
"""
Create and populate an elementData instance for an element.
This helper function gathers all necessary data for an element from the
global data structure, including:
- Node indices and coordinates
- Degree of freedom indices and values
- Current state vector and state increment
- Element properties and material properties
Args:
iElm (int): Element index in the group.
element (Any): The element object for which to gather data.
el_props (Properties): Properties object for the element's group.
globdat (Any): Global data/state object containing mesh, DOFs, and state.
Returns:
elementData: An instance populated with all element-specific information
needed for element computations.
"""
# Get element node indices
el_nodes = element.getNodes()
# Get element node coordinates from global node array
el_coords = globdat.nodes.getNodeCoords(el_nodes)
# Get element DOF indices based on node indices and DOF types
el_dofs = globdat.dofs.getForTypes(el_nodes, element.dofTypes)
# Extract current state and state increment for element DOFs
el_a = globdat.state[el_dofs]
el_Da = globdat.Dstate[el_dofs]
# Create elementData object with state vectors
elemdat = elementData(el_a, el_Da)
# Populate elementData with additional information
elemdat.coords = el_coords
elemdat.nodes = el_nodes
elemdat.props = el_props
elemdat.el_dofs = el_dofs
elemdat.iElm = iElm
# Attach global data to element for access if needed
element.globdat = globdat
# Add material properties if element has them
if hasattr(element, "matProps"):
elemdat.matprops = element.matProps
return elemdat