Author: | Francois Cuvelier <cuvelier@math.univ-paris13.fr> |
---|---|
Date: | 15/09/2013 |
Contains functions to build some finite element matrices using -Lagrange finite elements on a 2D mesh.
Each assembly matrix is computed by three different versions called base,
OptV1 and OptV2 (see here)
Contents
Let be a triangular mesh of
We denote by
a triangulation of
with the following data structure:
The -Lagrange basis functions associated to
are denoted by
for all
and defined by
We also define the global alternate basis by
and the global block basis by
Assembly of the Mass Matrix by -Lagrange finite elements
using base, OptV1 and OptV2 versions respectively (see report).
The Mass Matrix
is given by
Note
generic syntax:
M = MassAssembling2DP1<version>(nq,nme,me,areas)
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 :
Note
source code
Assembly of the Mass Matrix by -Lagrange finite elements using base version (see report).
Assembly of the Mass Matrix by -Lagrange finite elements using OptV1 version (see report).
Assembly of the Mass Matrix by -Lagrange finite elements using OptV2 version (see report).
Assembly of the Stiffness Matrix by -Lagrange finite elements using base,
OptV1 and OptV2 versions respectively (see report).
The Stiffness Matrix
is given by
Note
generic syntax:
M = StiffAssembling2DP1<version>(nq,nme,q,me,areas)
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
source code
Assembly of the Stiffness Matrix by -Lagrange finite elements using base version (see report).
Assembly of the Stiffness Matrix by -Lagrange finite elements using OptV1 version (see report).
Assembly of the Stiffness Matrix by -Lagrange finite elements using OptV2 version (see report).
Assembly of the Elastic Stiffness Matrix by -Lagrange finite elements using base,
OptV1 and OptV2 versions respectively (see report).
The Elastic Stiffness Matrix
is given by
where
and
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 ,
nme: total number of triangles, also denoted by ,
areas: (nme,) array of areas,
la: the first Lame coefficient in Hooke’s law, denoted by ,
mu: the second Lame coefficient in Hooke’s law, denoted by ,
returns a Scipy CSC sparse matrix of size
where <version> is base, OptV1 or OptV2
Benchmarks of theses functions are presented in Elastic Stiffness 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 (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 2: Sparsity of the Elastic Stiffness 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 (Num=1 or Num=3)
Note
source code
Assembly of the Elasticity Stiffness Matrix by -Lagrange finite elements using OptV2 version (see report).
Assembly of the Elasticity Stiffness Matrix by -Lagrange finite elements using OptV1 version (see report).
Assembly of the Elasticity Stiffness Matrix by -Lagrange finite elements using OptV2 version (see report).
Let be a triangle, of area
and with
,
and
its three vertices. We denote by
,
and
the
-Lagrange local basis functions such that
.
We also define the local alternate basis by
and the local block basis by
The elasticity tensor, , obtained from Hooke’s law with an isotropic material,
defined with the Lamé parameters
and
is given by
and, for a function the strain tensors is given by
The element stiffness matrix, ,
for the
is defined by
We have :
where ,
and
.
Let ,
and
.
The element elastic stiffness matrix, ,
for a given triangle
in the local alternate basis
is defined by
We also have
where is the elasticity tensor and
Note
Returns the element elastic stiffness matrix
for a given triangle
in the local alternate basis
Parameters: |
|
---|---|
Returns: |
|
Type : |
|
The element elastic stiffness matrix, ,
for a given triangle
in the local block basis
is defined by
We also have
where is the elasticity tensor and
Note
Returns the element elastic stiffness matrix
for a given triangle
in the local block basis
Parameters: |
|
---|---|
Returns: |
|
Type : |
|
By construction, the gradients of basis functions are constants on each element
So, we denote,
by
the
array defined,
by
On a triangle
we define
and
Then, we have
With these formulas, we obtain the vectorized algorithm given in Algorithm 12.
Algorithm 12
We have
Then with definition (see Section New Optimized assembly algorithm (OptV2 version)) , we obtain
We represent in figure 13 the corresponding row-wise operations.
So the vectorized algorithm for computation is simple and given in Algorithm 14.
Algorithm 14
Note
Computes all the element Mass matrices for
Parameters: | areas (![]() |
---|---|
Returns: | a one dimensional numpy array of size ![]() |
We have
Using vectorized algorithm function given in Algorithm 12, we obtain
the vectorized algorithm 15 for
computation for the Stiffness matrix in 2d.
Algorithm 15
Note
Computes all the element stiffness matrices for
Parameters: |
|
---|---|
Returns: | a one dimensional numpy array of size |
We define on the local alternate basis
by
where With notations of Presentation, we have
with,
(2)
For example, we can explicitly compute the first two terms in the first column of which are given by
and
Using vectorized algorithm function given in Algorithm 12, we obtain
the vectorized algorithm 15 for
computation for the Elastic Stiffness matrix in 2d.
Algorithm 16
Note
Computes all the element elastic stiffness matrices for
in local alternate basis.
Parameters: | |
---|---|
Returns: | a (36*nme,) numpy array of floats. |
We define on the local block basis
by
where
For example, using formula (2), we can explicitly compute the first two terms in the first column of which are given by
and
Using vectorized algorithm function given in Algorithm 12, we obtain
the vectorized algorithm 17 for
computation for the Elastic Stiffness matrix in 2d.
Algorithm 17
Note
Computes all the element elastic stiffness matrices
for
in local block basis.
Parameters: | |
---|---|
Returns: | a (36*nme,) numpy array of floats. |
Creates meshes of the unit square . Class attributes are :
nq, total number of mesh vertices (points), also denoted
.
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
-th vertex,
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
-th triangle in the array q of vertices coordinates,
.
areas, Array of mesh elements areas, (nme,) Numpy array.
areas[k] is the area of
-th triangle, k in range(0,nme)
Parameters: | N – number of points on each side 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)
Reads a FreeFEM++ mesh from file meshfile. Class attributes are :
nq, total number of mesh vertices (points), also denoted
.
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
-th vertex,
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
-th triangle in the array q of vertices coordinates,
.
areas, Array of mesh elements areas, (nme,) Numpy array.
areas[k] is the area of
-th triangle, k in range(0,nme)
Parameters: | N – number of points on each side of the square |
---|
optional parameter : version=0 or version=1
>>> from pyOptFEM.FEM2D import *
>>> Th=getMesh('mesh/disk4-1-5.msh')
>>> PlotMesh(Th)