CDMATH/tests/examples/Poisson3DTorusEF/FiniteElements3DPoissonTorus.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Laplace-Beltrami -\triangle u = f sur un tore
   4 # Author      : Michael Ndjinga
   5 # Copyright   : CEA Saclay 2018
   6 # Description : Utilisation de la méthode des éléménts finis P1 avec champs u et f discrétisés aux noeuds d'un maillage triangulaire
   7 #               Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
   8 #               Référence : M. A. Olshanskii, A. Reusken, and J. Grande. A finite element method for elliptic equations
   9 #                           on surfaces. SIAM J. Num. Anal., 47, p. 3355
  10 #               Résolution d'un système linéaire à matrice singulière : les vecteurs constants sont dans le noyau
  11 #================================================================================================================================
  12
  13 import cdmath
  14 from math import sin, cos, atan2, sqrt
  15 import PV_routines
  16 import VTK_routines
  17 import paraview.simple as pvs
  18
  19 #Chargement du maillage triangulaire du tore
  20 #=======================================================================================
  21 my_mesh = cdmath.Mesh("meshTorus.med")
  22 if(not my_mesh.isTriangular()) :
  23         raise ValueError("Wrong cell types : mesh is not made of triangles")
  24 if(my_mesh.getMeshDimension()!=2) :
  25         raise ValueError("Wrong mesh dimension : expected a surface of dimension 2")
  26 if(my_mesh.getSpaceDimension()!=3) :
  27         raise ValueError("Wrong space dimension : expected a space of dimension 3")
  28
  29 nbNodes = my_mesh.getNumberOfNodes()
  30 nbCells = my_mesh.getNumberOfCells()
  31
  32 print("Mesh building/loading done")
  33 print("nb of nodes=", nbNodes)
  34 print("nb of cells=", nbCells)
  35
  36 # Torus radii (calculation  will fail if the mesh is not correct)
  37 R=1 #Grand rayon
  38 r=0.6 #Petit rayon
  39
  40 #Discrétisation du second membre, de la solution exacte et détermination des noeuds intérieurs
  41 #======================================================================
  42 my_RHSfield = cdmath.Field("RHS field", cdmath.NODES, my_mesh, 1)
  43 exactSolField = cdmath.Field("Exact solution field", cdmath.NODES, my_mesh, 1)
  44
  45 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  46
  47 #parcours des noeuds pour discrétisation du second membre et extraction du nb max voisins d'un noeud
  48 for i in range(nbNodes):
  49         Ni=my_mesh.getNode(i)
  50         x = Ni.x()
  51         y = Ni.y()
  52         z = Ni.z()
  53
  54         theta=atan2(z,sqrt(x*x+y*y)-R)
  55         phi=atan2(y,x)
  56
  57         exactSolField[i] = sin(3*phi)*cos(3*theta+ phi) # for the exact solution we use the funtion given in the article of Olshanskii, Reusken 2009, page 19
  58         my_RHSfield[i] = 9*sin(3*phi)*cos(3*theta+ phi)/(r*r) + (10*sin(3*phi)*cos(3*theta+ phi) + 6*cos(3*phi)*sin(3*theta+ phi))/((R+r*cos(theta))*(R+r*cos(theta))) - 3*sin(theta)*sin(3*phi)*sin(3*theta+ phi)/(r*(R+r*cos(theta))) #for the right hand side we use the function given in the article of Olshanskii, Reusken 2009, page 19
  59         if my_mesh.isBorderNode(i): # Détection des noeuds frontière
  60                 raise ValueError("Mesh should not contain borders")
  61         else:
  62                 maxNbNeighbours = max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only for planar cells, otherwise use function Ni.getNumberOfEdges()
  63
  64 print("Right hand side discretisation done")
  65 print("Max nb of neighbours=", maxNbNeighbours)
  66 print("Integral of the RHS", my_RHSfield.integral(0))
  67
  68 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
  69 #=======================================================================================
  70 Rigidite=cdmath.SparseMatrixPetsc(nbNodes,nbNodes,maxNbNeighbours)
  71 RHS=cdmath.Vector(nbNodes)
  72
  73 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle
  74 GradShapeFunc0=cdmath.Vector(3)
  75 GradShapeFunc1=cdmath.Vector(3)
  76 GradShapeFunc2=cdmath.Vector(3)
  77
  78 normalFace0=cdmath.Vector(3)
  79 normalFace1=cdmath.Vector(3)
  80
  81 #On parcourt les triangles du domaine
  82 for i in range(nbCells):
  83
  84         Ci=my_mesh.getCell(i)
  85
  86         #Contribution à la matrice de rigidité
  87         nodeId0=Ci.getNodeId(0)
  88         nodeId1=Ci.getNodeId(1)
  89         nodeId2=Ci.getNodeId(2)
  90         N0=my_mesh.getNode(nodeId0)
  91         N1=my_mesh.getNode(nodeId1)
  92         N2=my_mesh.getNode(nodeId2)
  93
  94         #Build normal to cell Ci
  95         normalFace0[0]=Ci.getNormalVector(0,0)
  96         normalFace0[1]=Ci.getNormalVector(0,1)
  97         normalFace0[2]=Ci.getNormalVector(0,2)
  98         normalFace1[0]=Ci.getNormalVector(1,0)
  99         normalFace1[1]=Ci.getNormalVector(1,1)
 100         normalFace1[2]=Ci.getNormalVector(1,2)
 101
 102         normalCell = normalFace0.crossProduct(normalFace1)
 103         normalCell = normalCell/normalCell.norm()
 104
 105         cellMat=cdmath.Matrix(4)
 106         cellMat[0,0]=N0.x()
 107         cellMat[0,1]=N0.y()
 108         cellMat[0,2]=N0.z()
 109         cellMat[1,0]=N1.x()
 110         cellMat[1,1]=N1.y()
 111         cellMat[1,2]=N1.z()
 112         cellMat[2,0]=N2.x()
 113         cellMat[2,1]=N2.y()
 114         cellMat[2,2]=N2.z()
 115         cellMat[3,0]=normalCell[0]
 116         cellMat[3,1]=normalCell[1]
 117         cellMat[3,2]=normalCell[2]
 118         cellMat[0,3]=1
 119         cellMat[1,3]=1
 120         cellMat[2,3]=1
 121         cellMat[3,3]=0
 122
 123         #Formule des gradients voir EF P1 -> calcul déterminants
 124         GradShapeFunc0[0]= cellMat.partMatrix(0,0).determinant()/2
 125         GradShapeFunc0[1]=-cellMat.partMatrix(0,1).determinant()/2
 126         GradShapeFunc0[2]= cellMat.partMatrix(0,2).determinant()/2
 127         GradShapeFunc1[0]=-cellMat.partMatrix(1,0).determinant()/2
 128         GradShapeFunc1[1]= cellMat.partMatrix(1,1).determinant()/2
 129         GradShapeFunc1[2]=-cellMat.partMatrix(1,2).determinant()/2
 130         GradShapeFunc2[0]= cellMat.partMatrix(2,0).determinant()/2
 131         GradShapeFunc2[1]=-cellMat.partMatrix(2,1).determinant()/2
 132         GradShapeFunc2[2]= cellMat.partMatrix(2,2).determinant()/2
 133
 134         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
 135         GradShapeFuncs={nodeId0 : GradShapeFunc0}
 136         GradShapeFuncs[nodeId1]=GradShapeFunc1
 137         GradShapeFuncs[nodeId2]=GradShapeFunc2
 138
 139         # Remplissage de  la matrice de rigidité et du second membre
 140         for j in [nodeId0,nodeId1,nodeId2] :
 141                 #Ajout de la contribution de la cellule triangulaire i au second membre du noeud j
 142                 RHS[j]=Ci.getMeasure()/3*my_RHSfield[j]+RHS[j] # intégrale dans le triangle du produit f x fonction de base
 143                 #Contribution de la cellule triangulaire i à la ligne j du système linéaire
 144                 for k in [nodeId0,nodeId1,nodeId2] :
 145                         Rigidite.addValue(j,k,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 146
 147 print("Linear system matrix building done")
 148
 149 # Conditionnement de la matrice de rigidité
 150 #=================================
 151 cond = Rigidite.getConditionNumber(True)
 152 print("Condition number is ",cond)
 153
 154 # Résolution du système linéaire
 155 #=================================
 156 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"GMRES","ILU")#Remplacer CG par CHOLESKY pour solveur direct
 157 LS.setMatrixIsSingular()#En raison de l'absence de bord
 158 SolSyst=LS.solve()
 159 print( "Preconditioner used : ", LS.getNameOfPc() )
 160 print( "Number of iterations used : ", LS.getNumberOfIter() )
 161 print( "Final residual : ", LS.getResidu() )
 162 print("Linear system solved")
 163
 164 # Création du champ résultat
 165 #===========================
 166 my_ResultField = cdmath.Field("Numerical result field", cdmath.NODES, my_mesh, 1)
 167 for j in range(nbNodes):
 168     my_ResultField[j]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
 169 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 170 my_ResultField.writeVTK("FiniteElementsOnTorusPoisson")
 171
 172 print("Integral of the numerical solution", my_ResultField.integral(0))
 173 print("Numerical solution of Poisson equation on a torus using finite elements done")
 174
 175 #Calcul de l'erreur commise par rapport à la solution exacte
 176 #===========================================================
 177 max_sol_exacte=exactSolField.normMax()[0]
 178 erreur_max=(exactSolField - my_ResultField).normMax()[0]
 179 max_sol_num=my_ResultField.max()
 180 min_sol_num=my_ResultField.min()
 181
 182 print("Relative error =  max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_max/max_sol_exacte)
 183 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 184 print("Maximum exact solution = ", exactSolField.max(), " Minimum exact solution = ", exactSolField.min())
 185
 186 #Postprocessing :
 187 #================
 188 # Save 3D picture
 189 PV_routines.Save_PV_data_to_picture_file("FiniteElementsOnTorusPoisson"+'_0.vtu',"Numerical result field",'NODES',"FiniteElementsOnTorusPoisson")
 190 resolution=100
 191 VTK_routines.Clip_VTK_data_to_VTK("FiniteElementsOnTorusPoisson"+'_0.vtu',"Clip_VTK_data_to_VTK_"+ "FiniteElementsOnTorusPoisson"+'_0.vtu',[0.25,0.25,0.25], [-0.5,-0.5,-0.5],resolution )
 192 PV_routines.Save_PV_data_to_picture_file("Clip_VTK_data_to_VTK_"+"FiniteElementsOnTorusPoisson"+'_0.vtu',"Numerical result field",'NODES',"Clip_VTK_data_to_VTK_"+"FiniteElementsOnTorusPoisson")
 193
 194 # Plot  over slice circle
 195 finiteElementsOnTorus_0vtu = pvs.XMLUnstructuredGridReader(FileName=["FiniteElementsOnTorusPoisson"+'_0.vtu'])
 196 slice1 = pvs.Slice(Input=finiteElementsOnTorus_0vtu)
 197 slice1.SliceType.Normal = [0.5, 0.5, 0.5]
 198 renderView1 = pvs.GetActiveViewOrCreate('RenderView')
 199 finiteElementsOnTorus_0vtuDisplay = pvs.Show(finiteElementsOnTorus_0vtu, renderView1)
 200 pvs.ColorBy(finiteElementsOnTorus_0vtuDisplay, ('POINTS', 'Numerical result field'))
 201 slice1Display = pvs.Show(slice1, renderView1)
 202 pvs.SaveScreenshot("./FiniteElementsOnTorusPoisson"+"_Slice"+'.png', magnification=1, quality=100, view=renderView1)
 203 plotOnSortedLines1 = pvs.PlotOnSortedLines(Input=slice1)
 204 lineChartView2 = pvs.CreateView('XYChartView')
 205 plotOnSortedLines1Display = pvs.Show(plotOnSortedLines1, lineChartView2)
 206 plotOnSortedLines1Display.UseIndexForXAxis = 0
 207 plotOnSortedLines1Display.XArrayName = 'arc_length'
 208 plotOnSortedLines1Display.SeriesVisibility = ['Numerical result field (1)']
 209 pvs.SaveScreenshot("./FiniteElementsOnTorusPoisson"+"_PlotOnSortedLine_"+'.png', magnification=1, quality=100, view=lineChartView2)
 210 pvs.Delete(lineChartView2)
 211
 212 assert erreur_max/max_sol_exacte <1.