CDMATH/tests/examples/Poisson3DCubeSkinEF/FiniteElements3DPoissonCubeSkin.py

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