Fast program for the 20-node isoparametric hexahedral finite element

01

Introduction

brick3d_20_fast and brick3d_20_fast_full are accelerated, vectorized versions of the standard HEX20 program from Section 25.3. The approach continues the vectorization of the HEX8 program presented in Section 25.2.

The first version is vectorized over all finite elements while retaining a loop over the 27 Gauss points. The fully vectorized version processes all elements and all 27 Gauss points simultaneously.

02

Vectorization over all elements

brick3d_20_fast vectorizes the operations over all elements and retains only the short loop over the 27 Gauss points. At each Gauss point, shape20 evaluates the interpolation-function derivatives. The Jacobians, inverse Jacobians and strain–displacement matrices are then computed for all elements together.

multiprod performs the page-wise matrix products, including $\htmlClass{fea-math-bold-italic}{\boldsymbol B}^T\htmlClass{fea-math-bold-italic}{\boldsymbol D}\htmlClass{fea-math-bold-italic}{\boldsymbol B}$. Global assembly uses MATLAB sparse. This version provides a good compromise between computational speed and memory requirements.

03

Full vectorization

brick3d_20_fast_full evaluates the shape-function derivatives directly in vectorized form and processes all 27 Gauss points and all elements simultaneously, without a Gauss-point loop. Large multidimensional arrays are required. This version was developed mainly as a speed/memory benchmark.

A straightforward full vectorization would require forming the $60\times60$ matrices $\htmlClass{fea-math-bold-italic}{\boldsymbol B}^T\htmlClass{fea-math-bold-italic}{\boldsymbol D}\htmlClass{fea-math-bold-italic}{\boldsymbol B}$ simultaneously for all Gauss points and all elements, leading to very large temporary arrays. To reduce memory requirements, the present implementation first integrates weighted products of the Cartesian derivatives of the shape functions, forming $20\times20$ matrices, and then assembles the $60\times60$ element stiffness matrices from the derivative-product terms, only six of which need to be computed explicitly.

04

Benchmark example

The same cantilever problem as in Section 25.1 and Section 25.3 is used. Two meshes are considered:

Once the element operations and assembly are strongly accelerated, the solution of the global linear system becomes the dominant computational cost.

05

Preconditioned conjugate gradient solution

The same preconditioned conjugate-gradient approach as in Section 25.2 is used for the larger HEX20 system. When the iterative solution is used, the following ichol + pcg sequence replaces S=K\F; in the program:

opts.type    = 'ict';
opts.droptol = 1e-3;

L = ichol(K,opts);
[S,flag,relres,iter] = pcg(K,F,1e-6,200,L,L');

06

Memory requirements

brick3d_20_fast is the preferred general-purpose version because it provides a good compromise between speed and memory usage. brick3d_20_fast_full is faster, but requires substantially more memory because all elements and all Gauss points are processed simultaneously.

For the fine mesh with 12 000 HEX20 elements, arrays such as Bt, kel, H, xyzel_all, br and the global assembly index arrays may together require several gigabytes of working memory. The estimated peak working-memory requirement is approximately 2.5–4 GB, depending on MATLAB temporary arrays and sparse-assembly overhead.

07

Programs

Use either fast routine in place of brick3d_20 with the input and plotting routines from Section 25.3. brick3d_20_fast.m requires both shape20.m and multiprod.m. brick3d_20_fast_full.m does not require shape20.m; it requires only multiprod.m. Both downloadable routines use the direct solution S=K\F; by default; for large systems, this line may be replaced by the ichol + pcg sequence shown above.

multiprod.m is third-party code by Paolo de Leva and Jinhui Bai; see Section 9.2 and Section 9.5 for source and licensing information.

The two fast HEX20 programs presented in this section were developed with assistance from Chiara.