CDMATH/tests/examples/Poisson2DEF_DISK/FiniteElements2DPoisson_DISK.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Poisson 2D -\triangle u = f sur le disque unité avec conditions aux limites de Dirichlet u=0
   4 # Author      : Michaël 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 #               Comparaison de la solution numérique avec la solution exacte u=-(r-1)*r**2*sin(theta)
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import sin, sqrt, atan2
  13 import numpy as np
  14 import matplotlib.pyplot as plt
  15 import PV_routines
  16 import VTK_routines
  17
  18 #Chargement du maillage triangulaire du disque unité
  19 #=======================================================================================
  20 my_mesh = cdmath.Mesh("diskWithTriangles.med")
  21 if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  22     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  23 if(not my_mesh.isTriangular()) :
  24         raise ValueError("Wrong cell types : mesh is not made of triangles")
  25
  26 nbNodes = my_mesh.getNumberOfNodes()
  27 nbCells = my_mesh.getNumberOfCells()
  28
  29 print("Mesh loading done")
  30 print("Number of nodes=", nbNodes)
  31 print("Number of cells=", nbCells)
  32
  33 #Discrétisation du second membre et détermination des noeuds intérieurs
  34 #======================================================================
  35 my_RHSfield = cdmath.Field("RHS_field", cdmath.NODES, my_mesh, 1)
  36 my_ExactSol = cdmath.Field("Exact_field", cdmath.NODES, my_mesh, 1)
  37 nbInteriorNodes = 0
  38 nbBoundaryNodes = 0
  39 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  40 interiorNodes=[]
  41 boundaryNodes=[]
  42
  43 #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
  44 for i in range(nbNodes):
  45         Ni=my_mesh.getNode(i)
  46         x = Ni.x()
  47         y = Ni.y()
  48
  49         r=sqrt(x*x+y*y)
  50         theta=atan2(y,x)
  51
  52         my_RHSfield[i]=(8*r-3)*sin(theta)#mettre la fonction definie au second membre de l'edp
  53         my_ExactSol[i]=-(r-1)*r*r*sin(theta)#mettre la solution exacte de l'edp
  54
  55         if my_mesh.isBorderNode(i): # Détection des noeuds frontière
  56                 boundaryNodes.append(i)
  57                 nbBoundaryNodes=nbBoundaryNodes+1
  58         else: # Détection des noeuds intérieurs
  59                 interiorNodes.append(i)
  60                 nbInteriorNodes=nbInteriorNodes+1
  61                 maxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only in 2D, otherwise use function Ni.getNumberOfEdges()
  62
  63 print("Right hand side discretisation done")
  64 print("Number of interior nodes=", nbInteriorNodes)
  65 print("Number of boundary nodes=", nbBoundaryNodes)
  66 print("Maximum number of neighbours=", maxNbNeighbours)
  67
  68 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
  69 #=======================================================================================
  70 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
  71 RHS=cdmath.Vector(nbInteriorNodes)
  72
  73 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle (hypothèse 2D)
  74 GradShapeFunc0=cdmath.Vector(2)
  75 GradShapeFunc1=cdmath.Vector(2)
  76 GradShapeFunc2=cdmath.Vector(2)
  77
  78 #On parcourt les triangles du domaine
  79 for i in range(nbCells):
  80
  81         Ci=my_mesh.getCell(i)
  82
  83         #Contribution à la matrice de rigidité
  84         nodeId0=Ci.getNodeId(0)
  85         nodeId1=Ci.getNodeId(1)
  86         nodeId2=Ci.getNodeId(2)
  87
  88         N0=my_mesh.getNode(nodeId0)
  89         N1=my_mesh.getNode(nodeId1)
  90         N2=my_mesh.getNode(nodeId2)
  91
  92         #Formule des gradients voir EF P1 -> calcul déterminants
  93         GradShapeFunc0[0]= (N1.y()-N2.y())/2
  94         GradShapeFunc0[1]=-(N1.x()-N2.x())/2
  95         GradShapeFunc1[0]=-(N0.y()-N2.y())/2
  96         GradShapeFunc1[1]= (N0.x()-N2.x())/2
  97         GradShapeFunc2[0]= (N0.y()-N1.y())/2
  98         GradShapeFunc2[1]=-(N0.x()-N1.x())/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
 106         # Remplissage de  la matrice de rigidité et du second membre
 107         for j in [nodeId0,nodeId1,nodeId2] :
 108                 if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
 109                         j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
 110                         #Ajout de la contribution de la cellule triangulaire i au second membre du noeud j
 111                         RHS[j_int]=Ci.getMeasure()/3*my_RHSfield[j]+RHS[j_int] # intégrale dans le triangle du produit f x fonction de base
 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("Linear system matrix building done")
 120
 121 # Résolution du système linéaire
 122 #=================================
 123 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","ILU")#Remplacer CG par CHOLESKY pour solveur direct
 124 SolSyst=LS.solve()
 125
 126 print( "Preconditioner used : ", LS.getNameOfPc() )
 127 print( "Number of iterations used : ", LS.getNumberOfIter() )
 128 print( "Final residual : ", LS.getResidu() )
 129 print("Linear system solved")
 130
 131 # Création du champ résultat
 132 #===========================
 133 my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
 134 for j in range(nbInteriorNodes):
 135     my_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
 136 for j in range(nbBoundaryNodes):
 137     my_ResultField[boundaryNodes[j]]=0;#remplissage des valeurs pour les noeuds frontière (condition limite)
 138 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 139 my_ResultField.writeVTK("FiniteElements2DPoisson_DISK_ResultField")
 140
 141 # Postprocessing :
 142 #=================
 143 # save 2D picture
 144 PV_routines.Save_PV_data_to_picture_file("FiniteElements2DPoisson_DISK_ResultField"+'_0.vtu',"ResultField",'NODES',"FiniteElements2DPoisson_DISK_ResultField")
 145
 146 # extract and plot diagonal values
 147 resolution=100
 148 curv_abs=np.linspace(-1,1,resolution+1)
 149 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[-0.5,-0.5,0],[0.5,0.5,0], resolution)
 150 plt.plot(curv_abs, diag_data, label= '2D disk mesh with '+str(nbNodes) + ' nodes')
 151 plt.legend()
 152 plt.xlabel('Position on diagonal line')
 153 plt.ylabel('Value on diagonal line')
 154 plt.title('Plot over diagonal line for finite elements \n for Laplace operator on a 2D disk triangular mesh')
 155 plt.savefig("FiniteElements2DPoisson_DISK_ResultField_"+str(nbNodes) + '_nodes'+"_PlotOverDiagonalLine.png")
 156
 157 print("Numerical solution of 2D Poisson equation on a disk using finite elements done")
 158
 159 #Calcul de l'erreur commise par rapport à la solution exacte
 160 #===========================================================
 161 max_abs_sol_exacte=max(my_ExactSol.max(),-my_ExactSol.min())
 162 max_sol_num=my_ResultField.max()
 163 min_sol_num=my_ResultField.min()
 164 erreur_abs=0
 165 for i in range(nbNodes) :
 166     if  erreur_abs < abs(my_ExactSol[i] - my_ResultField[i]) :
 167         erreur_abs = abs(my_ExactSol[i] - my_ResultField[i])
 168
 169 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
 170 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 171 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 172 print("Maximum exact solution = ", my_ExactSol.max(), " Minimum exact solution = ", my_ExactSol.min() )
 173
 174 assert erreur_abs/max_abs_sol_exacte <1.