Navigation

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)

Assembly matrix (versions base, OptV1 and OptV2)

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 with \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 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 = 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,
  • return 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)
_images/Mass3DP1sparsity.png

Figure 140: Sparsity of Mass Matrix generated with command showSparsity(M)

Note

sources 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 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=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

sources code

pyOptFEM.FEM3D.assembly.StiffAssembling3DP1base(nq, nme, q, me, volumes)[source]

Assembly of the Stiff 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 Stiff 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 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} \Odv^t(\psi_m(q)) \Ocv(\psi_l(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 elasticity 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)
    _images/StiffElas3DP1_Num0_sparsity.png

    Figure 141: Sparsity of Stiffness Elasticity 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)
    _images/StiffElas3DP1_Num3_sparsity.png

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

Note

sources code

pyOptFEM.FEM3D.assembly.StiffElasAssembling3DP1base(nq, nme, q, me, volumes, la, mu, Num)[source]

Assembly of the Stiffness Elasticity 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 Stiffness Elasticity 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 Stiffness Elasticity Matrix by P_1-Lagrange finite elements using OptV2 version (see report).

Elementary matrix (used by versions base and OptV1)

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}\DP{u_1}{x}, & \DP{u_2}{y},& \DP{u_3}{z}, & \DP{u_2}{x}+\DP{u_1}{y},& \DP{u_2}{z}+\DP{u_3}{y},& \DP{u_1}{z}+\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]

Compute 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

sources code

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

Compute 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 Stiffness Elasticity Matrix

  • The element Stiffness Elasticity 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

    sources code

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

    Return the element Stiffness Elasticity 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 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} 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

    sources code

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

    Return the element Stiffness Elasticity 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 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}{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 148.

Algorithm 148

_images/GradientVec3D_algo.png

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

Vectorized elementary matrix (used by version OptV2)

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 assembling algorithm (version OptV2)) , 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 149.

Algorithm 149

_images/ElemMassMat3DP1_algo.png

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

Note

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

Compute all the elementaries 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 148, we obtain the vectorized algorithm 150 for \mathbb{K}_g computation of the Stiffness matrix in 3d.

Algorithm 150

_images/ElemStiffMat3DP1_algo.png

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

Note

pyOptFEM.FEM3D.elemMatrixVec.ElemStiffMat3DP1Vec(nme, q, me, volumes)[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 (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 Stiffness Elasticity Matrix

  • We define on tetrahedra 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 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( \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 148, we obtain the vectorized algorithm 151 for \mathbb{K}_g computation of the Elasticity Stiffness matrix in 3d.

    Algorithm 151

    _images/ElemStiffElasMatBa3DP1_algo.png

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

    Note

    pyOptFEM.FEM3D.elemMatrixVec.ElemStiffElasMatBa3DP1Vec(nme, q, me, volumes, la, mu)[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 ((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 (144*nme,) 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 148, we obtain the vectorized algorithm 152 for \mathbb{K}_g computation of the Elasticity Stiffness matrix in 3d.

    Algorithm 152

    _images/ElemStiffElasMatBb3DP1_algo.png

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

    Note

    pyOptFEM.FEM3D.elemMatrixVec.ElemStiffElasMatBb3DP1Vec(nme, q, me, volumes, 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 ((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 (144*nme,) numpy array of floats.

Mesh

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

Build 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 – points number on each edges of the cube

optional parameter : version=0 or version=1

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

Read 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

Table Of Contents

Previous topic

FEM2D module

Next topic

valid2D module

This Page

Navigation

© Copyright 2013, F. Cuvelier. Created using Sphinx 1.2b1.