CDMATH/tests/examples/Poisson3DSphereEF/FiniteElements3DPoissonSphere.py

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