Degree of Freedom (DOF) Manager and Numerics

The edelweissfe.numerics package provides core data structures and managers for global equation systems in EdelweissFE.

DofManager

This module contains important classes for describing the global equation system by means of a sparse system.

class edelweissfe.numerics.dofmanager.DofManager(nodeFields, scalarVariables=None, elements=None, constraints=None, nodeSets=None, initializeVIJPattern=True, initializeAccumulatedNodalFluxesFieldwise=True, determiningIndexToHostObjectMapping=True, **kwargs)

Bases: object

The DofManager coordinates the global equation system structure.

Specifically, DofManager:

• Analyzes the simulation domain (nodes, scalar variables, elements, constraints, and node sets).

• Collects and assigns global Degree-of-Freedom (DOF) indices for node fields and scalar variables.

• Manages entity-to-index mappings for fast scatter/gather operations.

• Generates sparse system matrix sparsity patterns (I and J vectors in COO/VIJ format).

• Calculates accumulated elemental/constraint fluxes fieldwise for convergence tests.

• Constructs domain-aware DofVector, ScatterDofVector, and VIJSystemMatrix instances.

Parameters:

• nodeFields (list of NodeField) – The list of NodeFields to be represented in the DofVector structure.

• scalarVariables (list of ScalarVariable, optional) – The list of ScalarVariables to be represented in the DofVector structure.

• elements (list of BaseNodeCouplingEntity, optional) – The list of Elements for which DOF index mappings will be created.

• constraints (list of BaseNodeCouplingEntity, optional) – The list of Constraints for which DOF index mappings will be created.

• nodeSets (list of NodeSet, optional) – The list of NodeSets for which DOF index mappings will be created.

• initializeVIJPattern (bool, optional) – Whether to initialize the sparse VIJ pattern (I and J vectors) during construction. Default is True.

• initializeAccumulatedNodalFluxesFieldwise (bool, optional) – Whether to compute accumulated nodal fluxes fieldwise during construction. Default is True.

• determiningIndexToHostObjectMapping (bool, optional) – Whether to determine the mapping from global DOF indices to their host objects (e.g. Node) during construction. Default is True.

Attributes:

dofNumberingDigest

A digest, a hexadecimal string, of what decides the index of every degree of freedom: for every field, its name, its dimension, its slice in the DofVector and the labels of its nodes in the order they are numbered; and the names of the scalar variables, which follow, in the order they are numbered. Two equation systems with the same digest give every degree of freedom the same index, so a DofVector of one is a DofVector of the other – e.g., before and after a rebuild caused by a contact search, which changes the coupling of the degrees of freedom but not their numbering. The scalar variables are named, not counted: they are those of the constraints switched on, which a constraint may change within a step (switched from outside, or by itself), so that two systems with as many scalar variables may number different ones.

The digest names the nodes by their labels, so it is the same in every run of a model: a restart checkpoint holds it together with the vectors it numbers. Computed at the first request; the numbering of a DofManager does not change.

Methods

constructDofVector()	Construct a vector with size=nDof and which has knowledge about the location of each entity.
constructVIJSystemMatrix()	Construct a VIJ (COO) Sparse System matrix object, which also has knowledge about the location of each entity.
getNodeForIndexInDofVector(index)	Find the node for a given index in the equuation system.
locateFieldOnNodeSet(field, nodeSet)	Where the degrees of freedom of a field on the nodes of a node set lie in the DofVector, e.g. those a Dirichlet boundary condition prescribes, or those a node force acts on.
refreshConstraintIndices(constraints)	Re-locate the constraints' degrees of freedom after their connectivity changed.
writeDofVectorToNodeField(dofVector, ...[, ...])	Write the values of a NodeField from the respective locations in a given DofVector; the counterpart of writeNodeFieldToDofVector().
writeDofVectorToScalarVariables(dofVector)	Set the values of the scalar variables of this equation system from their locations in a given DofVector; the counterpart of writeScalarVariablesToDofVector().
writeNodeFieldToDofVector(dofVector, ...[, ...])	Write the current values of an entire NodeField to the respective locations in a given DofVector.
writeScalarVariablesToDofVector(dofVector)	Write the values of the scalar variables of this equation system to their locations in a given DofVector.

constructDofVector()

Construct a vector with size=nDof and which has knowledge about the location of each entity.

Returns:

A DofVector.

Return type:

DofVector

constructVIJSystemMatrix()

Construct a VIJ (COO) Sparse System matrix object, which also has knowledge about the location of each entity.

Returns:

The system Matrix.

Return type:

VIJSystemMatrix

getNodeForIndexInDofVector(index)

Find the node for a given index in the equuation system.

Parameters:

index (int) – The index in the DofVector.

Returns:

The attached Node.

Return type:

Node

locateFieldOnNodeSet(field, nodeSet)

Where the degrees of freedom of a field on the nodes of a node set lie in the DofVector, e.g. those a Dirichlet boundary condition prescribes, or those a node force acts on.

Only the nodes whose field is part of this equation system count. In a finite element model, these are all nodes of the set that carry the field. A model whose equation system covers a part of its nodes only – the active domain of a material point method model, the nodes in the support of an RKPM model’s particles – leaves the other nodes of a set out. Values given per node of the set (the prescribed values of a boundary condition, the forces of a node force) are then selected by the positions of the remaining nodes in the set, which this method returns along with the indices.

The node sets given at construction are located up front, any other one at its first request.

Parameters:

• field (str) – The name of the field, one of this equation system.

• nodeSet (NodeSet) – The node set.

Returns:

• The indices of the degrees of freedom in the DofVector, node after node in the order of the set, on each node the components of the field.

• The positions of these nodes in the set: the rows of an array that holds one row per node of the set.

Return type:

tuple[np.ndarray, np.ndarray]

refreshConstraintIndices(constraints)

Re-locate the constraints’ degrees of freedom after their connectivity changed.

A contact search re-assigns which nodes a constraint couples. That changes the constraints’ DOF footprints and nothing else: no node, field, scalar variable or element is touched, so the DOF numbering and every element’s indices stay valid, and only the constraint-derived bookkeeping has to be recomputed – by the same methods, in the same order, as the constructor computed it. That is a small fraction of what constructing a DofManager for the same model costs, and it is the whole reason this method exists rather than a second constructor call.

This does NOT re-derive the degree-of-freedom layout itself: the node fields and the scalar variables keep the indices they were given. A change that moves those – a mesh refinement, a constraint gaining a scalar variable of its own – is not a connectivity change and needs a new DofManager.

The merged entity mapping is rebuilt as a new dict object on purpose: consumers cache plans keyed on its identity (the scatter template, the element loop’s gather plan), and a mapping mutated in place would leave those plans silently stale.

Parameters:

constraints (list) – The constraints, with their current connectivity.

writeDofVectorToNodeField(dofVector, nodeField, resultName, nodeSet=None)

Write the values of a NodeField from the respective locations in a given DofVector; the counterpart of writeNodeFieldToDofVector().

Parameters:

• dofVector (DofVector) – The source DofVector.

• nodeField (NodeField) – The NodeField to get the updated values. Written as a whole, it must be the one this DofManager numbers (see numberedNodeFields); written on a node set, any node field of a field this DofManager numbers, since its rows are then found by node.

• resultName (str) – The entry of the NodeField to write, one of NodeField.ENTRIES.

• nodeSet (NodeSet) – The NodeSet to write: its nodes whose field is part of this equation system (see locateFieldOnNodeSet()); the other rows of the NodeField keep their values. If None, all nodes of the NodeField are written.

Returns:

The updated NodeField.

Return type:

NodeField

Raises:

• KeyError – If the entry is not one a node field holds.

• ValueError – If this DofManager does not number the field, or, written as a whole, numbers another NodeField of that name.

writeDofVectorToScalarVariables(dofVector)

Set the values of the scalar variables of this equation system from their locations in a given DofVector; the counterpart of writeScalarVariablesToDofVector(). A scalar variable this DofManager does not number keeps its value.

Parameters:

dofVector (DofVector) – The source DofVector.

writeNodeFieldToDofVector(dofVector, nodeField, resultName, nodeSet=None)

Write the current values of an entire NodeField to the respective locations in a given DofVector.

Parameters:

• dofVector (DofVector) – The result DofVector.

• nodeField (NodeField) – The NodeField holding the values. Read as a whole, it must be the one this DofManager numbers (see numberedNodeFields); read on a node set, any node field of a field this DofManager numbers, since its rows are then found by node.

• resultName (str) – The entry of the NodeField to read, one of NodeField.ENTRIES.

• NodeSet – The NodeSet to consider: its nodes whose field is part of this equation system (see locateFieldOnNodeSet()). If None, all nodes of the NodeField are considered.

• nodeSet (NodeSet)

Returns:

The DofVector.

Return type:

DofVector

Raises:

• KeyError – If the entry is not one a node field holds.

• ValueError – If this DofManager does not number the field, or, read as a whole, numbers another NodeField of that name.

writeScalarVariablesToDofVector(dofVector)

Write the values of the scalar variables of this equation system to their locations in a given DofVector. Only the scalar variables this DofManager numbers are written: those of a constraint that is switched off are not part of the equation system.

Parameters:

dofVector (DofVector) – The result DofVector.

Returns:

The DofVector.

Return type:

DofVector

accumulatedConstraintNDof

Accumulated number of constraint DOFs (= sum of constraint vector sizes)

accumulatedElementNDof

Accumulated number of element DOFs (= sum of element vector sizes)

property dofNumberingDigest: str

for every field, its name, its dimension, its slice in the DofVector and the labels of its nodes in the order they are numbered; and the names of the scalar variables, which follow, in the order they are numbered. Two equation systems with the same digest give every degree of freedom the same index, so a DofVector of one is a DofVector of the other – e.g., before and after a rebuild caused by a contact search, which changes the coupling of the degrees of freedom but not their numbering. The scalar variables are named, not counted: they are those of the constraints switched on, which a constraint may change within a step (switched from outside, or by itself), so that two systems with as many scalar variables may number different ones.

The digest names the nodes by their labels, so it is the same in every run of a model: a restart checkpoint holds it together with the vectors it numbers. Computed at the first request; the numbering of a DofManager does not change.

Type:

A digest, a hexadecimal string, of what decides the index of every degree of freedom

fields

The list of field names in the DofVector

idcsOfConstraintsInDofVector

Dictionary mapping a constraint to its indices in the DofVector

idcsOfElementsInDofVector

Dictionary mapping an element to its indices in the DofVector

idcsOfFieldVariablesInDofVector

Dictionary mapping nodal field variables to their indices in the DofVector

idcsOfFieldsInDofVector

Dictionary mapping field names to all indices in the DofVector

idcsOfNodeFieldsInDofVector

Dictionary mapping a complete NodeField to all its indices in the DofVector

idcsOfScalarVariablesInDofVector

Dictionary mapping scalar variables to all indices in the DofVector

indexToHostObjectMapping

Reverse dictionary mapping an index to its Host object (e.g., a Node).

Warning

Not read within EdelweissFE (kept for external callers); building it costs a pass over every degree of freedom.

largestNumberOfConstraintNDof

Size of the largest constraint DOF vector

largestNumberOfElNDof

Size of the largest element DOF vector

nDof

The total number of degrees of freedom (and size of the DofVector)

DofVector

class edelweissfe.numerics.dofvector.DofVector(nDof, entitiesInDofVector=None)

Bases: ndarray

Represents a Degree-of-Freedom (DOF) vector with entity-aware indexing.

A DofVector is a 1D NumPy array subclass augmented with metadata mapping simulation entities (such as elements, constraints, or node sets) to their respective DOF indices. This enables direct indexing and assignment using entity objects as keys.

Parameters:

• nDof (int) – The total number of degrees of freedom (length of the vector).

• entitiesInDofVector (dict, optional) – A dictionary mapping entities (or entity keys) to their index arrays or slices in the DofVector.

Attributes:

T

View of the transposed array.

base

Base object if memory is from some other object.

ctypes

An object to simplify the interaction of the array with the ctypes module.

data

Python buffer object pointing to the start of the array’s data.

device

dtype

Data-type of the array’s elements.

flags

Information about the memory layout of the array.

flat

A 1-D iterator over the array.

imag

The imaginary part of the array.

itemsize

Length of one array element in bytes.

mT

View of the matrix transposed array.

nbytes

Total bytes consumed by the elements of the array.

ndim

Number of array dimensions.

real

The real part of the array.

shape

Tuple of array dimensions.

size

Number of elements in the array.

strides

Tuple of bytes to step in each dimension when traversing an array.

Methods

all([axis, out, keepdims, where])	Returns True if all elements evaluate to True.
any([axis, out, keepdims, where])	Returns True if any of the elements of a evaluate to True.
argmax([axis, out, keepdims])	Return indices of the maximum values along the given axis.
argmin([axis, out, keepdims])	Return indices of the minimum values along the given axis.
argpartition(kth[, axis, kind, order])	Returns the indices that would partition this array.
argsort([axis, kind, order, stable, descending])	Returns the indices that would sort this array.
asPlainArray()	The same buffer seen as a plain numpy.ndarray, built once and cached.
astype(dtype[, order, casting, subok, copy])	Copy of the array, cast to a specified type.
byteswap([inplace])	Swap the bytes of the array elements
choose(choices[, out, mode])	Use an index array to construct a new array from a set of choices.
clip([min, max, out])	Return an array whose values are limited to [min, max].
compress(condition[, axis, out])	Return selected slices of this array along given axis.
conj()	Complex-conjugate all elements.
conjugate()	Return the complex conjugate, element-wise.
copy([order])	Create a copy of this DofVector.
createScatterVector()	Create a scatter vector for ALL entities registered in this DofVector.
cumprod([axis, dtype, out])	Return the cumulative product of the elements along the given axis.
cumsum([axis, dtype, out])	Return the cumulative sum of the elements along the given axis.
diagonal([offset, axis1, axis2])	Return specified diagonals.
dot(other, /[, out])	Refer to numpy.dot() for full documentation.
dump(file)	Dump a pickle of the array to the specified file.
dumps()	Returns the pickle of the array as a string.
fill(value)	Fill the array with a scalar value.
flatten([order])	Return a copy of the array collapsed into one dimension.
getfield(dtype[, offset])	Returns a field of the given array as a certain type.
item(*args)	Copy an element of an array to a standard Python scalar and return it.
max([axis, out, keepdims, initial, where])	Return the maximum along a given axis.
mean([axis, dtype, out, keepdims, where])	Returns the average of the array elements along given axis.
min([axis, out, keepdims, initial, where])	Return the minimum along a given axis.
nonzero()	Return the indices of the elements that are non-zero.
partition(kth[, axis, kind, order])	Partially sorts the elements in the array in such a way that the value of the element in k-th position is in the position it would be in a sorted array.
prod([axis, dtype, out, keepdims, initial, ...])	Return the product of the array elements over the given axis
put(indices, values[, mode])	Set a.flat[n] = values[n] for all n in indices.
ravel([order])	Return a flattened array.
repeat(repeats[, axis])	Repeat elements of an array.
reshape(a.reshape)	Returns an array containing the same data with a new shape.
resize(a.resize)	Change shape and size of array in-place.
round([decimals, out])	Return a with each element rounded to the given number of decimals.
searchsorted(v[, side, sorter])	Find indices where elements of v should be inserted in a to maintain order.
setfield(val, dtype[, offset])	Put a value into a specified place in a field defined by a data-type.
setflags([write, align, uic])	Set array flags WRITEABLE, ALIGNED, WRITEBACKIFCOPY, respectively.
sort([axis, kind, order, stable, descending])	Sort an array in-place.
squeeze([axis])	Remove axes of length one from a.
std([axis, dtype, out, ddof, keepdims, ...])	Returns the standard deviation of the array elements along given axis.
sum([axis, dtype, out, keepdims, initial, where])	Return the sum of the array elements over the given axis.
swapaxes(axis1, axis2, /)	Return a view of the array with axis1 and axis2 interchanged.
take(indices[, axis, out, mode])	Return an array formed from the elements of a at the given indices.
to_device(device, /, *[, stream])	For Array API compatibility.
tobytes([order])	Construct Python bytes containing the raw data bytes in the array.
tofile(fid, /[, sep, format])	Write array to a file as text or binary (default).
tolist()	Return the array as an a.ndim-levels deep nested list of Python scalars.
trace([offset, axis1, axis2, dtype, out])	Return the sum along diagonals of the array.
transpose(*axes)	Returns a view of the array with axes transposed.
var([axis, dtype, out, ddof, keepdims, ...])	Returns the variance of the array elements, along given axis.
view([dtype][, type])	New view of array with the same data.

asPlainArray()

The same buffer seen as a plain numpy.ndarray, built once and cached.

Indexing through this bypasses the subclass wrapping that dominates hot-path access – see __getitem__(). It aliases this vector rather than copying it, so writes through it are writes to this vector. Callers that index a DofVector repeatedly in a loop should take this once outside the loop instead of calling view(np.ndarray) themselves.

Returns:

A plain-ndarray view of this vector’s buffer.

Return type:

np.ndarray

copy(order='C')

Create a copy of this DofVector.

Parameters:

order (str, optional) – The memory layout order (‘C’ for C-contiguous, ‘F’ for Fortran-contiguous). Default is ‘C’.

Returns:

A new DofVector instance with copied array data and an independent copy of entitiesInDofVector.

Return type:

DofVector

createScatterVector()

Create a scatter vector for ALL entities registered in this DofVector.

The underlying layout (entity lookup map and scatter indices) is computed once and cached, so repeated calls (e.g. once per Newton iteration) only allocate the zero-initialized data buffer.

Returns:

A ScatterDofVector initialized with this vector’s entity mapping and total DOF count.

Return type:

ScatterDofVector

Raises:

ValueError – If entitiesInDofVector is None or uninitialized.

ScatterDofVector

class edelweissfe.numerics.scatterdofvector.ScatterDofVector(template)

Bases: ndarray

A Scatter Vector that stores data for entities contiguously. Includes a fast lookup map to support random access by Entity.

Parameters:

template (ScatterDofVectorTemplate) – The precomputed layout of the scatter vector.

Attributes:

T

View of the transposed array.

base

Base object if memory is from some other object.

ctypes

An object to simplify the interaction of the array with the ctypes module.

data

Python buffer object pointing to the start of the array’s data.

device

dtype

Data-type of the array’s elements.

flags

Information about the memory layout of the array.

flat

A 1-D iterator over the array.

imag

The imaginary part of the array.

itemsize

Length of one array element in bytes.

mT

View of the matrix transposed array.

nbytes

Total bytes consumed by the elements of the array.

ndim

Number of array dimensions.

offsetMap

The entity-to-(offset, size) layout of this buffer.

real

The real part of the array.

shape

Tuple of array dimensions.

size

Number of elements in the array.

strides

Tuple of bytes to step in each dimension when traversing an array.

Methods

all([axis, out, keepdims, where])	Returns True if all elements evaluate to True.
any([axis, out, keepdims, where])	Returns True if any of the elements of a evaluate to True.
argmax([axis, out, keepdims])	Return indices of the maximum values along the given axis.
argmin([axis, out, keepdims])	Return indices of the minimum values along the given axis.
argpartition(kth[, axis, kind, order])	Returns the indices that would partition this array.
argsort([axis, kind, order, stable, descending])	Returns the indices that would sort this array.
assembleInto(targetDofVector[, absolute])	Scatter-add values into a global DOF vector.
astype(dtype[, order, casting, subok, copy])	Copy of the array, cast to a specified type.
byteswap([inplace])	Swap the bytes of the array elements
choose(choices[, out, mode])	Use an index array to construct a new array from a set of choices.
clip([min, max, out])	Return an array whose values are limited to [min, max].
compress(condition[, axis, out])	Return selected slices of this array along given axis.
conj()	Complex-conjugate all elements.
conjugate()	Return the complex conjugate, element-wise.
copy([order])	Create a copy of this ScatterDofVector.
cumprod([axis, dtype, out])	Return the cumulative product of the elements along the given axis.
cumsum([axis, dtype, out])	Return the cumulative sum of the elements along the given axis.
diagonal([offset, axis1, axis2])	Return specified diagonals.
dot(other, /[, out])	Refer to numpy.dot() for full documentation.
dump(file)	Dump a pickle of the array to the specified file.
dumps()	Returns the pickle of the array as a string.
fill(value)	Fill the array with a scalar value.
flatten([order])	Return a copy of the array collapsed into one dimension.
getfield(dtype[, offset])	Returns a field of the given array as a certain type.
item(*args)	Copy an element of an array to a standard Python scalar and return it.
max([axis, out, keepdims, initial, where])	Return the maximum along a given axis.
mean([axis, dtype, out, keepdims, where])	Returns the average of the array elements along given axis.
min([axis, out, keepdims, initial, where])	Return the minimum along a given axis.
nonzero()	Return the indices of the elements that are non-zero.
partition(kth[, axis, kind, order])	Partially sorts the elements in the array in such a way that the value of the element in k-th position is in the position it would be in a sorted array.
prod([axis, dtype, out, keepdims, initial, ...])	Return the product of the array elements over the given axis
put(indices, values[, mode])	Set a.flat[n] = values[n] for all n in indices.
ravel([order])	Return a flattened array.
repeat(repeats[, axis])	Repeat elements of an array.
reshape(a.reshape)	Returns an array containing the same data with a new shape.
resize(a.resize)	Change shape and size of array in-place.
round([decimals, out])	Return a with each element rounded to the given number of decimals.
searchsorted(v[, side, sorter])	Find indices where elements of v should be inserted in a to maintain order.
setfield(val, dtype[, offset])	Put a value into a specified place in a field defined by a data-type.
setflags([write, align, uic])	Set array flags WRITEABLE, ALIGNED, WRITEBACKIFCOPY, respectively.
sort([axis, kind, order, stable, descending])	Sort an array in-place.
squeeze([axis])	Remove axes of length one from a.
std([axis, dtype, out, ddof, keepdims, ...])	Returns the standard deviation of the array elements along given axis.
sum([axis, dtype, out, keepdims, initial, where])	Return the sum of the array elements over the given axis.
swapaxes(axis1, axis2, /)	Return a view of the array with axis1 and axis2 interchanged.
take(indices[, axis, out, mode])	Return an array formed from the elements of a at the given indices.
toDofVector([absolute])	Create a new global DofVector by scattering values from this scatter vector.
to_device(device, /, *[, stream])	For Array API compatibility.
tobytes([order])	Construct Python bytes containing the raw data bytes in the array.
tofile(fid, /[, sep, format])	Write array to a file as text or binary (default).
tolist()	Return the array as an a.ndim-levels deep nested list of Python scalars.
trace([offset, axis1, axis2, dtype, out])	Return the sum along diagonals of the array.
transpose(*axes)	Returns a view of the array with axes transposed.
var([axis, dtype, out, ddof, keepdims, ...])	Returns the variance of the array elements, along given axis.
view([dtype][, type])	New view of array with the same data.

assembleInto(targetDofVector, absolute=False)

Scatter-add values into a global DOF vector.

Parameters:

• targetDofVector (ndarray) – The target vector (e.g. DofVector) to assemble into.

• absolute (bool, optional) – If True, accumulate absolute values. Default is False.

Return type:

None

copy(order='C')

Create a copy of this ScatterDofVector.

Parameters:

order (str, optional) – The memory layout order (‘C’ for C-contiguous, ‘F’ for Fortran-contiguous). Default is ‘C’.

Returns:

A new ScatterDofVector instance with copied array data and metadata.

Return type:

ScatterDofVector

toDofVector(absolute=False)

Create a new global DofVector by scattering values from this scatter vector.

Parameters:

absolute (bool, optional) – If True, use absolute values. Default is False.

Returns:

A new DofVector containing the accumulated entries.

Return type:

DofVector

property offsetMap: dict

The entity-to-(offset, size) layout of this buffer.

For a caller that plans its own writes ahead of time – the explicit element loop, which assembles a whole chunk of elements into a buffer of its own and then places it here in one go – instead of taking a view per entity.

Returns:

Mapping of each entity to its (offset, size) within this buffer.

Return type:

dict

class edelweissfe.numerics.scatterdofvector.ScatterDofVectorTemplate(entitiesInDofVector, nDof)

Bases: object

The precomputed, immutable layout of a ScatterDofVector: the entity lookup map and the gather/scatter index array.

Since the layout only depends on the entities and their indices in the DofVector, it can be computed once and shared by all scatter vectors of the same equation system, e.g. across Newton iterations.

Parameters:

• entitiesInDofVector (dict) – The dictionary mapping entities to their global index arrays in the DofVector.

• nDof (int) – The total number of global degrees of freedom.

VIJSystemMatrix

class edelweissfe.numerics.vijsystemmatrix.VIJSystemMatrix(nDof, I, J, entitiesInVIJ)

Bases: ndarray

Represents the value vector V of a COO (VIJ triple) sparse matrix.

VIJSystemMatrix stores matrix entry values for sparse global systems (Coordinate format), maintaining associated row indices I, column indices J, and an entitiesInVIJ mapping. It supports entity-aware indexing to slice, reshape, and assign local entity matrices directly.

Parameters:

• nDof (int) – The dimension (number of rows/columns) of the system matrix.

• I (ndarray) – The 1D integer array of global row indices.

• J (ndarray) – The 1D integer array of global column indices.

• entitiesInVIJ (dict) – A dictionary mapping entities to their starting offset index in V.

Attributes:

T

View of the transposed array.

base

Base object if memory is from some other object.

ctypes

An object to simplify the interaction of the array with the ctypes module.

data

Python buffer object pointing to the start of the array’s data.

device

dtype

Data-type of the array’s elements.

flags

Information about the memory layout of the array.

flat

A 1-D iterator over the array.

imag

The imaginary part of the array.

itemsize

Length of one array element in bytes.

mT

View of the matrix transposed array.

nbytes

Total bytes consumed by the elements of the array.

ndim

Number of array dimensions.

real

The real part of the array.

shape

Tuple of array dimensions.

size

Number of elements in the array.

strides

Tuple of bytes to step in each dimension when traversing an array.

Methods

all([axis, out, keepdims, where])	Returns True if all elements evaluate to True.
any([axis, out, keepdims, where])	Returns True if any of the elements of a evaluate to True.
argmax([axis, out, keepdims])	Return indices of the maximum values along the given axis.
argmin([axis, out, keepdims])	Return indices of the minimum values along the given axis.
argpartition(kth[, axis, kind, order])	Returns the indices that would partition this array.
argsort([axis, kind, order, stable, descending])	Returns the indices that would sort this array.
astype(dtype[, order, casting, subok, copy])	Copy of the array, cast to a specified type.
byteswap([inplace])	Swap the bytes of the array elements
choose(choices[, out, mode])	Use an index array to construct a new array from a set of choices.
clip([min, max, out])	Return an array whose values are limited to [min, max].
compress(condition[, axis, out])	Return selected slices of this array along given axis.
conj()	Complex-conjugate all elements.
conjugate()	Return the complex conjugate, element-wise.
copy([order])	Create a copy of this VIJSystemMatrix.
cumprod([axis, dtype, out])	Return the cumulative product of the elements along the given axis.
cumsum([axis, dtype, out])	Return the cumulative sum of the elements along the given axis.
diagonal([offset, axis1, axis2])	Return specified diagonals.
dot(other, /[, out])	Refer to numpy.dot() for full documentation.
dump(file)	Dump a pickle of the array to the specified file.
dumps()	Returns the pickle of the array as a string.
fill(value)	Fill the array with a scalar value.
flatten([order])	Return a copy of the array collapsed into one dimension.
getfield(dtype[, offset])	Returns a field of the given array as a certain type.
item(*args)	Copy an element of an array to a standard Python scalar and return it.
max([axis, out, keepdims, initial, where])	Return the maximum along a given axis.
mean([axis, dtype, out, keepdims, where])	Returns the average of the array elements along given axis.
min([axis, out, keepdims, initial, where])	Return the minimum along a given axis.
nonzero()	Return the indices of the elements that are non-zero.
partition(kth[, axis, kind, order])	Partially sorts the elements in the array in such a way that the value of the element in k-th position is in the position it would be in a sorted array.
prod([axis, dtype, out, keepdims, initial, ...])	Return the product of the array elements over the given axis
put(indices, values[, mode])	Set a.flat[n] = values[n] for all n in indices.
ravel([order])	Return a flattened array.
repeat(repeats[, axis])	Repeat elements of an array.
reshape(a.reshape)	Returns an array containing the same data with a new shape.
resize(a.resize)	Change shape and size of array in-place.
round([decimals, out])	Return a with each element rounded to the given number of decimals.
searchsorted(v[, side, sorter])	Find indices where elements of v should be inserted in a to maintain order.
setfield(val, dtype[, offset])	Put a value into a specified place in a field defined by a data-type.
setflags([write, align, uic])	Set array flags WRITEABLE, ALIGNED, WRITEBACKIFCOPY, respectively.
sort([axis, kind, order, stable, descending])	Sort an array in-place.
squeeze([axis])	Remove axes of length one from a.
std([axis, dtype, out, ddof, keepdims, ...])	Returns the standard deviation of the array elements along given axis.
sum([axis, dtype, out, keepdims, initial, where])	Return the sum of the array elements over the given axis.
swapaxes(axis1, axis2, /)	Return a view of the array with axis1 and axis2 interchanged.
take(indices[, axis, out, mode])	Return an array formed from the elements of a at the given indices.
to_device(device, /, *[, stream])	For Array API compatibility.
tobytes([order])	Construct Python bytes containing the raw data bytes in the array.
tofile(fid, /[, sep, format])	Write array to a file as text or binary (default).
tolist()	Return the array as an a.ndim-levels deep nested list of Python scalars.
trace([offset, axis1, axis2, dtype, out])	Return the sum along diagonals of the array.
transpose(*axes)	Returns a view of the array with axes transposed.
var([axis, dtype, out, ddof, keepdims, ...])	Returns the variance of the array elements, along given axis.
view([dtype][, type])	New view of array with the same data.

copy(order='C')

Create a copy of this VIJSystemMatrix.

Parameters:

order (str, optional) – The memory layout order (‘C’ for C-contiguous, ‘F’ for Fortran-contiguous). Default is ‘C’.

Returns:

A new VIJSystemMatrix instance with copied values and metadata.

Return type:

VIJSystemMatrix

VIJEntityBase

class edelweissfe.numerics.vijentitybase.VIJEntityBase(nDof)

Bases: object

Initialize the VIJEntityBase.

Attributes:

couplesAllItsDofs

Whether the entity’s block couples every one of its degrees of freedom with every other one.

Parameters:

nDof (int)

Methods

getVIJContributionSize()	Return the number of entries this entity contributes to the VIJ (COO) system matrix.
initializeVIJContribution(idcs, I_, J_, offset)	Fill the global I and J index arrays for this constraint's VIJ contribution.
shapeVIJContribution(flat_view)	Shape the flat VIJ values slice for this entity.

getVIJContributionSize()

Return the number of entries this entity contributes to the VIJ (COO) system matrix.

By default this is nDof**2, which corresponds to a full dense block.

Returns:

Number of VIJ entries for this entity.

Return type:

int

initializeVIJContribution(idcs, I_, J_, offset)

Fill the global I and J index arrays for this constraint’s VIJ contribution.

The default implementation writes a full dense nDof × nDof block, identical to the behaviour used for finite elements. Entities that override getVIJContributionSize() to return a smaller value must also override this method so that exactly getVIJContributionSize() (I, J) pairs are written starting at offset.

Parameters:

• idcs (ndarray) – Global DOF indices for this entity (length = nDof).

• I – The global row-index array of the VIJ triple (written in-place).

• J – The global column-index array of the VIJ triple (written in-place).

• offset (int) – First position in I / J that belongs to this entity.

• I_ (ndarray)

• J_ (ndarray)

Return type:

None

shapeVIJContribution(flat_view)

Shape the flat VIJ values slice for this entity.

By default, if the contribution is dense (couplesAllItsDofs), it reshapes the flat view to a 2-D array of shape (nDof, nDof) using column-major (Fortran) order.

Parameters:

flat_view (ndarray)

Return type:

ndarray

property couplesAllItsDofs: bool

Whether the entity’s block couples every one of its degrees of freedom with every other one.

True (the default, every element) means the block is the dense nDof × nDof matrix, stored column-major, with the entity’s global DOF indices as its rows and its columns. False means the entity declared a sparse pattern of its own in initializeVIJContribution() (a constraint that couples only some of its degrees of freedom), and its block is the list of the values of these (I, J) pairs, in the order it declared them.

The entity says which, because the number of entries cannot tell: a sparse pattern can happen to have exactly nDof**2 entries. A rigid body contact with a single slave node in 3D couples the six degrees of freedom of the reference point with each other, the three of the node with each other, and the two groups both ways: 36 + 9 + 18 + 18 = 81 = 9², yet its entries are ordered block by block, not column by column.

Returns:

True for a dense block, False for a pattern of the entity’s own.

Return type:

bool

Parallelization Utilities

edelweissfe.numerics.parallelizationutilities.chunked_iterable(iterable, size)

Yield successive n-sized chunks from an iterable.

edelweissfe.numerics.parallelizationutilities.forEachSliceInParallel(work, length, sliceLength=1048576)

Call work(s) for consecutive slices s that together cover range(length), on the persistent thread pool, and return what the calls returned, in the order of the slices.

Meant for the simple array operations on arrays as long as the system matrix – zeroing it, adding a scaled mass matrix to it – which numpy runs on a single thread, and which on a large model take seconds per Newton iteration that way. Every call works on its own part of the arrays, so the result is bitwise the one of the same operation on the whole array. Short slices keep the temporaries of an expression like K[s] += a * M[s] small.

Runs serially in the calling thread if free threading is not available or the array fits into one slice.

When it returns or raises, no call is running. If a call raises, the exception is that of the first failing slice in the order of the slices, the one a serial loop meets first; the slices that have not started are not made, and the running ones are waited for.

Parameters:

• work – Called once per slice, with a slice of range(length); what it returns (e.g. the indices found in the slice) is collected.

• length (int) – The length of the arrays work operates on.

• sliceLength (int) – The length of the slices; the last one may be shorter.

Returns:

What work returned for each slice, in order.

Return type:

list

edelweissfe.numerics.parallelizationutilities.getNumberOfAvailableCpus()

Get the number of CPUs this process is actually permitted to run on.

This is not the machine’s core count. A batch scheduler’s cgroup, taskset, and an OpenMP runtime honouring OMP_PROC_BIND all narrow it – the last of them at the moment it is loaded, which in this framework is whenever the first compiled extension is imported.

Returns:

The number of CPUs available to this process, or 0 where the platform cannot report it.

Return type:

int

edelweissfe.numerics.parallelizationutilities.getNumberOfThreads()

Get the number of threads available for parallel computations.

EdelweissFE has a built-in mechanism to determine the number of threads to be used for parallel computations. It checks if the environment variable OMP_NUM_THREADS is set. If it is and can be converted to an integer, that value is used. If not, the function falls back to using a default value of 1. The result is clamped to at least 1, since a ThreadPoolExecutor requires at least one worker.

Returns:

Number of threads to be used.

Return type:

int

edelweissfe.numerics.parallelizationutilities.getNumberOfWorkerThreads()

The number of threads a parallel computation on the thread pool runs with: getNumberOfThreads() if free threading is available, else 1 – without free threading the workers would only take turns on the GIL, so the work runs serially in the calling thread.

Returns:

The number of worker threads, at least 1.

Return type:

int

edelweissfe.numerics.parallelizationutilities.getThreadPool(numThreads)

Get a persistent thread pool with the requested number of worker threads.

Pools are created lazily and reused for the lifetime of the process. This avoids the cost of spawning and joining worker threads for every parallel computation, which the solvers request once per Newton iteration.

Parameters:

numThreads (int) – The number of worker threads. Clamped to at least 1.

Returns:

The persistent thread pool.

Return type:

concurrent.futures.ThreadPoolExecutor

edelweissfe.numerics.parallelizationutilities.isFreeThreadingSupported()

Check if free threading is supported in the current build of EdelweissFE.

Free threading allows for parallel computations using multiple threads. This function checks if the current build of EdelweissFE supports this feature.

Returns:

True if free threading is supported and the GIL is disabled, False otherwise.

Return type:

bool

edelweissfe.numerics.parallelizationutilities.mapChunksInParallel(items, work, chunksPerThread=4)

Call work(chunk) for consecutive chunks of items on the persistent thread pool, and return what the calls returned, in the order of the chunks.

Meant for a loop over entities (the elements) whose iterations are independent: each call handles a whole chunk, so that the cost of handing out tasks stays negligible against the work, and there are several chunks per thread, so that the pool balances chunks of unequal cost. The results come in the order of items: what is assembled from them is what the serial loop assembles.

Runs serially in the calling thread, as one chunk holding all items, if free threading is not available or one thread is asked for.

When it returns or raises, no call is running. If a call raises, the exception is that of the first failing chunk in the order of the chunks, the one a serial loop meets first; the chunks that have not started are not made, and the running ones are waited for.

Parameters:

• items – The entities, a sized iterable (a list, a dictionary view).

• work – Called once per chunk, with a tuple of consecutive items; a worker must write only to what belongs to its items.

• chunksPerThread (int) – How many chunks each thread is handed.

Returns:

What work returned for each chunk, in order.

Return type:

list

edelweissfe.numerics.parallelizationutilities.mapChunksOnWorkers(items, work, numberOfWorkers, chunksPerWorker=4)

Call work(worker, chunk) for consecutive chunks of items, each chunk on a fixed worker thread, and return what the calls returned, in the order of the chunks.

For an entity loop whose workers write through memory of their own, which no other call may use at the same time – each worker its own scratch array for the blocks of the system matrix (see blocksOfWorker()). Unlike mapChunksInParallel(), which hands each chunk to whichever thread is free, the chunks are dealt out to the workers in turn: worker w handles chunks w, w + numberOfWorkers, …, in this order, so every worker meets the same entities in every run. Several chunks per worker, interleaved, spread a region of expensive entities (e.g., at a damage front, adjacent in the element order) over all workers.

Runs serially in the calling thread, as worker 0 with one chunk holding all items, for one worker.

When it returns or raises, no call is running. If a call raises, the exception is that of the first failing chunk in the order of the chunks, the one a serial loop meets first: a worker stops at its first failing chunk, so every chunk before the first failing one has run, and none of them failed.

Parameters:

• items – The entities, a sized iterable (a list, a dictionary view).

• work – Called once per chunk, with the number of the worker (0 <= worker < numberOfWorkers) and a tuple of consecutive items; a worker must write only to what belongs to its items, or to what belongs to the worker.

• numberOfWorkers (int) – The number of worker threads.

• chunksPerWorker (int) – How many chunks each worker is handed.

Returns:

What work returned for each chunk, in order.

Return type:

list

edelweissfe.numerics.parallelizationutilities.reportThreadAvailability(numThreads, journal, senderIdentification)

Report how many threads will be used, and warn if the process cannot actually use them.

A pool of numThreads workers only runs on numThreads cores if the process is permitted to use that many, and OMP_PROC_BIND/OMP_PLACES are the common reason for it not to be: they make the OpenMP runtime pin the thread that loads it to a single place, and every thread started afterwards – the pool’s workers included – inherits that one-core mask. The element loop is then serial in fact while still reporting the threads it asked for, and nothing fails, so this is worth saying out loud rather than leaving to be discovered by measurement.

Parameters:

• numThreads (int) – The number of threads the solver is about to use.

• journal – The journal to report to.

• senderIdentification (str) – The name of the reporting solver.