FEM3D 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 3D mesh. Each assembly matrix is computed by three differents versions called base, OptV1 and OptV2 (see here)

Contents

Assembly matrix (base, OptV1 and OptV2 versions)
Element matrix (used by base and OptV1 versions)
Vectorized tools (used by OptV2 version)
- Vectorized computation of basis functions gradients
Vectorized element matrix (used by OptV2 version)
Mesh

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

Let $\Omega_h$ be a tetrahedral mesh of $\Omega$ corresponding to 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 & $3 \times {\nq}$ & \begin{minipage}[t]{7.9cm} array of vertices coordinates. $\q(\al,j)$ is the $\nu$-th coordinate of the $j$-th vertex, $\al\in\{1,2,3\}$, $j\in\{1,\hdots,\nq\}.$ The $j$-th vertex will be also denoted by $\q^j$ \end{minipage}& \begin{minipage}[t]{3cm} \texttt{q} (transposed)\\ \texttt{q[j-1]} = $\q^j$ \end{minipage}\\ $\me$ & integer & $4 \times \nme$ & \begin{minipage}[t]{7.9cm} connectivity array. ${\me}(\jl,k)$ is the storage index of the $\beta$-th vertex of the $k$-th element, in the array~$q$, for $\jl\in\{1,\hdots,4\}$ and $k\in\{1,\hdots,{\nme}\}$ \end{minipage}&\texttt{me} (transposed)\\ $\rm volumes$ & double & $1\times {\nme}$ & \begin{minipage}[t]{7.9cm} array of volumes. ${\rm volumes}(k)$ is the $k$-th tetrahedron volume, $k\in\{1,\hdots,\nme\}$ \end{minipage}&\texttt{volumes}\\ \hline \end{tabular}}$

The $P_1$ -Lagrange basis functions associated to $\Omega_h$ are denoted by $\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=\{\FoncBaseDeuxD_1,\hdots,\FoncBaseDeuxD_{3\nq}\}=\left\{ \begin{pmatrix} \varphi_1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \varphi_1 \end{pmatrix}, \hdots, \begin{pmatrix} \varphi_\nq \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_\nq \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \varphi_\nq \end{pmatrix} \right\}$

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

$\mathcal{B}_b=\{\FoncBaseDeuxDB_1,\hdots,\FoncBaseDeuxDB_{3\nq}\}=\left\{ \begin{pmatrix} \varphi_1 \\ 0 \\ 0\end{pmatrix},\hdots, \begin{pmatrix} \varphi_\nq \\ 0 \\0\end{pmatrix}, \begin{pmatrix} 0 \\ \varphi_1 \\ 0\end{pmatrix},\hdots, \begin{pmatrix} 0 \\ \varphi_\nq \\0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0\\ \varphi_1 \end{pmatrix},\hdots, \begin{pmatrix} 0 \\ 0\\ \varphi_\nq \end{pmatrix} \right\}.$

Mass Matrix ¶

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

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

Note

generic syntax:

M = MassAssembling3DP1<version>(nq,nme,me,volumes)

nq: total number of nodes of the mesh, also denoted by $\nq$ ,
nme: total number of tetrahedra, also denoted by $\nme$ ,
me: Connectivity array, (nme,4) array,
volumes: Array of tetrahedra volumes, (nme,) array,
returns a Scipy CSC sparse matrix of size $\nq \times \nq$

where <version> is base, OptV1 or OptV2

>>> from pyOptFEM.FEM3D import *
>>> Th=CubeMesh(5)
>>> Mbase =  MassAssembling3DP1base(Th.nq,Th.nme,Th.me,Th.volumes)
>>> MOptV1= MassAssembling3DP1OptV1(Th.nq,Th.nme,Th.me,Th.volumes)
>>> print(" NormInf(Mbase-MOptV1)=%e " % NormInf(Mbase-MOptV1))
NormInf(Mbase-MOptV1)=1.734723e-18
>>> MOptV2= MassAssembling3DP1OptV2(Th.nq,Th.nme,Th.me,Th.volumes)
>>> print(" NormInf(Mbase-MOptV2)=%e " % NormInf(Mbase-MOptV2))
NormInf(Mbase-MOptV2)=1.734723e-18

We can show sparsity of the Mass matrix :

>>> from pyOptFEM.FEM3D import *
>>> Th=CubeMesh(5)
>>> M=MassAssembling3DP1OptV2(Th.nq,Th.nme,Th.me,Th.areas)
>>> showSparsity(M)
Figure 18: Sparsity of Mass Matrix generated with command showSparsity(M)

Note

source code

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

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

pyOptFEM.FEM3D.assembly.MassAssembling3DP1OptV2(nq, nme, me, volumes)[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 base, OptV1 and OptV2 versions respectively (see report). The Stiffness Matrix $\mathbb{S}$ is given by

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

Note

generic syntax

M=StiffAssembling3DP1<version>(nq,nme,q,me,volumes)

Compute the stiffness sparse matrix where <version> is base, OptV1 or OptV2

Parameters:

nq – total number of nodes of the mesh, also denoted by $\nq$ ,
nme – total number of tetrahedra, also denoted by $\nme$ ,
q (numpy array of float) –
array of vertices coordinates,
- (nq,3) array for base and OptV1 versions,
- (3,nq) array for OptV2 version,
me (numpy array of int) –
Connectivity array,
- (nme,4) array for base and OptV1 versions,
- (4,nme) array for OptV2 version,
volumes (numpy array of floats,) – (nme,) array of tetrahedra volumes,

Returns:

a Scipy CSC sparse matrix of size $2\nq \times 2\nq$

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

>>> from pyOptFEM.FEM3D import *
>>> Th=CubeMesh(5)
>>> Sbase =  StiffAssembling3DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.volumes)
>>> SOptV1= StiffAssembling3DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.volumes)
>>> print(" NormInf(Sbase-SOptV1)=%e " % NormInf(Sbase-SOptV1))
NormInf(Sbase-SOptV1)=2.220446e-15
>>> SOptV2= StiffAssembling3DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.volumes)
>>> print(" NormInf(Sbase-SOptV2)=%e " % NormInf(Sbase-SOptV2))
NormInf(Sbase-SOptV2)=2.220446e-15

Note

source code

pyOptFEM.FEM3D.assembly.StiffAssembling3DP1base(nq, nme, q, me, volumes)[source]: Assembly of the Stiffness Matrix by $P_1$ -Lagrange finite elements using base version (see report).

pyOptFEM.FEM3D.assembly.StiffAssembling3DP1OptV1(nq, nme, q, me, volumes)[source]: Assembly of the Stiffness Matrix by $P_1$ -Lagrange finite elements using OptV1 version (see report).

pyOptFEM.FEM3D.assembly.StiffAssembling3DP1OptV2(nq, nme, q, me, volumes)[source]: Assembly of the Stiffness Matrix by $P_1$ -Lagrange finite elements using OptV2 version (see report).

Elastic Stiffness Matrix ¶

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

$\mathbb{K}_{m,l}=\int_{\Omega_h} \Odv^t(\psi_l(q)) \Ocv(\psi_m(q)), \ \ \forall (m,l) \in \{1,...,3\nq\}^2.$

where $\Ocv=(\Occ_{xx},\Occ_{yy},\Occ_{zz},\Occ_{xy},\Occ_{yz},\Occ_{xz})^t$ and $\Odv=(\Odc_{xx},\Odc_{yy},\Odc_{zz},2\Odc_{xy},2\Odc_{yz},2\Odc_{xz})^t$ are the elastic stress and strain tensors respectively.

Note

generic syntax

M=StiffElasAssembling3DP1<version>(nq,nme,q,me,volumes,la,mu,Num)

Compute the elastic stiffness sparse matrix where <version> is base, OptV1 or OptV2

Parameters:

nq – total number of nodes of the mesh,
nme – total number of tetrahedrons,
q (numpy array of float) –
vertices coordinates,
- (nq,3) array for base and OptV1 versions,
- (3,nq) array for OptV2 version,
me (numpy array of int) –
connectivity array,
- (nme,4) array for base and OptV1 versions,
- (4,nme) array for OptV2 version,
volumes (numpy array of floats,) – (nme,) array of tetrahedra volumes,
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.

Returns:

a Scipy CSC sparse matrix of size $3nq\times 3nq$

>>> from pyOptFEM.FEM3D import *
>>> Th=CubeMesh(5)
>>> Kbase =  StiffElasAssembling3DP1base(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,2,0.5,0)
>>> KOptV1= StiffElasAssembling3DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,2,0.5,0)
>>> print(" NormInf(Kbase-KOptV1)=%e " % NormInf(Kbase-KOptV1))
NormInf(Kbase-KOptV1)=1.332268e-15
>>> KOptV2= StiffElasAssembling3DP1OptV2(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,2,0.5,0)
>>> print(" NormInf(Kbase-KOptV2)=%e " % NormInf(Kbase-KOptV2))
NormInf(Kbase-KOptV2)=1.332268e-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.FEM3D import *
>>> Th=CubeMesh(5)
>>> K0=StiffElasAssembling3DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,2,0.5,0)
>>> showSparsity(K0)

Figure 19: Sparsity of the Elastic Stiffness Matrix generated with global alternate numbering (Num=0 or 2)

global block basis $\mathcal{B}_a$ (Num=1 or Num=3)
```
>>> K3=StiffElasAssembling3DP1OptV1(Th.nq,Th.nme,Th.q,Th.me,Th.volumes,2,0.5,3)
>>> showSparsity(K3)
```
Figure 20: Sparsity of the Elastic Stiffness Matrix generated with global block numbering (Num=1 or 3)

Note

source code

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

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

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

Element matrix (used by `base` and `OptV1` versions)¶

Let $T$ be a tetrahedron of volume $|T|$ and with $\tilde{\q}^1$ , $\tilde{\q}^2$ , $\tilde{\q}^3$ and $\tilde{\q}^4$ its four vertices. We denote by $\tilde{\varphi}_1$ , $\tilde{\varphi}_2$ , $\tilde{\varphi}_3$ and $\tilde{\varphi}_4$ the $P_1$ -Lagrange local basis functions such that $\tilde{\varphi}_i(\pmb{q}^j)=\delta_{i,j},\ \forall(i,j)\in\{1,\hdots,4\}^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}}_{12}\}=\left\{ \begin{pmatrix} \tilde{\varphi}_1 \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_1 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_1 \end{pmatrix}, \hdots, \begin{pmatrix} \tilde{\varphi}_4 \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \tilde{\varphi}_4 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_4 \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}}_{12}\}=\left\{ \begin{pmatrix} \tilde{\varphi}_1 \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_2 \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_3 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \tilde{\varphi}_4 \\ 0 \\ 0\end{pmatrix},\hdots, \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_1 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_2\end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_3 \end{pmatrix} \begin{pmatrix} 0 \\ 0 \\ \tilde{\varphi}_4 \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 & \lambda & 0 & 0 & 0\\ \lambda & \lambda+2\mu & \lambda & 0 & 0 & 0\\ \lambda & \lambda & \lambda+2\mu & 0 & 0 & 0\\ 0 & 0 & 0 & \mu & 0 & 0\\ 0 & 0 & 0 & 0 & \mu & 0\\ 0 & 0 & 0 & 0 & 0 & \mu \end{pmatrix}$

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

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

Element Mass Matrix ¶

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

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

We obtain :

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

Note

source code

pyOptFEM.FEM3D.elemMatrix.ElemMassMat3DP1(V)[source]

Computes the element mass matrix $\mathbb{M}^e(T)$ for a given tetrahedron $T$ of volume $|T|$

Parameters:	V (float) – volume of the tetrahedron.
Returns:	$4 \times 4$ 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,4\}^2$

Note

source code

pyOptFEM.FEM3D.elemMatrix.ElemStiffMat3DP1(ql, volume)[source]

Computes the element stiffness matrix $\mathbb{S}^e(T)$ for a given tetrahedron $T$

Parameters:	ql ( $2 \times 4$ numpy array) – the four vertices of the tetrahedron, volume (float) – volume of the tetrahedron.
Returns:
Type :	$4 \times 4$ numpy array of floats.

Element Elastic Stiffness Matrix ¶

The element elastic stiffness matrix, $\mathbb{K}^e(T)$ , for a given tetrahedron $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,12\}^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}$ is a $6 \times 12$ matrix defined by

$\mathbb{B}=\begin{pmatrix} \Odv(\tilde{\pmb{\psi}}_1)&| & \Odv(\tilde{\pmb{\psi}}_2)& \hdots & \Odv(\tilde{\pmb{\psi}}_{11})&| & \Odv(\tilde{\pmb{\psi}}_{12}) \end{pmatrix}$

So in $\tilde{\mathcal{B}}_a$ basis we obtain

$\mathbb{B}=\begin{pmatrix} \DP{\tilde{\varphi}_1}{x}& 0 &0 & & \DP{\tilde{\varphi}_4}{x}& 0 &0 \\ 0& \DP{\tilde{\varphi}_1}{y} & 0& & 0& \DP{\tilde{\varphi}_4}{y} & 0\\ 0& 0 & \DP{\tilde{\varphi}_1}{z} & & 0& 0 & \DP{\tilde{\varphi}_4}{z}\\ \DP{\tilde{\varphi}_1}{y}& \DP{\tilde{\varphi}_1}{x} & 0 & \hdots\hdots& \DP{\tilde{\varphi}_4}{y}& \DP{\tilde{\varphi}_4}{x} & 0\\ 0&\DP{\tilde{\varphi}_1}{z}& \DP{\tilde{\varphi}_1}{y} & & 0&\DP{\tilde{\varphi}_4}{z}& \DP{\tilde{\varphi}_4}{y}\\ \DP{\tilde{\varphi}_1}{z}& 0 & \DP{\tilde{\varphi}_1}{x} & & \DP{\tilde{\varphi}_4}{z}& 0 & \DP{\tilde{\varphi}_4}{x} \end{pmatrix}$

Note

source code

pyOptFEM.FEM3D.elemMatrix.ElemStiffElasMatBa3DP1(ql, V, C)[source]¶

Returns the element elastic stiffness matrix $\mathbb{K}^e(T)$ for a given tetrahedron $T$ in the local alternate basis $\mathcal{B}_a$

Parameters:	ql ( $4 \times 2$ numpy array) – contains the four vertices of the tetrahedron, V (float) – volume of the tetrahedron H ( $6 \times 6$ numpy array) – Elasticity tensor, $\mathbb{H}$ .
Returns:	$\mathbb{K}^e(T)$ in $\mathcal{B}_a$ basis.
Type :	$12 \times 12$ numpy array of floats.

The element elastic stiffness 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}$ is a $6 \times 12$ matrix defined by

$\mathbb{B}=\begin{pmatrix} \Odv(\tilde{\pmb{\phi}}_1)&\vdots & \Odv(\tilde{\pmb{\phi}}_2)& \hdots & \Odv(\tilde{\pmb{\phi}}_{11})&\vdots & \Odv(\tilde{\pmb{\phi}}_{12}) \end{pmatrix}$

So in $\tilde{\mathcal{B}}_b$ basis we obtain

$\mathbb{B}=\begin{pmatrix} \DP{\tilde{\varphi}_1}{x}& \hdots &\DP{\tilde{\varphi}_4}{x} & 0& \hdots & 0 &0& \hdots & 0 \\ 0& \hdots & 0 & \DP{\tilde{\varphi}_1}{y}& \hdots & \DP{\tilde{\varphi}_4}{y}& 0 & \hdots & 0\\ 0& \hdots & 0 &0 & \hdots & 0 & \DP{\tilde{\varphi}_1}{z} & \hdots &\DP{\tilde{\varphi}_4}{z}\\ \DP{\tilde{\varphi}_1}{y}& \hdots & \DP{\tilde{\varphi}_4}{y}& \DP{\tilde{\varphi}_1}{x}& \hdots & \DP{\tilde{\varphi}_4}{x}& 0& \hdots & 0\\ 0& \hdots & 0 & \DP{\tilde{\varphi}_1}{z}& \hdots & \DP{\tilde{\varphi}_4}{z}& \DP{\tilde{\varphi}_1}{y}& \hdots &\DP{\tilde{\varphi}_4}{y}\\ \DP{\tilde{\varphi}_1}{z}& \hdots & \DP{\tilde{\varphi}_4}{z} & 0& \hdots & 0 & \DP{\tilde{\varphi}_1}{x}& \hdots & \DP{\tilde{\varphi}_4}{x} \end{pmatrix}$

Note

source code

pyOptFEM.FEM3D.elemMatrix.ElemStiffElasMatBb3DP1(ql, V, C)[source]¶

Returns the element elastic stiffness matrix $\mathbb{K}^e(T)$ for a given tetrahedron $T$ in the local block basis $\mathcal{B}_b$

Parameters:	ql ( $4 \times 2$ numpy array) – contains the four vertices of the tetrahedron, V (float) – volume of the tetrahedron H ( $6 \times 6$ numpy array) – Elasticity tensor, $\mathbb{H}$ .
Returns:	$\mathbb{K}^e(T)$ in $\mathcal{B}_b$ basis.
Type :	$12 \times 12$ numpy array of floats.

Vectorized tools (used by `OptV2` version)¶

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}{4},$ by $\vecb{G}^\il$ the $3 \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 $T_k$ tetrahedra we set

$\begin{array}{lcllcl} \vecb{D}^{12}&=&\q^{\me(1,k)}-\q^{\me(2,k)},\ & \vecb{D}^{23}&=&\q^{\me(2,k)}-\q^{\me(3,k)}\\ \vecb{D}^{13}&=&\q^{\me(1,k)}-\q^{\me(3,k)},\ & \vecb{D}^{24}&=&\q^{\me(2,k)}-\q^{\me(4,k)}\\ \vecb{D}^{14}&=&\q^{\me(1,k)}-\q^{\me(4,k)},\ & \vecb{D}^{34}&=&\q^{\me(3,k)}-\q^{\me(4,k)} \end{array}$

Then, we have

$\begin{array}{ll} \GRAD\BasisFunc_{1}^k(\q)=\frac{1}{6|T_k|} \begin{pmatrix} -\vecb{D}^{23}_y \vecb{D}^{24}_z + \vecb{D}^{23}_z \vecb{D}^{24}_y\\ \vecb{D}^{23}_x \vecb{D}^{24}_z - \vecb{D}^{23}_z \vecb{D}^{24}_x\\ -\vecb{D}^{23}_x \vecb{D}^{24}_y + \vecb{D}^{23}_y \vecb{D}^{24}_x \end{pmatrix},&\GRAD\BasisFunc_{2}^k(\q)=\frac{1}{6|T_k|} \begin{pmatrix} \vecb{D}^{13}_y \vecb{D}^{14}_z - \vecb{D}^{13}_z \vecb{D}^{14}_y\\ -\vecb{D}^{13}_x \vecb{D}^{14}_z + \vecb{D}^{13}_z \vecb{D}^{14}_x\\ \vecb{D}^{13}_x \vecb{D}^{14}_y - \vecb{D}^{13}_y \vecb{D}^{14}_x \end{pmatrix}\\ \GRAD\BasisFunc_{3}^k(\q)=\frac{1}{6|T_k|} \begin{pmatrix} -\vecb{D}^{12}_y \vecb{D}^{14}_z + \vecb{D}^{12}_z \vecb{D}^{14}_y\\ \vecb{D}^{12}_x \vecb{D}^{14}_z - \vecb{D}^{12}_z \vecb{D}^{14}_x\\ -\vecb{D}^{12}_x \vecb{D}^{14}_y + \vecb{D}^{12}_y \vecb{D}^{14}_x \end{pmatrix},& \GRAD\BasisFunc_{4}^k(\q)=\frac{1}{6|T_k|} \begin{pmatrix} \vecb{D}^{12}_y \vecb{D}^{13}_z - \vecb{D}^{12}_z \vecb{D}^{13}_y\\ -\vecb{D}^{12}_x \vecb{D}^{13}_z + \vecb{D}^{12}_z \vecb{D}^{13}_x\\ \vecb{D}^{12}_x \vecb{D}^{13}_y - \vecb{D}^{12}_y \vecb{D}^{13}_x \end{pmatrix} \end{array}$

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

Algorithm 26

Figure 21: Vectorized algorithm for computation of basis functions gradients in 3d

Vectorized element matrix (used by `OptV2` version)¶

Element Mass Matrix ¶

We have

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

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

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

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

Algorithm 27

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

Note

pyOptFEM.FEM3D.elemMatrixVec.ElemMassMat3DP1Vec(nme, volumes)[source]

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

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

Element Stiffness Matrix ¶

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

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

Using vectorized algorithm function $\FNametxt{GradientVec3D}$ given in Algorithm 26, we obtain the vectorized algorithm 28 for $\mathbb{K}_g$ computation for the Stiffness matrix in 3d.

Algorithm 28

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

Note

pyOptFEM.FEM3D.elemMatrixVec.ElemStiffMat3DP1Vec(nme, q, me, volumes)[source]

Computes all the element stiffness matrices $\mathbb{S}^e(T_k)$ for $k\in\{0,\hdots,\nme-1\}$

Parameters:	nme (int) – number of mesh elements, q ( $3\times \nq$ numpy array of floats) – mesh vertices, me ( $4 \times\nme$ 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 Elastic Stiffness Matrix ¶

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

$\begin{array}{c} \mathcal{B}_a^k=\{\BasisFuncTwoD_1^k,\hdots,\BasisFuncTwoD_{12}^k\}\\=\\ \left\{\tiny \begin{pmatrix} \BasisFunc_1^k \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_1^k \\0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_1^k \end{pmatrix}, \begin{pmatrix} \BasisFunc_2^k \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_2^k \\0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_2^k \end{pmatrix}, \begin{pmatrix} \BasisFunc_3^k \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_3^k \\0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_3^k \end{pmatrix}, \begin{pmatrix} \BasisFunc_4^k \\ 0 \\ 0\end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_4^k \\0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_4^k \end{pmatrix} \right\} \end{array}$

where $\BasisFunc_\il^k=\BasisFunc_{\me(\il,k)}.$ With notations of Presentation, we have $\forall (\il,\jl) \in \ENS{1}{12}^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,u_3)\in\HUnD{\DOM}^3,$ $\forall \vecb{v}=(v_1,v_2,v_3)\in\HUnD{\DOM}^3,$ by

$\begin{array}{c} \mathcal{H}(\vecb{u},\vecb{v})\\=\\ \tiny{ \DOT{\begin{pmatrix} \gamma & 0 &0\\ 0 & \mu &0\\ 0 & 0 &\mu\end{pmatrix}\GRAD u_1 }{\GRAD v_1} +\DOT{\begin{pmatrix} 0 & \lambda & 0\\ \mu & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\GRAD u_2 }{\GRAD v_1} +\DOT{\begin{pmatrix} 0 & 0 & \lambda\\ 0 & 0 & 0 \\ \mu & 0 & 0 \end{pmatrix}\GRAD u_3 }{\GRAD v_1} }\\ %&+& \tiny{+ \DOT{\begin{pmatrix} 0 & \mu &0\\ \lambda & 0 &0\\ 0 & 0 &0\end{pmatrix}\GRAD u_1 }{\GRAD v_2} +\DOT{\begin{pmatrix} \mu & 0 & 0\\ 0 & \gamma & 0 \\ 0 & 0 & \mu \end{pmatrix}\GRAD u_2 }{\GRAD v_2} +\DOT{\begin{pmatrix} 0 & 0 & 0\\ 0 & 0 & \lambda \\ 0 & \mu & 0 \end{pmatrix}\GRAD u_3 }{\GRAD v_2} }\\ %&+& \tiny{+ \DOT{\begin{pmatrix} 0 & 0 &\mu\\ 0 & 0 &0\\ \lambda & 0 & 0\end{pmatrix}\GRAD u_1 }{\GRAD v_3} +\DOT{\begin{pmatrix} 0 & 0 & 0\\ 0 & 0 & \mu \\ 0 & \lambda & 0 \end{pmatrix}\GRAD u_2 }{\GRAD v_3} +\DOT{\begin{pmatrix} \mu & 0 & 0\\ 0 & \mu & 0 \\ 0 & 0 & \gamma \end{pmatrix}\GRAD u_3 }{\GRAD v_3} } \end{array}$

where $\lambda$ and $\mu$ are the Lame coefficients and $\gamma=\lambda+2\mu.$

For example, 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}(\BasisFuncTwoD^k_{1},\BasisFuncTwoD^k_{1})(\q)d\q\\ &=&\int_{T_k} \mathcal{H}\left( \tiny\begin{pmatrix} \BasisFunc^k_{1}\\ 0\\ 0 \end{pmatrix}, \tiny\begin{pmatrix} \BasisFunc^k_{1}\\ 0\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\tiny\begin{pmatrix} \gamma & 0 &0\\ 0 & \mu &0\\ 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}+\DP{\BasisFunc^k_{1}}{z}\DP{\BasisFunc^k_{1}}{z}) \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( \tiny\begin{pmatrix} \BasisFunc^k_{1}\\ 0\\0 \end{pmatrix}, \begin{pmatrix} 0\\ \BasisFunc^k_{1}\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\tiny\begin{pmatrix} 0 & \mu &0\\ \lambda & 0 &0\\ 0 & 0 &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{GradientVec3D}$ given in Algorithm 26, we obtain the vectorized algorithm 29 for $\mathbb{K}_g$ computation for the Elastic Stiffness matrix in 3d.

Algorithm 29

Figure 24: Vectorized algorithm for $\mathbb{K}_g$ associated to 3d Elastic Stiffness matrix

Note

pyOptFEM.FEM3D.elemMatrixVec.ElemStiffElasMatBa3DP1Vec(nme, q, me, volumes, la, mu)[source]

Computes all the element elastic stiffness 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 (`(3,nq)` numpy array of floats) – mesh vertices, me (`(4,nme)` numpy array of integers) – mesh connectivity, volumes (`(nme,)` numpy array of floats) – volumes of all the mesh elements. la (float) – the $\\lambda$ Lame parameter, mu (float) – the $\\mu$ Lame parameter.
Returns:	a `(144nme,)` numpy* array of floats.

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

$\begin{array}{c} \mathcal{B}_b^k=\{\BasisFuncTwoDB_1^k,\hdots,\BasisFuncTwoDB_{12}^k\} \\ = \\ \left\{\tiny \begin{pmatrix} \BasisFunc_1^k \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc_2^k \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc_3^k \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc_4^k \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_1^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_2^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_3^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ \BasisFunc_4^k \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_1^k \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_2^k \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_3^k \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ \BasisFunc_4^k \end{pmatrix} \right\} \end{array}$

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

For example, using formula (?), 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( \tiny\begin{pmatrix} \BasisFunc^k_{1}\\ 0\\0 \end{pmatrix}, \begin{pmatrix} \BasisFunc^k_{1}\\ 0\\0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\tiny\begin{pmatrix} \gamma & 0 &0\\ 0 & \mu &0\\ 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} +\DP{\BasisFunc^k_{1}}{z}\DP{\BasisFunc^k_{1}}{z})\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( \tiny \begin{pmatrix} \BasisFunc^k_{1}\\ 0\\ 0 \end{pmatrix}, \begin{pmatrix} \BasisFunc^k_{2}\\ 0\\ 0 \end{pmatrix} \right)(\q)d\q\\ &=&|T_k| \DOT{\tiny\begin{pmatrix} \gamma & 0 &0\\ 0 & \mu &0\\ 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}+\DP{\BasisFunc^k_{1}}{z}\DP{\BasisFunc^k_{2}}{z}) \right). \end{array}$

Using vectorized algorithm function $\FNametxt{GradientVec3D}$ given in Algorithm 26, we obtain the vectorized algorithm 30 for $\mathbb{K}_g$ computation for the Elastic Stiffness matrix in 3d.

Algorithm 30

Figure 25: Vectorized algorithm for $\mathbb{K}_g$ associated to 3d Elastic Stiffness matrix

Note

pyOptFEM.FEM3D.elemMatrixVec.ElemStiffElasMatBb3DP1Vec(nme, q, me, volumes, L, M)[source]

Compute all the element elastic stiffness 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 (`(3,nq)` numpy array of floats) – mesh vertices, me (`(4,nme)` numpy array of integers) – mesh connectivity, volumes (`(nme,)` numpy array of floats) – volumes of all the mesh elements. la (float) – the $\\lambda$ Lame parameter, mu (float) – the $\\mu$ Lame parameter.
Returns:	a `(144nme,)` numpy* array of floats.

Mesh ¶

class pyOptFEM.FEM3D.mesh.CubeMesh(N, **kwargs)[source]¶

Creates meshes of the unit cube $[0,1]^3$ . Class attributes are :

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

nme, total number of mesh elements (tetrahedra in 3d),

version, mesh structure version,

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

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

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

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

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

volumes[k] is the volume of $k$ -th tetrahedron, k in range(0,nme)

Parameters:	N – number of points on each edge of the cube

optional parameter : version=0 or version=1

class pyOptFEM.FEM3D.mesh.getMesh(filename, **kwargs)[source]¶

Reads a medit mesh from file meshfile. Class attributes are :

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

nme, total number of mesh elements (tetrahedra in 3d),

version, mesh structure version,

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

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

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

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

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

volumes[k] is the volume of $k$ -th tetrahedron, k in range(0,nme)

Parameters:	meshfile – medit mesh file

optional parameter : version=0 or version=1

FEM3D module¶

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

Mass Matrix ¶

Stiffness Matrix ¶

Elastic Stiffness Matrix ¶

Element matrix (used by `base` and `OptV1` versions)¶

Element Mass Matrix ¶

Element Stiffness Matrix ¶

Element Elastic Stiffness Matrix ¶

Vectorized tools (used by `OptV2` version)¶

Vectorized computation of basis functions gradients ¶

Vectorized element matrix (used by `OptV2` version)¶

Element Mass Matrix ¶

Element Stiffness Matrix ¶

Element Elastic Stiffness Matrix ¶

Mesh ¶

Table Of Contents

Previous topic

Next topic

This Page

Navigation

FEM3D module¶

Previous topic

Next topic

This Page

Quick search

Navigation