CDMATH/tests/examples/SpectrumLaplace2DEF/SpectrumLaplace2DEF_SQUARE.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Calcul EF du spectre de l'opérateur de Laplace 2D -\triangle avec conditions aux limites de Dirichlet u=0
   4 # Author      : Michaël 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 if len(sys.argv) >2 :#load a mesh file
  14     my_mesh = cdmath.Mesh(sys.argv[1])
  15     mesh_name=sys.argv[2]
  16 else :   #rectangular mesh split into triangles
  17     xmin=0
  18     xmax=1
  19     ymin=0
  20     ymax=1
  21
  22     nx=15
  23     ny=15
  24
  25     my_mesh = cdmath.Mesh(xmin,xmax,nx,ymin,ymax,ny,0)
  26     mesh_name="RightTriangles"
  27
  28 if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  29     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  30 if(not my_mesh.isTriangular()) :
  31         raise ValueError("Wrong cell types : mesh is not made of triangles")
  32 eps=1e-6
  33 my_mesh.setGroupAtPlan(0.,0,eps,"DirichletBorder")#Bord GAUCHE
  34 my_mesh.setGroupAtPlan(1.,0,eps,"DirichletBorder")#Bord DROIT
  35 my_mesh.setGroupAtPlan(0.,1,eps,"DirichletBorder")#Bord BAS
  36 my_mesh.setGroupAtPlan(1.,1,eps,"DirichletBorder")#Bord HAUT
  37
  38 nbNodes = my_mesh.getNumberOfNodes()
  39 nbCells = my_mesh.getNumberOfCells()
  40
  41 print("Mesh loading done")
  42 print("Number of nodes=", nbNodes)
  43 print("Number of cells=", nbCells)
  44
  45 #Détermination des noeuds intérieurs
  46 #===================================
  47 nbInteriorNodes = 0
  48 nbBoundaryNodes = 0
  49 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  50 interiorNodes=[]
  51 boundaryNodes=[]
  52
  53 #parcours des noeuds pour discrétisation du second membre et extraction 1) des noeuds intérieur 2) des noeuds frontière 3) du nb max voisins d'un noeud
  54 for i in range(nbNodes):
  55         Ni=my_mesh.getNode(i)
  56         if my_mesh.isBorderNode(i): # Détection des noeuds frontière
  57                 boundaryNodes.append(i)
  58                 nbBoundaryNodes=nbBoundaryNodes+1
  59         else: # Détection des noeuds intérieurs
  60                 interiorNodes.append(i)
  61                 nbInteriorNodes=nbInteriorNodes+1
  62                 maxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only in 2D, otherwise use function Ni.getNumberOfEdges()
  63
  64 print("nb of interior nodes=", nbInteriorNodes)
  65 print("nb of boundary nodes=", nbBoundaryNodes)
  66 print("Max nb of neighbours=", maxNbNeighbours)
  67
  68 # Construction de la matrice de rigidité
  69 #========================================
  70 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
  71
  72 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle (hypothèse 2D)
  73 GradShapeFunc0=cdmath.Vector(2)
  74 GradShapeFunc1=cdmath.Vector(2)
  75 GradShapeFunc2=cdmath.Vector(2)
  76
  77 nodal_volumes=cdmath.Vector(nbInteriorNodes)
  78
  79 #On parcourt les triangles du domaine
  80 for i in range(nbCells):
  81
  82         Ci=my_mesh.getCell(i)
  83
  84         #Contribution à la matrice de rigidité
  85         nodeId0=Ci.getNodeId(0)
  86         nodeId1=Ci.getNodeId(1)
  87         nodeId2=Ci.getNodeId(2)
  88
  89         N0=my_mesh.getNode(nodeId0)
  90         N1=my_mesh.getNode(nodeId1)
  91         N2=my_mesh.getNode(nodeId2)
  92
  93         #Formule des gradients voir EF P1 -> calcul déterminants
  94         GradShapeFunc0[0]= (N1.y()-N2.y())/2
  95         GradShapeFunc0[1]=-(N1.x()-N2.x())/2
  96         GradShapeFunc1[0]=-(N0.y()-N2.y())/2
  97         GradShapeFunc1[1]= (N0.x()-N2.x())/2
  98         GradShapeFunc2[0]= (N0.y()-N1.y())/2
  99         GradShapeFunc2[1]=-(N0.x()-N1.x())/2
 100
 101         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
 102         GradShapeFuncs={nodeId0 : GradShapeFunc0}
 103         GradShapeFuncs[nodeId1]=GradShapeFunc1
 104         GradShapeFuncs[nodeId2]=GradShapeFunc2
 105
 106
 107         # Remplissage de  la matrice de rigidité et du second membre
 108         for j in [nodeId0,nodeId1,nodeId2] :
 109                 if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
 110                         j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
 111                         nodal_volumes[j_int]+=Ci.getMeasure()/3
 112                         #Contribution de la cellule triangulaire i à la ligne j_int du système linéaire
 113                         for k in [nodeId0,nodeId1,nodeId2] :
 114                                 if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
 115                                         k_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
 116                                         Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 117                                 #else: si condition limite non nulle au bord, ajouter la contribution du bord au second membre de la cellule j
 118
 119 print("Stiffness matrix construction done")
 120 quit()
 121 # Conditionnement de la matrice de rigidité
 122 #=================================
 123 cond = Rigidite.getConditionNumber()
 124 print("Condition number is ",cond)
 125
 126 # Spectre de la matrice de rigidité
 127 #==================================
 128 #Homogénéisation de la matrice de rigidité (sinon elle tend vers zero de meme que ses valeurs propres)
 129 for i in range(nbInteriorNodes):
 130         nodal_volumes[i]=1/nodal_volumes[i]
 131 Rigidite.leftDiagonalScale(nodal_volumes)
 132
 133 nev=10
 134 d=Rigidite.getEigenvectorsDataArrayDouble(nev)
 135 my_eigenfield = cdmath.Field("Eigenvectors field", cdmath.NODES, my_mesh, nev)
 136 for j in range(nbInteriorNodes):
 137     for k in range(nev):
 138       my_eigenfield[interiorNodes[j],k]=d[j,k];#remplissage des valeurs pour les noeuds intérieurs
 139 for j in range(nbBoundaryNodes):
 140     for k in range(nev):
 141       my_eigenfield[boundaryNodes[j],k]=0;#remplissage des valeurs pour les noeuds frontière (condition limite)
 142 for k in range(nev):
 143     my_eigenfield.setInfoOnComponent(k,d.getInfoOnComponent(k))
 144
 145 # Sauvegarde du champ résultat
 146 #===========================
 147 my_eigenfield.writeVTK("spectrumFiniteElementsOn"+mesh_name+"Laplace")