CDMATH/tests/examples/SpectrumLaplaceBeltrami3DEF/SpectrumFiniteElements3DLaplace-Beltrami.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Calcul EF du spectre de l'opérateur de Laplace-Beltrami -\triangle sur une surface en 3D
   4 # Author      : Michael Ndjinga
   5 # Copyright   : CEA Saclay 2020
   6 # Description : Utilisation de la méthode des éléménts finis P1 avec champs discrétisés aux noeuds d'un maillage triangulaire
   7 #               Création et sauvegarde des champs résultant en utilisant la librairie CDMATH
   8 #================================================================================================================================
   9
  10 import cdmath
  11 import sys
  12
  13 #Chargement du maillage triangulaire de la surface
  14 #=================================================
  15 my_mesh = cdmath.Mesh(sys.argv[1])
  16 mesh_name=sys.argv[2]
  17 if(not my_mesh.isTriangular()) :
  18         raise ValueError("Wrong cell types : mesh is not made of triangles")
  19 if(my_mesh.getMeshDimension()!=2) :
  20         raise ValueError("Wrong mesh dimension : expected a surface of dimension 2")
  21 if(my_mesh.getSpaceDimension()!=3) :
  22         raise ValueError("Wrong space dimension : expected a space of dimension 3")
  23
  24 nbNodes = my_mesh.getNumberOfNodes()
  25 nbCells = my_mesh.getNumberOfCells()
  26
  27 print("Mesh building/loading done")
  28 print("nb of nodes=", nbNodes)
  29 print("nb of cells=", nbCells)
  30
  31 maxNbNeighbours = my_mesh.getMaxNbNeighbours(cdmath.NODES)+1#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  32
  33 # Construction de la matrice de rigidité
  34 #=======================================
  35 Rigidite=cdmath.SparseMatrixPetsc(nbNodes,nbNodes,maxNbNeighbours)# warning : third argument is number of non zero coefficients per line
  36
  37 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle
  38 GradShapeFunc0=cdmath.Vector(3)
  39 GradShapeFunc1=cdmath.Vector(3)
  40 GradShapeFunc2=cdmath.Vector(3)
  41
  42 normalFace0=cdmath.Vector(3)
  43 normalFace1=cdmath.Vector(3)
  44
  45 nodal_volumes=cdmath.Vector(nbNodes)
  46
  47 #On parcourt les triangles du domaine
  48 for i in range(nbCells):
  49
  50         Ci=my_mesh.getCell(i)
  51
  52         #Contribution à la matrice de rigidité
  53         nodeId0=Ci.getNodeId(0)
  54         nodeId1=Ci.getNodeId(1)
  55         nodeId2=Ci.getNodeId(2)
  56         N0=my_mesh.getNode(nodeId0)
  57         N1=my_mesh.getNode(nodeId1)
  58         N2=my_mesh.getNode(nodeId2)
  59
  60         #Build normal to cell Ci
  61         normalFace0[0]=Ci.getNormalVector(0,0)
  62         normalFace0[1]=Ci.getNormalVector(0,1)
  63         normalFace0[2]=Ci.getNormalVector(0,2)
  64         normalFace1[0]=Ci.getNormalVector(1,0)
  65         normalFace1[1]=Ci.getNormalVector(1,1)
  66         normalFace1[2]=Ci.getNormalVector(1,2)
  67
  68         normalCell = normalFace0.crossProduct(normalFace1)
  69         normalCell = normalCell/normalCell.norm()
  70
  71         cellMat=cdmath.Matrix(4)
  72         cellMat[0,0]=N0.x()
  73         cellMat[0,1]=N0.y()
  74         cellMat[0,2]=N0.z()
  75         cellMat[1,0]=N1.x()
  76         cellMat[1,1]=N1.y()
  77         cellMat[1,2]=N1.z()
  78         cellMat[2,0]=N2.x()
  79         cellMat[2,1]=N2.y()
  80         cellMat[2,2]=N2.z()
  81         cellMat[3,0]=normalCell[0]
  82         cellMat[3,1]=normalCell[1]
  83         cellMat[3,2]=normalCell[2]
  84         cellMat[0,3]=1
  85         cellMat[1,3]=1
  86         cellMat[2,3]=1
  87         cellMat[3,3]=0
  88
  89         #Formule des gradients voir EF P1 -> calcul déterminants
  90         GradShapeFunc0[0]= cellMat.partMatrix(0,0).determinant()/2
  91         GradShapeFunc0[1]=-cellMat.partMatrix(0,1).determinant()/2
  92         GradShapeFunc0[2]= cellMat.partMatrix(0,2).determinant()/2
  93         GradShapeFunc1[0]=-cellMat.partMatrix(1,0).determinant()/2
  94         GradShapeFunc1[1]= cellMat.partMatrix(1,1).determinant()/2
  95         GradShapeFunc1[2]=-cellMat.partMatrix(1,2).determinant()/2
  96         GradShapeFunc2[0]= cellMat.partMatrix(2,0).determinant()/2
  97         GradShapeFunc2[1]=-cellMat.partMatrix(2,1).determinant()/2
  98         GradShapeFunc2[2]= cellMat.partMatrix(2,2).determinant()/2
  99
 100         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
 101         GradShapeFuncs={nodeId0 : GradShapeFunc0}
 102         GradShapeFuncs[nodeId1]=GradShapeFunc1
 103         GradShapeFuncs[nodeId2]=GradShapeFunc2
 104
 105         # Remplissage de  la matrice de rigidité et du second membre
 106         for j in [nodeId0,nodeId1,nodeId2] :
 107                 nodal_volumes[j]+=Ci.getMeasure()/3
 108                 #Contribution de la cellule triangulaire i à la ligne j du système linéaire
 109                 for k in [nodeId0,nodeId1,nodeId2] :
 110                         Rigidite.addValue(j,k,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 111
 112 print("Linear system matrix building done")
 113
 114 # Conditionnement de la matrice de rigidité
 115 #==========================================
 116 cond = Rigidite.getConditionNumber(True)
 117 print("Condition number is ",cond)
 118
 119 # Spectre de la matrice de rigidité
 120 #==================================
 121 #Homogénéisation de la matrice de rigidité (sinon elle tend vers zero de meme que ses valeurs propres)
 122 for i in range(nbNodes):
 123         nodal_volumes[i]=1/nodal_volumes[i]
 124 Rigidite.leftDiagonalScale(nodal_volumes)
 125
 126 nev=9
 127 d=Rigidite.getEigenvectorsDataArrayDouble(nev)
 128 my_eigenfield = cdmath.Field("Eigenvectors field", cdmath.NODES, my_mesh, nev)
 129 my_eigenfield.setFieldByDataArrayDouble(d)
 130
 131 # Sauvegarde du champ résultat
 132 #===========================
 133 my_eigenfield.writeVTK("spectrumFiniteElementsOn"+mesh_name+"LaplaceBeltrami")