Assembling a model#

bfscplate2d computes the contributions of one element; the mesh, the global matrices, the boundary conditions and the loads are written by the user with NumPy and SciPy, as in the tests of the repository, see Linear static analysis.

Mesh#

A rectangular plate of sides \(a\) and \(b\) is divided into nx - 1 by ny - 1 elements of equal sides lex = a/(nx - 1) and ley = b/(ny - 1). The tests number the nodes with \(x\) varying slowest, nids_mesh = nids.reshape(nx, ny), and take the four nodes of each element from the corners of this array, in the order of The element and its degrees of freedom:

n1s = nids_mesh[:-1, :-1].flatten()
n2s = nids_mesh[1:, :-1].flatten()
n3s = nids_mesh[1:, 1:].flatten()
n4s = nids_mesh[:-1, 1:].flatten()

with nid_pos mapping each node identification number to its index, such that c1 = DOF*nid_pos[n1], and so on, and the model has N = DOF*nx*ny degrees of freedom.

Element loop and sparse arrays#

Each function update_K* and update_M() writes the entries of one element in three 1-D arrays, the row indices, the column indices and the values, at the positions init_k to init_k + SPARSE_SIZE - 1, with init_k the attribute of the element for that matrix and SPARSE_SIZE the number of entries of one element:

Function	Position attribute	Entries per element
`update_KC0()`	`init_k_KC0`	`KC0_SPARSE_SIZE`
`update_KCNL()`	`init_k_KCNL`	`KCNL_SPARSE_SIZE`
`update_KG()`, `update_KG_cte_N()`	`init_k_KG`	`KG_SPARSE_SIZE`
`update_M()`	`init_k_M`	`M_SPARSE_SIZE`
`update_KCfollower()`	`init_k_KCfollower`	`KCfollower_SPARSE_SIZE`

For num_elements elements the arrays have SPARSE_SIZE*num_elements entries, of type INT for the indices and DOUBLE for the values, and the element of index k gets init_k = k*SPARSE_SIZE, e.g. by increasing a counter by SPARSE_SIZE in the element loop:

KC0r = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
KC0c = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=INT)
KC0v = np.zeros(KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE)
init_k_KC0 = 0
for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s):
    plate = BFSCPlate2D()
    plate.c1 = DOF*nid_pos[n1]
    ...
    plate.init_k_KC0 = init_k_KC0
    update_KC0(plate, points, weights, KC0r, KC0c, KC0v)
    init_k_KC0 += KC0_SPARSE_SIZE

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

scipy.sparse.coo_matrix sums the values at repeated positions, the entries of the nodes shared by several elements. INT, DOUBLE, DOF and the sizes are given by the bfscplate2d namespace. The element writes all the pairs of degrees of freedom of each matrix, not only an upper triangle; the matrices that are symmetric come out symmetric after the sum.

The functions overwrite the indices but add the values to the existing ones. Arrays reused for a new configuration, e.g. KCNL and KG at each iteration of a nonlinear analysis, are zeroed first, KCNLv *= 0; the indices do not change for a fixed mesh.

The vectors are assembled without position attributes: update_fint() and update_fext_follower() add the vector of one element directly to a global vector of size N, at the positions of its nodes, and the functions that depend on the displacements read them from a global vector u of size N.

Loads#

A concentrated load is a value of the global force vector at a degree of freedom of a node, e.g. f[6::DOF][check] = fmid for a transverse force at the nodes selected by the boolean array check. Distributed loads are integrated with the interpolation functions of the element, the consistent load vector: a pressure \(P\) with the rows of BFSCPlate2D.update_Sw() over the integration points of the element, see Linear static analysis, and a force per unit length on an edge with BFSCPlate2D.update_Su() at \(\xi = \pm 1\), see Linear buckling analysis. The force vector of a pressure that follows the deformation is given by update_fext_follower(), see Follower pressure.

Boundary conditions#

The tests constrain degrees of freedom by partitioning the global system: a boolean array bk marks the constrained degrees of freedom, selected with the slices of The element and its degrees of freedom and the coordinates of the nodes, and bu = ~bk the unknown ones:

bk = np.zeros(N, dtype=bool)
check = isclose(x, 0) | isclose(x, a) | isclose(y, 0) | isclose(y, b)
bk[6::DOF] = check            # w = 0 at the nodes of the four edges
bk[0::DOF] = isclose(x, 0)    # u = 0 at x = 0
bk[3::DOF] = isclose(y, 0)    # v = 0 at y = 0
bu = ~bk

Kuu = KC0[bu, :][:, bu]
u = np.zeros(N)
u[bu] = spsolve(Kuu, f[bu])

The solution of the reduced system is expanded to the size of the model, with zeros at the constrained degrees of freedom, before it is given to the functions that depend on the displacements, e.g. update_KG().

Since \(w\) along an edge is interpolated from \(w\) and its derivative along the edge at the nodes, see The element and its degrees of freedom, \(w\) vanishes along the whole edge \(x\) = constant when \(w\) and \(w_{,y}\) are constrained at its nodes, positions 6 and 8, and along an edge \(y\) = constant with \(w\) and \(w_{,x}\), positions 6 and 7. A clamped edge \(x\) = constant also constrains the normal slope \(w_{,x}\) and its derivative along the edge \(w_{,xy}\), positions 7 and 9. The function boundary_flags of tests/test_follower_pressure.py builds bk in this way, for the value and the derivative along the edge of each constrained displacement, see Follower pressure.