MATLAB Implementation of C1 finite elements: Bogner-Fox-Schmit rectangle

TL;DR

This paper presents a MATLAB implementation of Bogner-Fox-Schmit C1 finite elements in 2D, emphasizing detailed construction, vectorized assembly, and visual demonstrations of energy evaluations.

Contribution

It extends MATLAB codes to implement Bogner-Fox-Schmit elements using vectorization, providing detailed construction and visualization tools.

Findings

Efficient MATLAB implementation of C1 finite elements

Demonstrations of energy evaluations and visualizations

Extension of vectorization techniques to finite element assembly

Abstract

Rahman and Valdman (2013) introduced a new vectorized way to assemble finite element matrices. We utilize underlying vectorization concepts and extend MATLAB codes to implementation of Bogner-Fox-Schmit C1 rectangular elements in 2D. Our focus is on the detailed construction of elements and simple computer demonstrations including energies evaluations and their visualizations.

Equations

$$\begin{array}{rl}& {\hat{H}}_1(\hat{x}):=2{\hat{x}}^3-3{\hat{x}}^2+1,\\ & {\hat{H}}_2(\hat{x}):=-2{\hat{x}}^3+3{\hat{x}}^2,\\ & {\hat{H}}_3(\hat{x}):={\hat{x}}^3-2{\hat{x}}^2+\hat{x},\\ & {\hat{H}}_4(\hat{x}):={\hat{x}}^3-{\hat{x}}^2\end{array}$$

$$\begin{array}{rl}& {\hat{H}}_1(0)=1,{\hat{H}}_1(1)=0,{\hat{H}}_1^{\prime }(0)=0,{\hat{H}}_1^{\prime }(0)=0,\\ & {\hat{H}}_2(0)=0,{\hat{H}}_2(1)=1,{\hat{H}}_2^{\prime }(0)=0,{\hat{H}}_2^{\prime }(0)=0,\\ & {\hat{H}}_3(0)=0,{\hat{H}}_3(1)=0,{\hat{H}}_3^{\prime }(0)=1,{\hat{H}}_3^{\prime }(0)=0,\\ & {\hat{H}}_4(0)=0,{\hat{H}}_4(1)=0,{\hat{H}}_4^{\prime }(0)=0,{\hat{H}}_4^{\prime }(0)=1,\end{array}$$

$$\begin{array}{r}{H}_1(x):={\hat{H}}_1(\hat{x}(x)),\\ H_2(x):={\hat{H}}_2(\hat{x}(x)),\\ H_3(\hat{x}):=h{\hat{H}}_3(\hat{x}(x)),\\ H_4(\hat{x}):=h{\hat{H}}_4(\hat{x}(x))\end{array}$$

$$h:=b-a.$$

$$\begin{array}{rl}& H_1^{\prime }(x)={\hat{H}}_1^{\prime }(\hat{x})/h,H_1^{\prime \prime }(x)={\hat{H}}_1^{\prime \prime }(\hat{x})/h^2,\\ & H_2^{\prime }(x)={\hat{H}}_2^{\prime }(\hat{x})/h,H_2^{\prime \prime }(x)={\hat{H}}_2^{\prime \prime }(\hat{x})/h^2,\\ & H_3^{\prime }(x)={\hat{H}}_3^{\prime }(\hat{x}),H_3^{\prime \prime }(x)={\hat{H}}_3^{\prime \prime }(\hat{x})/h,\\ & H_4^{\prime }(x)={\hat{H}}_4^{\prime }(\hat{x}),H_4^{\prime \prime }(x)={\hat{H}}_4^{\prime \prime }(\hat{x})/h.\end{array}$$

$${\tilde{\phi }}_{j,k}(\hat{x},\hat{y}):={\hat{H}}_j(\hat{x}){\hat{H}}_k(\hat{y}),j,k=1,\dots ,4$$

$${\hat{\phi }}_i(\hat{x},\hat{y}):={\tilde{\phi }}_{j_i,k_i}(\hat{x},\hat{y}),i=1,\dots ,16,$$

$$\begin{array}{r}(j_i,k_i{)}_{i=1}^{16}=\begin{array}{c}\{(1,1),(2,1),(2,2),(1,2),\\ (3,1),(4,1),(4,2),(3,2),\\ (1,3),(2,3),(2,4),(1,4),\\ (3,3),(4,3),(4,4),(3,4)\}.\end{array}\end{array}$$

$$v(\hat{x},\hat{y})=\sum _{i=1}^{16}{v}_i{\hat{\phi }}_i(\hat{x},\hat{y}),$$

$${\hat{N}}_1:=[0,0],{\hat{N}}_2:=[1,0],{\hat{N}}_3:=[1,1],{\hat{N}}_4:=[0,1].$$

$$(\hat{x},\hat{y})(x,y):=((x-a)/h_x,(y-c)/h_y),$$

$$h_x:=b-a,h_y:=d-c.$$

$${\phi }_i(x,y):={\hat{H}}_{j_i}\left(\frac{x-a}{h_x}\right){\hat{H}}_{k_i}\left(\frac{y-c}{h_y}\right),i=1,\dots ,16,$$

$$\frac{\partial {\phi }_i}{\partial x},\frac{\partial {\phi }_i}{\partial y},\frac{{\partial }^2{\phi }_i}{\partial x^2},\frac{{\partial }^2{\phi }_i}{\partial y^2},\frac{{\partial }^2{\phi }_i}{\partial x\partial y},i=1,\dots ,16$$

$$VC1=\left(\begin{array}{cccc}v(N_1),& \frac{\partial v}{\partial x}(N_1),& \frac{\partial v}{\partial y}(N_1),& \frac{{\partial }^{2}v}{\partial x\partial y}(N_1)\\ ⋮& ⋮& ⋮& ⋮\\ v(N_n),& \frac{\partial v}{\partial x}(N_n),& \frac{\partial v}{\partial y}(N_n),& \frac{{\partial }^{2}v}{\partial x\partial y}(N_n)\end{array}\right)$$

$$v,\frac{\partial v}{\partial x},\frac{\partial v}{\partial y},\frac{{\partial }^{2}v}{\partial x\partial y}$$

$$\frac{{\partial }^{2}v}{\partial x^2},\frac{{\partial }^{2}v}{\partial y^2}.$$

$$v(x,y)=(1-x^2{)}^2(1-y^2{)}^2$$

$$\frac{\partial v}{\partial x}(x,y)=-4x(1-x^2)(1-y^2{)}^2$$

$$\frac{\partial v}{\partial x}(x,y)=-4y(1-x^2{)}^2(1-y^2)$$

$$\frac{{\partial }^{2}v}{\partial x\partial y}(x,y)=16xy(1-x^2)(1-y^2)$$

$$||v|{|}^2:$$

$$||\nabla v|{|}^2:$$

$$||{\nabla }^{2}v|{|}^2:$$

$$(f,v):$$

$$(||v|{|}^2+||\nabla v|{|}^2+||{\nabla }^{2}v|{|}^2{)}^{1/2}$$

$$||v|{|}^2$$

$$||\nabla v|{|}^2$$

$$||{\nabla }^{2}v|{|}^2$$

$$(f,v)$$

Full text

1 Introduction

Boundary problems with fourth order elliptic operators [3] appear in many applications including thin beams and plates and strain gradient elasticity [5]. Weak formulations and implementations of these problems require H2-conforming finite elements, leading to $C^{1}$ continuity (of functions as well as of their gradients) of approximations over elements edges. This continuity condition is generally technical to achieve and few types of finite elements are known to quarantee it. We consider probably the simplest of them, the well known Bogner-Fox-Schmit rectangle [2], i.e., a rectangular $C^{1}$ element in two space dimensions.

We are primarily interested in explaining the construction of BFS elements, their practical visualization and evalutions. Our MATLAB implementation is based on codes from [1, 9, 6]. The main focus of these papers were assemblies of finite element matrices and local element matrices were computed all at once by array operations and stored in multi-dimentional arrays (matrices). Here, we utilize underlying vectorization concepts without the particular interest in corresponding FEM matrices. More details on our recent implementations of $C^{1}$ models in nonlinear elastic models of solids can be found in [7, 4, 8]. The complementary software to this paper is available at

https://www.mathworks.com/matlabcentral/fileexchange/71346

for download and testing.

2 Construction of $C^{1}$ finite elements

2.1 Hermite elements in 1D

We define four cubic polynomials

[TABLE]

over a reference interval $\hat{I}:=[0,1]$ and can easily check the conditions:

[TABLE]

so only one value or derivative is equal to 1 and all other three values are equal to 0. These cubic functions create a finite element basis on $\hat{I}$.

More generally, we define

[TABLE]

for $x\in I:=[a,b]$, where $\hat{x}(x):=(x-a)/h$ is an affine mapping from $I$ to $\hat{I}$ and $h$ denotes the interval $I$ size

[TABLE]

These functions are also cubic polynomials and satisfy again the conditions (2) with function arguments $0,1$ replaced by $a,b$. They create actual finite element basis which ensures $C^{1}$ continuity of finite element approximations. The chain rule provides higher order derivatives:

[TABLE]

Example 1

Example of basis functions defined on reference and actual intervals are shown in Figure 1 and pictures can be reproduced by

draw_C1basis_1D

script located in the main folder.

2.2 Bogner-Fox-Schmit rectangular element in 2D

Products of functions

[TABLE]

define 16 Bogner-Fox-Schmit (BFS) basis functions on a reference rectangle $\hat{R}:=[0,1]\times[0,1]$. For practical implementations, we reorder them as

[TABLE]

where sub-indices are ordered in a sequence

[TABLE]

With this ordering, a finite element approximation $v\in C^{1}(\hat{R})$ rewrites as a linear combination

[TABLE]

where:

• •

coefficients $v_{1},\dots,v_{4}$ specify values of $v$,

• •

coefficients $v_{5},\dots,v_{8}$ specify values of $\frac{\partial v}{\partial x}$,

• •

coefficients $v_{9},\dots,v_{12}$ specify values of $\frac{\partial v}{\partial y}$,

• •

coefficients $v_{13},\dots,v_{16}$ specify values of $\frac{\partial^{2}v}{\partial x\partial y}$

at nodes

[TABLE]

Example 2

Four (out of 16) reference basis functions corresponding to the node $\hat{N}_{2}$ are shown in Figure 2 and pictures can be reproduced by

draw_C1basis_2D

script located in the main folder.

2.2.1 More implementation details on functions evaluations

For a general rectangle $R:=[a,b]\times[c,d]$, we define an affine mapping

[TABLE]

from $R$ to $\hat{R}$, where the rectangular lengths are

[TABLE]

It enables us to define BFS basis functions on $R$ as

[TABLE]

Based on (4), higher order derivatives up to the second order,

[TABLE]

can be derived as well. All basis functions (7) are evaluated by the function

shapefun(points’,etype,h)

and their derivatives (8) by the function

shapeder(points’,etype,h)

For BFS elements, we have to set etype=’C1’ and a vector of rectangular lengths h=[hx,hy]. The matrix points then contains a set of points $\hat{x}\in\hat{R}$ in a reference element at which functions are evaluated. Both functions are vectorized and their outputs are stored as vectors, matrices or three-dimensional arrays.

Example 3

The command

[shape]=shapefun([0.5 0.5]’,’C1’,[1 1])

returns a (column) vector shape $\in\mathbb{R}^{16\times 1}$ of all BFS basis function defined on $\hat{R}:=[0,1]\times[0,1]$ and evaluated in the rectangular midpoint $[0.5,0.5]\in\hat{R}$. The command

[dshape]=shapeder([0.5 0.5]’,’C1’,[2 3])

returns a three-dimensional array dshape $\in\mathbb{R}^{16\times 1\times 5}$ of all derivatives up to the second order of all BFS basis function defined on a general rectangle with lengths $h_{x}=2$ and $h_{y}=3$ and evaluated in the rectangular midpoint $[0.5,0.5]\in\hat{R}$.

For instance, if $R:=[1,3]\times[2,5]$, values of all derivatives are evaluated in the rectangular midpoint $[2,3.5]\in R$.

More generally, if points $\in\mathbb{R}^{np\times 2}$ consists of $np>1$ points, then shapefun return a matrix of size $\mathbb{R}^{16\times np}$ and shapeder returns a three-dimensional array of size $\mathbb{R}^{16\times np\times 5}$.

2.3 Representation of $C^{1}$ functions and their visualization

Let us assume a triangulation $\mathcal{T}(\Omega)$ into rectangles of a domain $\Omega$. In correspondence with our implementation, we additionally assume that all rectangles are of the same size, i.e., with lengths $h_{x},h_{y}>0$. Examples of $\mathcal{T}(\Omega)$ are given in Figure 3.

Let $\mathcal{N}$ denotes the set of all rectangular nodes and $|\mathcal{N}|:=n$ the number of them. A $C^{1}$ function $v\in\mathcal{T}(\Omega)$ is represented in BSF basis by a matrix

[TABLE]

containing values of $v,\frac{\partial v}{\partial x},\frac{\partial v}{\partial y},\frac{\partial^{2}v}{\partial x\partial y}$ in all nodes of $\mathcal{T}(\Omega)$. In the spirit of the finite element method, values of $v$ on each rectangle $T\in\mathcal{T}(\Omega)$ are obtained by an affine mapping to the reference element $\hat{R}$. Our implementations allows to evaluate and visualize continuous fields

[TABLE]

and also two additional (generally discontinuous) fields

[TABLE]

It is easy to evaluate a $C^{1}$ function in a particular element point for all elements (rectangles) at once. A simple matrix-matrix MATLAB multiplication

Vfun=VC1_elems*shape

where a matrix VC1_elems $\in\mathbb{R}^{ne\times 16}$ contains in each row all 16 coefficients (taken from VC1) corresponding to each element ($ne$ denotes a number of elements) returns a matrix Vfun $\in\mathbb{R}^{ne\times np}$ containing all function values in all elements and all points. Alternate multiplications

V1=VC1_elemssqueeze(dshape(:,1,:)); % Dx V2=VC1_elemssqueeze(dshape(:,2,:)); % Dy V11=VC1_elemssqueeze(dshape(:,3,:)); % Dxx V22=VC1_elemssqueeze(dshape(:,4,:)); % Dyy V12=VC1_elems*squeeze(dshape(:,5,:)); % Dxy

return matrices V1, V2, V11, V22, V12 $\in\mathbb{R}^{ne\times np}$ containing values of all derivatives up to the second order in all elements and all points. A modification for evaluation of function values at particular edges points is also available and essential for instance for models with energies formulated on boundary edges [8].

Example 4

We consider a function

[TABLE]

on the domain $\Omega=(-1,1)^{2}$. This function was also used in [4] as an initial vertical displacement in a time-dependent simulation of viscous von Kármán plates.

To represent $v$ in terms of BFS elements, we additionally need to know values of

[TABLE]

in nodes of a rectangular mesh $\mathcal{T}(\Omega)$. The function and its derivatives up to the second order are shown in Figure 4 and its values in elements and edges midpoints in Figure 5.

All pictures can be reproduced by

draw_C1example_2D

script located in the main folder.

2.4 Evaluation of $C^{1}$ functions and their numerical integration

Various energy formulations include evaluations of integrals of the types

[TABLE]

where $v\in H^{2}(\Omega)$ and $f\in L^{2}(\Omega)$ is given. The expression

[TABLE]

then defines the full norm in the Sobolev space $H^{2}(\Omega)$. For $v$ represented in the BFS basis we can evaluate above mentioned integrals numerically by quadrature rules. Our implementation provides three different rules with 1, 4 or 9 Gauss points. Each quadrature rule is defined by coordinates of all Gauss points and their weights.

Example 5

Gauss points are displayed in Figure 6 and all pictures can be reproduced by

draw_ips

script located in the main folder.

Example 6

An analytical integration for the function $v$ from (9) and $f=x^{2}y^{2}$ reveals that

[TABLE]

for a domain $\Omega=(-1,1)^{2}$. We consider a sequence of uniformly refined meshes with levels 1-10:

• •

the coarsest (level=1) mesh with 9 nodes and 4 elements is shown in Figure 6,

• •

a finer (level=4) mesh with 289 nodes and 256 elements is shown in Figure 3 (left),

• •

the finest (level=10) mesh consists of 1.050.625 nodes and 1.048.576 elements, not shown here.

Figure 7 the depicts the convergence of numerical quadratures to the exact values above. We notice that all three quadrature rules provide the same rates of convergence. The only exception is the evaluation of the second gradient integral $||\nabla^{2}v||^{2}$, where the numerical quadrature using 1 Gauss point deteriorates the convergence.

All pictures can be reproduced by

start_integrate

script located in the main folder.