FEM2D module¶

Author:	Francois Cuvelier <cuvelier@math.univ-paris13.fr>
Date:	15/09/2013

Contains functions to build some finite element matrices using $P_1$ -Lagrange finite elements on a 2D mesh. Each assembly matrix is computed by three differents versions called base, OptV1 and OptV2 (see here)

Contents

Assembly matrices (versions base, OptV1 and OptV2)
Elementary matrices (used by versions base and OptV1)
Vectorized tools (used by version OptV2)
- Vectorized computation of basis functions gradients
Vectorized elementary matrices (used by version OptV2)
Mesh

Assembly matrices (versions `base`, `OptV1` and `OptV2`)¶

Let ${\cal T}_h$ be a triangular mesh of $\Omega.$ We note $\Omega_h=\bigcup_{T_k \in {\cal T}_h} T_k$ a triangulation of $\Omega$ with the following structure data:

$\mbox{\begin{tabular}{lccll} \hline \textbf{name} & \textbf{type} & \textbf{dimension} & \textbf{description}& \textbf{Python}\\ \hline $\nq$ & integer & 1 & number of vertices& \texttt{nq}\\ $\nme$ & integer & 1 & number of elements& \texttt{nme}\\ $\q$ & double & $2 \times \nq$ & \begin{minipage}[t]{7.9cm} array of vertices coordinates. $\q(\nu,j)$ is the $\nu$-th coordinate of the $j$-th vertex, $\nu\in\{1,2\}$, $j\in\{1,\hdots,\rm{n_q}\}.$ The $j$-th vertex will be also denoted by $\rm{q}^j$ \end{minipage}& \begin{minipage}[t]{3cm} \texttt{q} (transposed)\\ \texttt{q[j-1]} = $\q^j$ \end{minipage}\\ $\me$ & integer & $3 \times \nme$ & \begin{minipage}[t]{7.9cm} connectivity array. $\me(\beta,k)$ is the storage index of the $\beta$-th vertex of the $k$-th element, in the array~$q$, for $\beta\in\{1,2,3\}$ and $k\in\{1,\hdots,{\nme}\}$ \end{minipage}&\texttt{me} (transposed)\\ $\rm areas$ & double & $1\times {\nme}$ & \begin{minipage}[t]{7.9cm} array of areas. ${\rm areas}(k)$ is the $k$-th triangle area, $k\in\{1,\hdots,{\nme}\}$ \end{minipage}&\texttt{areas}\\ \hline \end{tabular}}$

The $P_1$ -Lagrange basis functions associated with $\Omega_h$ are noted $\varphi_i$ for all $i\in\{1,\hdots,\rm{n_q}\}$ and are defined by

$\varphi_i(\rm{q}^j)=\delta_{i,j},\ \ \forall (i,j)\in\{1,\hdots,\rm{n_q}\}^2$

We also define the global alternate basis $\mathcal{B}_a$ by

$\mathcal{B}_a=\{\pmb{\psi}_1,\hdots,\pmb{\psi}_{2\rm{n_q}}\}=\left\{ \begin{pmatrix} \varphi_1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_1 \end{pmatrix}, \begin{pmatrix} \varphi_2 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_2\end{pmatrix}, \hdots, \begin{pmatrix} \varphi_{\rm n_q} \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_{\rm n_q} \end{pmatrix} \right\}$

and the global block basis $\mathcal{B}_b$ by

$\mathcal{B}_b=\{\pmb{{\phi}}_1,\hdots,\pmb{{\phi}}_{2\rm{n_q}}\}=\left\{ \begin{pmatrix} \varphi_1 \\ 0 \end{pmatrix}, \begin{pmatrix} \varphi_2 \\ 0 \end{pmatrix},\hdots, \begin{pmatrix} \varphi_{\rm n_q} \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_1 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_2\end{pmatrix},\hdots, \begin{pmatrix} 0 \\ \varphi_{\rm n_q} \end{pmatrix} \right\}.$

Mass Matrix ¶

Assembly of the Mass Matrix by $P_1$ -Lagrange finite elements using respectively version base, OptV1 and OptV2 (see report). The Mass Matrix $\mathbb{M}$ is given by

$\mathbb{M}_{i,j}=\int_{\Omega_h} \varphi_i(q)\ \varphi_j(q)dq,\ \ \forall (i,j) \in \{1,...,\rm{n_q}\}^2.$

Note

generic syntax:

M = MassAssembling2DP1<version>(nq,nme,me,areas)

nq: total number of nodes of the mesh, also denoted by $\nq$ ,
nme: total number of triangles, also denoted by $\nme$ ,
me: Connectivity array, (nme,3) array,
areas: Array of areas, (nme,) array,
return a Scipy CSC sparse matrix of size $\nq \times \nq$

where <version> is base, OptV1 or OptV2

Benchmarks of theses functions are presented in Mass Matrix. We give a simple usage :

>>> from pyOptFEM.FEM2D import *
>>> Th=SquareMesh(5)
>>> Mbase =  MassAssembling2DP1base(Th.nq,Th.nme,Th.me,Th.areas)
>>> MOptV1= MassAssembling2DP1OptV1(Th.nq,Th.nme,Th.me,Th.areas)
>>> print(" NormInf(Mbase-MOptV1)=%e " % NormInf(Mbase-MOptV1))
NormInf(Mbase-MOptV1)=6.938894e-18
>>> MOptV2= MassAssembling2DP1OptV2(Th.nq,Th.nme,Th.me,Th.areas)
>>> print(" NormInf(Mbase-MOptV2)=%e " % NormInf(Mbase-MOptV2))
NormInf(Mbase-MOptV2)=6.938894e-18

We can show sparsity of the Mass matrix :

>>> from pyOptFEM.FEM2D import *
>>> Th=SquareMesh(20)
>>> M=MassAssembling2DP1base(Th.nq,Th.nme,Th.me,Th.areas)
>>> showSparsity(M)
Figure 123: Sparsity of Mass Matrix generated with command showSparsity(M)

Note

sources code

pyOptFEM.FEM2D.assembly.MassAssembling2DP1base(nq, nme, me, areas)[source]: Assembly of the Mass Matrix by $P_1$ -Lagrange finite elements using base version (see report).

pyOptFEM.FEM2D.assembly.MassAssembling2DP1OptV1(nq, nme, me, areas)[source]: Assembly of the Mass Matrix by $P_1$ -Lagrange finite elements using OptV1 version (see report).

pyOptFEM.FEM2D.assembly.MassAssembling2DP1OptV2(nq, nme, me, areas)[source]: Assembly of the Mass Matrix by $P_1$ -Lagrange finite elements using OptV2 version (see report).

Stiffness Matrix ¶

Assembly of the Stiffness Matrix by $P_1$ -Lagrange finite elements using respectively version base, OptV1 and OptV2 (see report). The Stiff Matrix $\mathbb{S}$ is given by

$\mathbb{S}_{i,j}=\int_{\Omega_h} \nabla\varphi_i(q).\nabla \varphi_j(q)dq,\ \ \forall (i,j) \in \{1,...,\rm{n_q}\}^2.$

Note

generic syntax:

M = StiffAssembling2DP1<version>(nq,nme,q,me,areas)

nq: total number of nodes of the mesh, also denoted by $\nq$ ,
nme: total number of triangles, also denoted by $\nme$ ,
q: Array of vertices coordinates, (nq,2) array
me: Connectivity array, (nme,3) array,
areas: Array of areas, (nme,) array,
return a Scipy CSC sparse matrix of size $\nq \times \nq$

where <version> is base, OptV1 or OptV2

Benchmarks of theses functions are presented in Stiffness Matrix. We give a simple usage :

>>> pyOptFEM.FEM2D import *
>>> Th=SquareMesh(5)
>>> Sbase =  StiffAssembling2DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.areas)
>>> SOptV1= StiffAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas)
>>> print(" NormInf(Sbase-SOptV1)=%e " % NormInf(Sbase-SOptV1))
NormInf(Sbase-SOptV1)=0.000000e+00
>>> SOptV2= StiffAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas)
>>> print(" NormInf(Sbase-SOptV2)=%e " % NormInf(Sbase-SOptV2))
NormInf(Sbase-SOptV1)=4.440892e-16

Note

sources code

pyOptFEM.FEM2D.assembly.StiffAssembling2DP1base(nq, nme, q, me, areas)[source]: Assembly of the Stiff Matrix by $P_1$ -Lagrange finite elements using base version (see report).

pyOptFEM.FEM2D.assembly.StiffAssembling2DP1OptV1(nq, nme, q, me, areas)[source]: Assembly of the Stiff Matrix by $P_1$ -Lagrange finite elements using OptV1 version (see report).

pyOptFEM.FEM2D.assembly.StiffAssembling2DP1OptV2(nq, nme, q, me, areas)[source]: Assembly of the Stiff Matrix by $P_1$ -Lagrange finite elements using OptV2 version (see report).

Stiffness Elasticity Matrix ¶

Assembly of the Stiffness Elasticity Matrix by $P_1$ -Lagrange finite elements using respectively version base, OptV1 and OptV2 (see report). The Stiffness Elasticity Matrix $\mathbb{K}$ is given by

$\mathbb{K}_{m,l}=\int_{\Omega_h} \underline{\pmb{\epsilon}}^t(\psi_m(q)) \underline{\pmb{\sigma}}(\psi_l(q)), \ \ \forall (m,l) \in \{1,...,2\rm{n_q}\}^2.$

where $\underline{\pmb{\sigma}}=(\sigma_{xx},\sigma_{yy},\sigma_{xy})^t$ and $\underline{\pmb{\epsilon}}=(\epsilon_{xx},\epsilon_{yy},\epsilon_{xy})^t$ are the elastic stress and strain tensors respectively.

Note

generic syntax:

M = StiffElasAssembling2DP1<version>(nq,nme,q,me,areas,la,mu,Num)

nq: total number of nodes of the mesh, also denoted by $\nq$ ,
nme: total number of triangles, also denoted by $\nme$ ,
q: array of vertices coordinates,
- (nq,2) array for base and OptV1,
- (2,nq) array for OptV2 version,
me: Connectivity array,
- (nme,3) array for base and OptV1,
- (3,nme) array for OptV2 version,
areas: (nme,) array of areas,
la: the first Lame coefficient in Hooke’s law, denoted by $\lambda$ ,
mu: the second Lame coefficient in Hooke’s law, denoted by $\mu$ ,
Num:
- 0: global alternate numbering with local alternate numbering (classical method),
- 1: global block numbering with local alternate numbering,
- 2: global alternate numbering with local block numbering,
- 3: global block numbering with local block numbering.
return a Scipy CSC sparse matrix of size $2\nq \times 2\nq$

where <version> is base, OptV1 or OptV2

Benchmarks of theses functions are presented in Stiffness Elasticity Matrix. We give a simple usage :

>>> from pyOptFEM.FEM2D import *
>>> Th=SquareMesh(5)
>>> Kbase =  StiffElasAssembling2DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.areas,2,0.5,0)
>>> KOptV1= StiffElasAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas,2,0.5,0)
>>> print(" NormInf(Kbase-KOptV1)=%e " % NormInf(Kbase-KOptV1))
NormInf(Kbase-KOptV1)=8.881784e-16
>>> KOptV2= StiffElasAssembling2DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.areas,2,0.5,0)
>>> print(" NormInf(Kbase-KOptV2)=%e " % NormInf(Kbase-KOptV2))
NormInf(Kbase-KOptV2)=1.776357e-15

We now illustrate the consequences of the choice of the global basis on matrix sparsity

global alternate basis $\mathcal{B}_a$ (Num=0 or Num=2)

>>> from pyOptFEM.FEM2D import *
>>> Th=SquareMesh(15)
>>> K0=StiffElasAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas,2,0.5,0)
>>> showSparsity(K0)

Figure 124: Sparsity of Stiffness Elasticity Matrix generated with global alternate numbering (Num=0 or 2)

>>> K3=StiffElasAssembling2DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.areas,2,0.5,3)
>>> showSparsity(K3)

global block basis $\mathcal{B}_a$ (Num=1 or Num=3)

Figure 125: Sparsity of Stiffness Elasticity Matrix generated with global block numbering (Num=1 or 3)

Note

sources code

pyOptFEM.FEM2D.assembly.StiffElasAssembling2DP1base(nq, nme, q, me, areas, la, mu, Num)[source]: Assembly of the Stiffness Elasticity Matrix by $P_1$ -Lagrange finite elements using OptV2 version (see report).

pyOptFEM.FEM2D.assembly.StiffElasAssembling2DP1OptV1(nq, nme, q, me, areas, la, mu, Num)[source]: Assembly of the Stiffness Elasticity Matrix by $P_1$ -Lagrange finite elements using OptV1 version (see report).

pyOptFEM.FEM2D.assembly.StiffElasAssembling2DP1OptV2(nq, nme, q, me, areas, la, mu, Num)[source]: Assembly of the Stiffness Elasticity Matrix by $P_1$ -Lagrange finite elements using OptV2 version (see report).

Elementary matrices (used by versions `base` and `OptV1`)¶

Let $T$ be a triangle, of area $|T|$ and with $\pmb{q}^1$ , $\pmb{q}^2$ and $\pmb{q}^3$ its three vertices. We denote by $\tilde{\varphi}_1$ , $\tilde{\varphi}_2$ and $\tilde{\varphi}_3$ the $P_1$ -Lagrange local basis functions such that $\tilde{\varphi}_i(\pmb{q}^j)=\delta_{i,j},\ \forall(i,j)\in\{1,2,3\}^2$ .

We also define the local alternate basis $\tilde{\mathcal{B}}_a$ by

$\tilde{\mathcal{B}}_a=\{\pmb{\tilde{\psi}}_1,\hdots,\pmb{\tilde{\psi}}_6\}=\left\{ \begin{pmatrix} \tilde{\varphi}_1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_1 \end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_2 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_2\end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_3 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_3 \end{pmatrix} \right\}$

and the local block basis $\tilde{\mathcal{B}}_b$ by

$\tilde{\mathcal{B}}_b=\{\pmb{\tilde{\phi}}_1,\hdots,\pmb{\tilde{\phi}}_6\}=\left\{ \begin{pmatrix} \tilde{\varphi}_1 \\ 0 \end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_2 \\ 0 \end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_3 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_1 \end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_2\end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_3 \end{pmatrix} \right\}.$

The elasticity tensor, $\mathbb{H}$ , obtained from Hooke’s law with an isotropic material, defined with the Lamé parameters $\lambda$ and $\mu$ is given by

$\mathbb{H} =\begin{pmatrix} \lambda+2\mu & \lambda & 0\\ \lambda & \lambda+2\mu & 0\\ 0 & 0 & \mu \end{pmatrix}$

and, for a function $\vecb{u}=(u_1,u_2)$ the strain tensors is given by

$\Odv(\vecb{u})=\begin{pmatrix}\DP{u_1}{x}, & \DP{u_2}{y}, \DP{u_2}{x}+\DP{u_1}{y} \end{pmatrix}^t$

Element Mass Matrix ¶

The element Mass matrix, $\mathbb{M}^e(T)$ , for the triangle $T$ , is defined by

$\mathbb{M}^e_{i,j}(T)=\int_T \tilde{\varphi}_i(\pmb{q}) \tilde{\varphi}_j(\pmb{q})\ d\pmb{q},\ \ \forall(i,j)\in\{1,2,3\}^2$

We obtain :

(1) $\mathbb{M}^e(T) =\frac{|T|}{12}\begin{pmatrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{pmatrix}$

Note

pyOptFEM.FEM2D.elemMatrix.ElemMassMat2DP1(area)[source]

Compute element Mass matrix, $\mathbb{M}^e(T)$ for a given triangle $T$ of area $|T|$

Parameters:	area (float) – area of the triangle.
Returns:	$3 \times 3$ numpy array of floats.

Element Stiffness Matrix ¶

The element Stiffness matrix, $\mathbb{S}^e(T)$ , for the $T$ is defined by

$\mathbb{S}^e_{i,j}(T)=\int_T <\nabla \tilde{\varphi}_i(\pmb{q}),\nabla \tilde{\varphi}_j(\pmb{q})>\ d\pmb{q},\ \ (i,j)\in\{1,2,3\}^2$

We have :

$\mathbb{S}^e(T)= \frac{1}{4|T|} \begin{pmatrix} <\pmb{u},\pmb{u}> & <\pmb{u},\pmb{v}> &<\pmb{u},\pmb{w}>\\ <\pmb{v},\pmb{u}> & <\pmb{v},\pmb{v}> &<\pmb{v},\pmb{w}>\\ <\pmb{w},\pmb{u}> & <\pmb{w},\pmb{v}> &<\pmb{w},\pmb{w}> \end{pmatrix}.$

where $\pmb{u}=\pmb{q}^2-\pmb{q}^3$ , $\pmb{v}=\pmb{q}^3-\pmb{q}^1$ and $\pmb{w}=\pmb{q}^1-\pmb{q}^2$ .

Note

pyOptFEM.FEM2D.elemMatrix.ElemStiffMat2DP1(q1, q2, q3, area)[source]

Compute element Stiffness matrix, $\mathbb{S}^e(T)$ for a given triangle $T$

Parameters:	q1,q2,q3 ( $2 \times 1$ numpy array) – the three vertices of the triangle, area (float) – area of the triangle.
Returns:
Type :	$3 \times 3$ numpy array of floats.

Element Stiffness Elasticity Matrix ¶

Let $\pmb{u}=\pmb{q}^2-\pmb{q}^3$ , $\pmb{v}=\pmb{q}^3-\pmb{q}^1$ and $\pmb{w}=\pmb{q}^1-\pmb{q}^2$ .

The element Stiffness Elasticity matrix, $\mathbb{K}^e(T)$ , for a given triangle $T$ in the local alternate basis $\tilde{\mathcal{B}}_a$ is defined by

$\mathbb{K}^e_{i,j}(T) = \int_T \underline{\pmb{\epsilon}}(\tilde{\pmb{\psi}}_j)^t(\pmb{q}) \mathbb{H} \underline{\pmb{\epsilon}}(\tilde{\pmb{\psi}}_i)(\pmb{q})\, d\pmb{q}, \ \ \ \forall (i,j)\in\{1,\hdots,6\}^2.$

We also have

$\mathbb{K}^e(T)=\frac{1}{4|T|} \mathbb{B}^t\mathbb{H}\mathbb{B}$

where $\mathbb{H}$ is the elasticity tensor and

$\mathbb{B}=\begin{pmatrix} u_2 & 0 & v_2 & 0 & w_2 & 0\\ 0 & -u_1 & 0 & -v_1 & 0 & -w_1\\ -u_1 & u_2 & -v_1 & v_2 & -w_1 & w_2 \end{pmatrix}$

Note

pyOptFEM.FEM2D.elemMatrix.ElemStiffElasMat2DP1Ba(ql, area, H)[source]

Return the element Stiffness Elasticity matrix, $\mathbb{K}^e(T)$ , for a given triangle $T$ in the local alternate basis $\mathcal{B}_a$

Parameters:	ql ( $3 \times 2$ numpy array) – contains the three vertices of the triangle : `ql[0]`, `ql[1]` and `ql[2]`, area (float) – area of the triangle , H ( $3 \times 3$ numpy array) – Elasticity tensor, $\mathbb{H}$ .
Returns:	$\mathbb{K}^e(T)$ in $\mathcal{B}_a$ basis.
Type :	$6 \times 6$ numpy array of floats.

The element Stiffness Elasticity matrix, $\mathbb{K}^e(T)$ , for a given triangle $T$ in the local block basis $\tilde{\mathcal{B}}_b$ is defined by

$\mathbb{K}^e_{i,j}(T) = \int_T \underline{\pmb{\epsilon}}(\tilde{\pmb{\phi}}_j)^t(\pmb{q}) \mathbb{H} \underline{\pmb{\epsilon}}(\tilde{\pmb{\phi}}_i)(\pmb{q})\, d\pmb{q}, \ \ \ \forall (i,j)\in\{1,\hdots,6\}^2.$

We also have

$\mathbb{K}^e(T)=\frac{1}{4|T|} \mathbb{B}^t\mathbb{H}\mathbb{B}$

where $\mathbb{H}$ is the elasticity tensor and

$\mathbb{B}=\begin{pmatrix} u_2 & v_2 & w_2 & 0 & 0 & 0\\ 0 & 0 & 0 & -u_1 & -v_1 & -w_1\\ -u_1 & -v_1 & -w_1 & u_2 & v_2 & w_2 \end{pmatrix}$

Note

pyOptFEM.FEM2D.elemMatrix.ElemStiffElasMat2DP1Bb(ql, area, H)[source]

Return the element Stiffness Elasticity matrix, $\mathbb{K}^e(T)$ , for a given triangle $T$ in the local block basis $\mathcal{B}_b$

Parameters:	ql ( $3 \times 2$ numpy array) – contains the three vertices of the triangle : `ql[0]`, `ql[1]` and `ql[2]`, area (float) – area of the triangle, H ( $3 \times 3$ numpy array) – Elasticity tensor, $\mathbb{H}$ .
Returns:	$\mathbb{K}^e(T)$ in $\mathcal{B}_b$ basis.
Type :	$6 \times 6$ numpy array of floats

Vectorized tools (used by version `OptV2`)¶

Vectorized computation of basis functions gradients ¶

By construction, the gradients of basis functions are constants on each element $T_k.$ So, we denote, $\forall \il\in\ENS{1}{3},$ by $\vecb{G}^\il$ the $2 \times \nme$ array defined, $\forall k\in\ENS{1}{\nme},$ by

$\vecb{G}^\il(:,k)= \GRAD\BasisFunc_{\me(\il,k)}(\q),\ \forall \q\in T_k.$

On a triangle $T_k,$ we denote by $\vecb{D}^{12}=\q^{\me(1,k)}-\q^{\me(2,k)},$ $\vecb{D}^{13}=\q^{\me(1,k)}-\q^{\me(3,k)}$ and $\vecb{D}^{23}=\q^{\me(2,k)}-\q^{\me(3,k)}.$ Then, we have

$\GRAD\BasisFunc_{1}^k(\q)=\frac{1}{2|T_k|} \begin{pmatrix} \vecb{D}^{23}_y\\-\vecb{D}^{23}_x \end{pmatrix},\ \GRAD\BasisFunc_{2}^k(\q)=\frac{1}{2|T_k|} \begin{pmatrix} -\vecb{D}^{13}_y\\\vecb{D}^{13}_x \end{pmatrix},\ \GRAD\BasisFunc_{3}^k(\q)=\frac{1}{2|T_k|} \begin{pmatrix} \vecb{D}^{12}_y\\-\vecb{D}^{12}_x \end{pmatrix}.$

With these formulas, we obtain the vectorized algorithm given in Algorithm 134.

Algorithm 134

Figure 126: Vectorized algorithm for computation of basis functions gradients in 2d

Vectorized elementary matrices (used by version `OptV2`)¶

Element Mass Matrix ¶

We have

$\mathbb{M}^e(T_k) =\frac{|T_k|}{12}\begin{pmatrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{pmatrix}$

Then with $\MAT{K}_g$ definition (see Section New Optimized assembling algorithm (version OptV2)) , we obtain

$\MAT{K}_g(3(i-1)+j,k)=|T_k|\frac{1+\delta_{i,j}}{12} \quad 1\le i,j \le 3,$

We represent in figure 135 the corresponding row-wise operations.

Figure 127: Build of $\mathbb{K}_g$ associated to 2d Mass matrix in

So the vectorized algorithm for $\mathbb{K}_g$ computation is simple and given in Algorithm 136.

Algorithm 136

Figure 128: Vectorized algorithm for $\mathbb{K}_g$ associated to 2d Mass matrix

Note

pyOptFEM.FEM2D.elemMatrixVec.ElemMassMat2DP1Vec(areas)[source]

Compute all the elementaries Mass matrices, $\mathbb{M}^e(T_k)$ for $k\in\{0,\hdots,\nme-1\}$

Parameters:	areas ( $\nme$ numpy array of floats) – areas of all the mesh elements.
Returns:	a one dimensional numpy array of size $9 \nme$

Element Stiffness Matrix ¶

We have $\forall (\il,\jl)\in\ENS{1}{3}^2$

$\mathbb{S}_{\il,\jl}^e(T_k)= |T_k| \DOT{\GRAD\BasisFunc_{\jl}^k}{\GRAD\BasisFunc_{\il}^k}.$

Using vectorized algorithm function $\FNametxt{GradientVec2D}$ given in Algorithm 134, we obtain the vectorized algorithm 137 for $\mathbb{K}_g$ computation of the Stiffness matrix in 2d.

Algorithm 137

Figure 129: Vectorized algorithm for $\mathbb{K}_g$ associated to 2d Stiffness matrix

Note

pyOptFEM.FEM2D.elemMatrixVec.ElemStiffMat2DP1Vec(nme, q, me, areas)[source]

Compute all the elementaries Stiff matrices, $\mathbb{S}^e(T_k)$ for $k\in\{0,\hdots,\nme-1\}$

Parameters:	nme (int) – number of mesh elements, q ( $\nq\times 2$ numpy array of floats) – mesh vertices, me ( $\nme\times 3$ numpy array of integers) – mesh connectivity, areas ( $\nme$ numpy array of floats) – areas of all the mesh elements.
Returns:	a one dimensional numpy array of size $9 \nme$

Element Stiffness Elasticity Matrix ¶

We define on $T_k$ the local alternate basis $\mathcal{B}_a^k$ by

$\mathcal{B}_a^k=\{\BasisFuncTwoD_1^k,\hdots,\BasisFuncTwoD_6^k\}=\left\{ \begin{pmatrix} \BasisFunc_1^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_1^k \end{pmatrix}, \begin{pmatrix} \BasisFunc_2^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_2^k\end{pmatrix}, \begin{pmatrix} \BasisFunc_3^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_3^k \end{pmatrix} \right\}$

where $\BasisFunc_\il^k=\BasisFunc_{\me(\il,k)}.$ With notations of Presentation, we have $\forall (\il,\jl) \in \ENS{1}{6}^2$

$\StiffElasElem_{\il,\jl}(T_k)= \int_{T_k} \Odv^t(\BasisFuncTwoD^k_\jl) \mathbb{C}\Odv(\BasisFuncTwoD^k_\il)d\q= \int_{T_k} \mathcal{H}(\BasisFuncTwoD^k_\jl,\BasisFuncTwoD^k_\il)(\q)d\q$

with, $\forall \vecb{u}=(u_1,u_2)\in\HUnD{\DOM}^2,$ $\forall \vecb{v}=(v_1,v_2)\in\HUnD{\DOM}^2,$

(2) $\begin{array}{lcl} \mathcal{H}(\vecb{u},\vecb{v})&=&{\footnotesize \DOT{\begin{pmatrix} \gamma & 0\\ 0 & \mu \end{pmatrix}\GRAD u_1 }{\GRAD v_1} +\DOT{\begin{pmatrix} 0 & \lambda\\ \mu & 0 \end{pmatrix}\GRAD u_2 }{\GRAD v_1}}\\ &+&{\footnotesize\DOT{\begin{pmatrix} 0 & \mu\\ \lambda & 0 \end{pmatrix}\GRAD u_1 }{\GRAD v_2}+ \DOT{\begin{pmatrix} \mu & 0\\ 0 &\gamma\end{pmatrix}\GRAD u_2 }{\GRAD v_2}} \end{array}$

For example, we can compute explicitely the first two terms in the first column of $\StiffElasElem(T_k)$ which are given by

$\begin{array}{lcl} \StiffElasElem_{1,1}(T_k)&=&\int_{T_k} \mathcal{H}(\BasisFuncTwoD^k_{1},\BasisFuncTwoD^k_{1})(\q)d\q\\ &=&\int_{T_k} \mathcal{H}\left( \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\begin{pmatrix} \gamma & 0\\ 0 & \mu \end{pmatrix}\GRAD \BasisFunc^k_{1} }{\GRAD \BasisFunc^k_{1}} =|T_k|\left(\gamma\DP{\BasisFunc^k_{1}}{x}\DP{\BasisFunc^k_{1}}{x}+\mu\DP{\BasisFunc^k_{1}}{y}\DP{\BasisFunc^k_{1}}{y}\right). \end{array}$

and

$\begin{array}{lcl} \StiffElasElem_{2,1}(T_k)&=&\int_{T_k} \mathcal{H}(\BasisFuncTwoD^k_{1},\BasisFuncTwoD^k_{2})(\q)d\q\\ &=&\int_{T_k} \mathcal{H}\left( \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix}, \begin{pmatrix} 0\\ \BasisFunc^k_{1} \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\begin{pmatrix} 0 & \mu\\ \lambda & 0 \end{pmatrix}\GRAD \BasisFunc^k_{1} }{\GRAD \BasisFunc^k_{1}} =|T_k|(\lambda+\mu)\DP{\BasisFunc^k_{1}}{x}\DP{\BasisFunc^k_{1}}{y}. \end{array}$

Using vectorized algorithm function $\FNametxt{GradientVec2D}$ given in Algorithm 134, we obtain the vectorized algorithm 137 for $\mathbb{K}_g$ computation of the Elasticity Stiffness matrix in 2d.

Algorithm 138

Figure 130: Vectorized algorithm for $\mathbb{K}_g$ associated to 2d Elasticity Stiffness matrix

Note

pyOptFEM.FEM2D.elemMatrixVec.ElemStiffElasMatBaVec2DP1(nme, q, me, areas, L, M)[source]

Compute all the elementaries Stiffness elasticity matrices, $\mathbb{K}^e(T_k)$ for $k\in\{0,\hdots,\nme-1\}$ in local alternate basis.

Parameters:	nme (int) – number of mesh elements, q (`(2,nq)` numpy array of floats) – mesh vertices, me (`(3,nme)` numpy array of integers) – mesh connectivity, areas (`(nme,)` numpy array of floats) – areas of all the mesh elements. L – the $\lambda$ Lame parameter, L – the $\mu$ Lame parameter.
Returns:	a `(36nme,)` numpy* array of floats.

We define on $T_k$ the local block basis $\mathcal{B}_b^k$ by

$\mathcal{B}_b^k=\{\BasisFuncTwoDB_1^k,\hdots,\BasisFuncTwoDB_6^k\}=\left\{ \begin{pmatrix} \BasisFunc_1^k \\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc_2^k \\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc_3^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_1^k \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_2^k\end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_3^k \end{pmatrix} \right\}$

where $\BasisFunc_\il^k=\BasisFunc_{\me(\il,k)}.$

For example, using formula (2), we can explicitly compute the first two terms in the first column of $\StiffElasElem(T_k)$ which are given by

$\begin{array}{lcl} \StiffElasElem_{1,1}(T_k)&=&\int_{T_k} \mathcal{H}(\BasisFuncTwoDB^k_{1},\BasisFuncTwoDB^k_{1})(\q)d\q\\ &=&\int_{T_k} \mathcal{H}\left( \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\begin{pmatrix} \gamma & 0\\ 0 & \mu \end{pmatrix}\GRAD \BasisFunc^k_{1} }{\GRAD \BasisFunc^k_{1}} =|T_k|\left(\gamma\DP{\BasisFunc^k_{1}}{x}\DP{\BasisFunc^k_{1}}{x}+\mu\DP{\BasisFunc^k_{1}}{y}\DP{\BasisFunc^k_{1}}{y}\right). \end{array}$

and

$\begin{array}{lcl} \StiffElasElem_{2,1}(T_k)&=&\int_{T_k} \mathcal{H}(\BasisFuncTwoDB^k_{1},\BasisFuncTwoDB^k_{2})(\q)d\q\\ &=&\int_{T_k} \mathcal{H}\left( \begin{pmatrix} \BasisFunc^k_{1}\\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc^k_{2}\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\begin{pmatrix} \gamma & 0\\ 0 & \mu \end{pmatrix}\GRAD \BasisFunc^k_{1} }{\GRAD \BasisFunc^k_{2}} =|T_k|\left(\gamma\DP{\BasisFunc^k_{1}}{x}\DP{\BasisFunc^k_{2}}{x}+\mu\DP{\BasisFunc^k_{1}}{y}\DP{\BasisFunc^k_{2}}{y} \right). \end{array}$

Using vectorized algorithm function $\FNametxt{GradientVec2D}$ given in Algorithm 134, we obtain the vectorized algorithm 139 for $\mathbb{K}_g$ computation of the Elasticity Stiffness matrix in 2d.

Algorithm 139

Figure 131: Vectorized algorithm for $\mathbb{K}_g$ associated to 2d Elasticity Stiffness matrix

Note

pyOptFEM.FEM2D.elemMatrixVec.ElemStiffElasMatBbVec2DP1(nme, q, me, areas, L, M)[source]

Compute all the elementaries Stiffness elasticity matrices, $\mathbb{K}^e(T_k)$ for $k\in\{0,\hdots,\nme-1\}$ in local block basis.

Parameters:	nme (int) – number of mesh elements, q (`(2,nq)` numpy array of floats) – mesh vertices, me (`(3,nme)` numpy array of integers) – mesh connectivity, areas (`(nme,)` numpy array of floats) – areas of all the mesh elements. L – the $\lambda$ Lame parameter, L – the $\mu$ Lame parameter.
Returns:	a `(36nme,)` numpy* array of floats.

Mesh ¶

class pyOptFEM.FEM2D.mesh.SquareMesh(N, **kwargs)[source]¶

Build meshes of the unit square $[0,1]\times [0,1]$ . Class attributes are :

nq, total number of mesh vertices (points), also denoted $\nq$ .

nme, total number of mesh elements (triangles in 2d),

version, mesh structure version,

q, Numpy array of vertices coordinates, dimension (nq,2) (version 0) or (2,nq) (version 1).

q[j] (version 0) or q[:,j] (version 1) are the two coordinates of the $j$ -th vertex, $j\in\{0,..,nq-1\}$

me, Numpy connectivity array, dimension (nme,3) (version 0) or (3,nme) (version 1).

me[k] (version 0) or me[:,k] (version 1) are the storage index of the three vertices of the $k$ -th triangle in the array q of vertices coordinates, $k\in\{0,...,nme-1\}$ .

areas, Array of mesh elements areas, (nme,) Numpy array.

areas[k] is the area of $k$ -th triangle, k in range(0,nme)

Parameters:	N – points number on each sides of the square

optional parameter : version=0 or version=1

>>> from pyOptFEM.FEM2D import *
>>> Th=SquareMesh(3)
>>> Th.nme,Th.nq
(18, 16)
>>> Th.q
array([[ 0.        ,  0.        ],
       [ 0.33333333,  0.        ],
       [ 0.66666667,  0.        ],
       [ 1.        ,  0.        ],
       [ 0.        ,  0.33333333],
       [ 0.33333333,  0.33333333],
       [ 0.66666667,  0.33333333],
       [ 1.        ,  0.33333333],
       [ 0.        ,  0.66666667],
       [ 0.33333333,  0.66666667],   
       [ 0.66666667,  0.66666667],
       [ 1.        ,  0.66666667],
       [ 0.        ,  1.        ],
       [ 0.33333333,  1.        ],
       [ 0.66666667,  1.        ],
       [ 1.        ,  1.        ]])
>>> PlotMesh(Th)

Figure 132: SquareMesh(3) visualisation

class pyOptFEM.FEM2D.mesh.getMesh(meshfile, **kwargs)[source]¶

Read a FreeFEM++ mesh from file meshfile. Class attributes are :

nq, total number of mesh vertices (points), also denoted $\nq$ .

nme, total number of mesh elements (triangles in 2d),

version, mesh structure version,

q, Numpy array of vertices coordinates, dimension (nq,2) (version 0) or (2,nq) (version 1).

q[j] (version 0) or q[:,j] (version 1) are the two coordinates of the $j$ -th vertex, $j\in\{0,..,nq-1\}$

me, Numpy connectivity array, dimension (nme,3) (version 0) or (3,nme) (version 1).

me[k] (version 0) or me[:,k] (version 1) are the storage index of the three vertices of the $k$ -th triangle in the array q of vertices coordinates, $k\in\{0,...,nme-1\}$ .

areas, Array of mesh elements areas, (nme,) Numpy array.

areas[k] is the area of $k$ -th triangle, k in range(0,nme)

Parameters:	N – points number on each sides of the square

optional parameter : version=0 or version=1

>>> from pyOptFEM.FEM2D import *
>>> Th=getMesh('mesh/disk4-1-5.msh')
>>> PlotMesh(Th)

Figure 133: Visualisation of a FreeFEM++ mesh (disk unit)

FEM2D module¶

Assembly matrices (versions `base`, `OptV1` and `OptV2`)¶

Mass Matrix ¶

Stiffness Matrix ¶

Stiffness Elasticity Matrix ¶

Elementary matrices (used by versions `base` and `OptV1`)¶

Element Mass Matrix ¶

Element Stiffness Matrix ¶

Element Stiffness Elasticity Matrix ¶

Vectorized tools (used by version `OptV2`)¶

Vectorized computation of basis functions gradients ¶

Vectorized elementary matrices (used by version `OptV2`)¶

Element Mass Matrix ¶

Element Stiffness Matrix ¶

Element Stiffness Elasticity Matrix ¶

Mesh ¶

Table Of Contents

Previous topic

Next topic

This Page

Navigation

FEM2D module¶

Previous topic

Next topic

This Page

Quick search

Navigation