Overview and workflow

pyfe3d gives finite elements, not a pre-processor: the model is built in Python, with NumPy arrays for the nodes and the boundary conditions, and each element writes its contributions into sparse arrays that SciPy turns into the structural matrices. The same arrays serve linear static, linear buckling, natural frequency and geometrically non-linear analyses, solved with pyfe3d.solver or with any SciPy solver.

The pieces of a model

Object	Role
nodal coordinates	1-D array x of size 3*M, \(X_1, Y_1, Z_1, X_2, ...\), see Conventions
property	ShellProp for the shells, created by isotropic_plate() or laminated_plate(), see Shell and beam properties; BeamProp for the beams and the truss; the stiffnesses of a Spring are attributes of the element
*Data	e.g. Quad4Data: the number of sparse entries that one element writes for each matrix, KC0_SPARSE_SIZE, KCNL_SPARSE_SIZE, KG_SPARSE_SIZE, M_SPARSE_SIZE, …
*Probe	e.g. Quad4Probe: the buffers of one element in element coordinates, xe (coordinates), ue (displacements), finte (internal forces), shared by all the elements of a type and overwritten by each update_probe_* call
element	e.g. Quad4: the connectivity n1, n2, ..., the positions c1, c2, ... of its nodes in the global vectors, the positions init_k_KC0, init_k_KG, ... of its entries in the sparse arrays, the rotation matrix r11, ..., r33, and the update_* methods

The elements are Cython extension types with compiled update_* methods, such that a loop over the elements in Python costs one call per element and matrix.

Assembly of the structural matrices

The structural matrices are assembled in the coordinate (COO) format: one array of row indices, one of column indices and one of values, each with *_SPARSE_SIZE entries per element, of the types pyfe3d.INT and pyfe3d.DOUBLE. An element writes its entries at the positions init_k_KC0 to init_k_KC0 + KC0_SPARSE_SIZE - 1, so consecutive elements get init_k_KC0 incremented by KC0_SPARSE_SIZE. The entries of different elements at the same row and column are summed when the arrays are converted with scipy.sparse.coo_matrix, which is the assembly

\[\boldsymbol{K} = \sum_e \boldsymbol{L}_e^\top \boldsymbol{T}_e \boldsymbol{K}_e \boldsymbol{T}_e^\top \boldsymbol{L}_e\]

with \(\boldsymbol{K}_e\) the element matrix in element coordinates, \(\boldsymbol{T}_e\) the transformation matrix of the element, block diagonal with the matrix \(\boldsymbol{T}_{3\times 3}\) of the element, whose entries are r11, ..., r33, one block per node and per translations or rotations, see Conventions, and \(\boldsymbol{L}_e\) the Boolean matrix that extracts the dofs of the element, given by c1, c2, .... The internal force vector is assembled in place, by update_fint, into a dense vector of size N.

Entries written per element (*Data attributes)

Element	KC0	KCNL	KG	M	Others
Quad4, Quad4R	576	576	144	480	KA_BETA, KA_GAMMA, CA: 144
Tria3DSG	324	324	81	324
Tria3R	324	324	81	270
BeamC	144	144	144	144
BeamLR	144	144	36	144
Truss	72			144
Spring	72		0	0

The matrices are:

• KC0, the linear constitutive stiffness matrix \(\boldsymbol{K}_{C0}\);

• KCNL, the non-linear constitutive stiffness matrix \(\boldsymbol{K}_{CNL}(\boldsymbol{u})\), of the terms of \(\boldsymbol{B}_{NL}\), plus the geometric stiffness matrix \(\boldsymbol{K}_{GNL}(\boldsymbol{u})\) of the stress resultants of the non-linear strains, \(\boldsymbol{K}_{CNL}(\boldsymbol{u}) + \boldsymbol{K}_{GNL}(\boldsymbol{u})\), such that \(\boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL} + \boldsymbol{K}_{GNL} + \boldsymbol{K}_G\) is the tangent stiffness matrix of the geometrically non-linear analyses;

• KG, the geometric stiffness matrix \(\boldsymbol{K}_G(\boldsymbol{u})\), from the stress resultants of the displacements in the probe, or, for the shells, of given membrane stress resultants with update_KG_given_stress;

• M, the mass matrix, consistent, lumped or, for some shells, with reduced integration, mtype 0, 1 or 2 of update_M for every element, see Choosing an element;

• KA_beta, KA_gamma and CA, the aerodynamic stiffness and damping matrices of the piston theory, for the quadrilaterals.

The methods of each element are listed in Choosing an element.

Steps of an analysis

The tests of the repository, e.g. Linear static analysis of a plate, follow these steps:

1. Nodes: the array x of the coordinates and a map from the node identification numbers to their positions pos.

2. Properties: isotropic_plate(), laminated_plate() or a BeamProp.

3. Elements: one probe and one *Data per element type; for each element, set n1, n2, ..., c1 = DOF*pos1, ... and the init_k_* of the matrices needed, then call, in this order, update_rotation_matrix, update_probe_xe and the update_K* or update_M methods. The order matters: each method reads the attributes updated by the previous ones, see the class documentation of each element.

4. Assembly: KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc().

5. Boundary conditions: a Boolean array bk of the constrained dofs, e.g. bk[2::DOF] = edge for \(w = 0\) at the nodes of an edge, and bu = ~bk for the free ones. The matrices are partitioned, KC0uu = KC0[bu, :][:, bu], and the solution is expanded, u[bu] = uu. With Analysis, bu is passed instead and the partition is done by the solver, see Non-linear static analyses. The solvers of pyfe3d.solver do not remove null rows and columns: a dof left without stiffness, e.g. an unattached node, an unstiffened drilling rotation or a missing boundary condition, raises a numpy.linalg.LinAlgError naming the dofs, see check_null_rows_cols().

6. Solution: the linear static system \(\boldsymbol{K}_{{C0}_{uu}} \boldsymbol{u}_u = \boldsymbol{f}_u\), with solve() or scipy.sparse.linalg.spsolve(); the eigenvalue problems with linear_buckling(), natural_frequency() and complex_frequency(), see Linear buckling and natural frequencies; the non-linear analyses with Analysis.

7. Post-processing: update_probe_ue(u) brings the displacements of an element to element coordinates, in probe.ue; update_fint adds its internal forces, in global coordinates, to a global vector, whose values at the constrained dofs are the reactions, \(\boldsymbol{F}_{{int}_k} - \boldsymbol{f}_k\); and the strains of Quad4 are recovered with Quad4Probe.update_BL().

Since the probe is shared, update_probe_xe and update_probe_ue must be called for an element right before the methods that read the probe, e.g. update_KG, update_KCNL and update_fint, as in the loops of the examples.

Linear and non-linear problems

The linear static problem, the eigenvalue problems and the non-linear static problem use the same element matrices:

\[\boldsymbol{K}_{C0} \boldsymbol{u} = \boldsymbol{f}\]

\[(\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G(\boldsymbol{u}_0)) \boldsymbol{u} = \boldsymbol{0}, \qquad (\boldsymbol{K}_{C0} - \omega_n^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\]

\[\boldsymbol{R}(\boldsymbol{u}, \lambda) = \boldsymbol{F}_{int}(\boldsymbol{u}) - \lambda \boldsymbol{F}_{ext} = \boldsymbol{0}, \qquad \boldsymbol{K}_T = \frac{\partial \boldsymbol{R}}{\partial \boldsymbol{u}} = \boldsymbol{K}_{C0} + \boldsymbol{K}_{CNL}(\boldsymbol{u}) + \boldsymbol{K}_{GNL}(\boldsymbol{u}) + \boldsymbol{K}_G(\boldsymbol{u})\]

where \(\boldsymbol{u}_0\) is the linear static solution of the reference load of the buckling analysis and \(\boldsymbol{F}_{int}(\boldsymbol{u})\) is the vector of update_fint(..., nonlinear=1), whose exact derivative is \(\boldsymbol{K}_T\) for every element with update_KCNL (verified in tests/test_tangent_consistency.py). The linear internal force vector, nonlinear=0, is \(\boldsymbol{K}_{C0} \boldsymbol{u}\).