Isoparametric 4-Node Quadrilateral: Fast Solver

01

Purpose and benchmark model

plane2d_fast is the vectorized alternative to the standard plane2d routine presented in Section 9.4. It follows the same general strategy as the CST fast solver in Section 9.2: simultaneous element computations with multiprod and direct sparse assembly with MATLAB sparse.

The benchmark uses the plane-stress cantilever of Section 9.4. The mesh has $n_H$ Q4 elements through the height and $n_L$ elements along the length, so that

02

Four Gauss-point contributions

Each Q4 finite element is integrated with four Gauss points. The number of vectorized integration contributions is therefore:

The strain-displacement and constitutive matrices are stored for every Gauss-point contribution in three-dimensional arrays:

Thus, the four terms $\boldsymbol{B}_g^T\boldsymbol{D}_g\boldsymbol{B}_g$ belonging to each finite element are evaluated vectorially. During sparse assembly, their entries have the same element connectivity and are naturally summed at the corresponding global matrix positions.

03

Vectorized element computation

The nodal coordinate arrays are first repeated for the four Gauss points of every element. The Jacobian matrices, their inverses, and the derivatives with respect to the Cartesian coordinates are then evaluated simultaneously for all $4n_{el}$ contributions.

The Q4 strain-displacement matrix has eight columns. Its public implementation therefore uses the explicit preallocation

B=zeros(3,8,nel4);

The constitutive matrices include the Gauss-point Jacobian determinant and the element thickness. The simultaneous matrix products are

DB=multiprod(DHooke,B);
kelt=reshape(multiprod(permute(B,[2 1 3]),DB),ngel^2,nel4);

Here DB stores the products $\boldsymbol{D}_g\boldsymbol{B}_g$ for all Gauss points, while each column of kelt contains the 64 entries of one $8\times8$ Gauss-point stiffness contribution.

04

Sparse global assembly

The connectivity arrays are expanded to the four Gauss-point contributions. MATLAB sparse then assembles the global stiffness matrix directly and sums coincident entries:

K=sparse(ipos1,ipos2(:),kelt(:),neq,neq);

This avoids the element-by-element stiffness-assembly loop. After applying the constraints and loads, the global linear system is solved in the standard form

S=K\F;

05

Use with Section 9.4

Start from the complete MATLAB programs provided in Section 9.4, add plane2d_fast.m and multiprod.m, and replace the call plane2d in main.m with plane2d_fast. The geometry, material, loading, constraints, and result-plotting routines remain unchanged.

06

MATLAB R2026a benchmark

The assembly interval includes element-matrix generation, loads, constraints, and global stiffness assembly. The linear solution $\boldsymbol{K}\backslash\boldsymbol{F}$ was measured separately. Total time is the sum of these two directly measured intervals.

The vectorized Gauss-point calculations and direct sparse assembly reduce the pre-solve interval by more than two orders of magnitude for these meshes. The overall acceleration is smaller because the sparse linear-system solution is essentially the same operation in both implementations.

07

Numerical equivalence

The nodal displacement vectors and recovered element-node stresses were compared before accepting the timings.

The differences are consistent with floating-point round-off caused by the different operation ordering.

08

Large-mesh comparison

For a much larger mesh with $n_H=500$ and $n_L=1250$, corresponding to 1,253,502 global equations, plane2d_fast required 10.251 s for generation and assembly and 21.082 s for the solution of the global system.

For the same number of global equations, the 4-node quadrilateral formulation requires a longer generation-and-assembly time than the CST formulation because each Q4 element contributes an $8\times8$ stiffness matrix rather than a $6\times6$ matrix. This produces more nonzero couplings and therefore a denser global stiffness matrix. The Q4 generation-and-assembly interval is approximately $10.251/2.553\approx4.02$ times longer, while the linear-system solution times for this case are of the same order.

09

Programs and downloads

multiprod.m is by Paolo de Leva and is redistributed under the included BSD-style license. Its authoritative distribution page is the MATLAB Central File Exchange entry.