CDMATH/tests/examples/Poisson2DEF_DISK_StiffBC/FiniteElements2DPoisson_DISK_StiffBC.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Poisson 2D -\triangle u = 0 sur le disque unité avec conditions aux limites de Dirichlet discontinues
   4 # Author      : Michaël Ndjinga
   5 # Copyright   : CEA Saclay 2019
   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=atan(2*x/(x**2+y**2-1))=atan(2 r cos(theta)/(r*2-1))
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import atan, pi
  13 from numpy import sign, linspace
  14 import matplotlib.pyplot as plt
  15 import PV_routines
  16 import VTK_routines
  17
  18 # Fonction qui remplace successivement les colonnes d'une matrices par un vecteur donné et retourne la liste des déterminants
  19 def gradientNodal(M, values):
  20     matrices=[0]*(len(values)-1)
  21     for i in range(len(values)-1):
  22         matrices[i] = M.deepCopy()
  23         for j in range(len(values)):
  24             matrices[i][j,i] = values[j]
  25
  26     result = cdmath.Vector(len(values)-1)
  27     for i in range(len(values)-1):
  28         result[i] = matrices[i].determinant()
  29
  30     return result
  31
  32
  33 #Chargement du maillage triangulaire du disque unité
  34 #=======================================================================================
  35 meshName="diskWithTriangles"
  36 my_mesh = cdmath.Mesh(meshName+".med")
  37 if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  38     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  39 if(not my_mesh.isTriangular()) :
  40     raise ValueError("Wrong cell types : mesh is not made of triangles")
  41
  42 nbNodes = my_mesh.getNumberOfNodes()
  43 nbCells = my_mesh.getNumberOfCells()
  44
  45 print("Mesh loading done")
  46 print("Number of nodes=", nbNodes)
  47 print("Number of cells=", nbCells)
  48
  49 #Détermination des noeuds intérieurs
  50 #======================================================================
  51 my_ExactSol = cdmath.Field("Exact_field", cdmath.NODES, my_mesh, 1)
  52 nbInteriorNodes = 0
  53 nbBoundaryNodes = 0
  54 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  55 interiorNodes=[]
  56 boundaryNodes=[]
  57 eps=1e-10
  58
  59 #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
  60 for i in range(nbNodes):
  61     Ni=my_mesh.getNode(i)
  62     x = Ni.x()
  63     y = Ni.y()
  64
  65     #Robust calculation of atan(2x/(x**2+y**2-1)
  66     #my_ExactSol[i]=atan2(2*x*sign(x**2+y**2-1),abs(x**2+y**2-1))#mettre la solution exacte de l'edp
  67     if x**2+y**2-1 > eps :
  68         print("!!! Warning Mesh ", meshName," !!! Node is not in the unit disk.",", eps=",eps, ", x**2+y**2-1=",x**2+y**2 - 1)
  69         #raise ValueError("x**2+y**2 > 1 !!! Domain should be the unit disk.")
  70     if x**2+y**2-1 < -eps :
  71         my_ExactSol[i] = atan(2*x/(x**2+y**2-1))
  72     elif x>0 : #x**2+y**2-1>=0
  73         my_ExactSol[i] = -pi/2
  74     elif x<0 : #x**2+y**2-1>=0
  75         my_ExactSol[i] =  pi/2
  76     else : #x=0
  77         my_ExactSol[i] = 0
  78
  79     if my_mesh.isBorderNode(i): # Détection des noeuds frontière
  80         boundaryNodes.append(i)
  81         nbBoundaryNodes=nbBoundaryNodes+1
  82     else: # Détection des noeuds intérieurs
  83         interiorNodes.append(i)
  84         nbInteriorNodes=nbInteriorNodes+1
  85         maxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only in 2D, otherwise use function Ni.getNumberOfEdges()
  86
  87 print("Right hand side discretisation done")
  88 print("Number of interior nodes=", nbInteriorNodes)
  89 print("Number of boundary nodes=", nbBoundaryNodes)
  90 print("Maximum number of neighbours=", maxNbNeighbours)
  91
  92 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
  93 #=======================================================================================
  94 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
  95 RHS=cdmath.Vector(nbInteriorNodes)
  96
  97 M=cdmath.Matrix(3,3)
  98 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle (hypothèse 2D)
  99 GradShapeFunc0=cdmath.Vector(2)
 100 GradShapeFunc1=cdmath.Vector(2)
 101 GradShapeFunc2=cdmath.Vector(2)
 102
 103 #On parcourt les triangles du domaine
 104 for i in range(nbCells):
 105
 106     Ci=my_mesh.getCell(i)
 107
 108     #Contribution à la matrice de rigidité
 109     nodeId0=Ci.getNodeId(0)
 110     nodeId1=Ci.getNodeId(1)
 111     nodeId2=Ci.getNodeId(2)
 112
 113     N0=my_mesh.getNode(nodeId0)
 114     N1=my_mesh.getNode(nodeId1)
 115     N2=my_mesh.getNode(nodeId2)
 116
 117     M[0,0]=N0.x()
 118     M[0,1]=N0.y()
 119     M[0,2]=1
 120     M[1,0]=N1.x()
 121     M[1,1]=N1.y()
 122     M[1,2]=1
 123     M[2,0]=N2.x()
 124     M[2,1]=N2.y()
 125     M[2,2]=1
 126
 127     #Values of each shape function at each node
 128     values0=[1,0,0]
 129     values1=[0,1,0]
 130     values2=[0,0,1]
 131
 132     GradShapeFunc0 = gradientNodal(M,values0)/2
 133     GradShapeFunc1 = gradientNodal(M,values1)/2
 134     GradShapeFunc2 = gradientNodal(M,values2)/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
 142     # Remplissage de  la matrice de rigidité et du second membre
 143     for j in [nodeId0,nodeId1,nodeId2] :
 144         if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
 145             j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
 146             #Pas de contribution au second membre car pas de terme source
 147             boundaryContributionAdded=False#Needed in case j is a border cell
 148             #Contribution de la cellule triangulaire i à la ligne j_int du système linéaire
 149             for k in [nodeId0,nodeId1,nodeId2] :
 150                 if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
 151                     k_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
 152                     Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 153                 elif boundaryContributionAdded == False: #si condition limite non nulle au bord (ou maillage non recouvrant), ajouter la contribution du bord au second membre de la cellule j
 154                     if boundaryNodes.count(nodeId0)!=0 :
 155                         u0=my_ExactSol[nodeId0]
 156                     else:
 157                         u0=0
 158                     if boundaryNodes.count(nodeId1)!=0 :
 159                         u1=my_ExactSol[nodeId1]
 160                     else:
 161                         u1=0
 162                     if boundaryNodes.count(nodeId2)!=0 :
 163                         u2=my_ExactSol[nodeId2]
 164                     else:
 165                         u2=0
 166                     boundaryContributionAdded=True#Contribution from the boundary to matrix line j is done in one step
 167                     GradGh = gradientNodal(M,[u0,u1,u2])/2
 168                     RHS[j_int] += -(GradGh*GradShapeFuncs[j])/Ci.getMeasure()
 169
 170 print("Linear system matrix building done")
 171
 172 # Résolution du système linéaire
 173 #=================================
 174 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","ILU")#Remplacer CG par CHOLESKY pour solveur direct
 175 SolSyst=LS.solve()
 176
 177 print "Preconditioner used : ", LS.getNameOfPc()
 178 print "Number of iterations used : ", LS.getNumberOfIter()
 179 print "Final residual : ", LS.getResidu()
 180 print("Linear system solved")
 181
 182 # Création du champ résultat
 183 #===========================
 184 my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
 185 for j in range(nbInteriorNodes):
 186     my_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
 187 for j in range(nbBoundaryNodes):
 188     my_ResultField[boundaryNodes[j]]=my_ExactSol[boundaryNodes[j]];#remplissage des valeurs pour les noeuds frontière (condition limite)
 189 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 190 my_ResultField.writeVTK("FiniteElements2DPoissonStiffBC_DISK_ResultField")
 191 my_ExactSol.writeVTK("ExactSol2DPoissonStiffBC_DISK")
 192
 193 # Postprocessing :
 194 #=================
 195 # save 2D picture
 196 PV_routines.Save_PV_data_to_picture_file("FiniteElements2DPoissonStiffBC_DISK_ResultField"+'_0.vtu',"ResultField",'NODES',"FiniteElements2DPoissonStiffBC_DISK_ResultField")
 197 PV_routines.Save_PV_data_to_picture_file("ExactSol2DPoissonStiffBC_DISK"+'_0.vtu',"Exact_field",'NODES',"ExactSol2DPoissonStiffBC_DISK")
 198
 199 # extract and plot diagonal values
 200 resolution=100
 201 curv_abs=linspace(-1,1,resolution+1)
 202 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,-1,0],[0,1,0], resolution)
 203 plt.plot(curv_abs, diag_data, label= '2D disk mesh with '+str(nbNodes) + ' nodes')
 204 plt.legend()
 205 plt.xlabel('Position on diagonal line')
 206 plt.ylabel('Value on diagonal line')
 207 plt.title('Plot over diagonal line for finite elements \n for Laplace operator on a 2D disk triangular mesh')
 208 plt.savefig("FiniteElements2DPoissonStiffBC_DISK_ResultField_"+str(nbNodes) + '_nodes'+"_PlotOverDiagonalLine.png")
 209
 210 print("Numerical solution of 2D Poisson equation on a disk using finite elements done")
 211
 212 #Calcul de l'erreur commise par rapport à la solution exacte
 213 #===========================================================
 214 l2_norm_sol_exacte=my_ExactSol.normL2()[0]
 215 l2_error = (my_ExactSol - my_ResultField).normL2()[0]
 216
 217 print("L2 absolute error = norm( exact solution - numerical solution ) = ",l2_error )
 218 print("L2 relative error = norm( exact solution - numerical solution )/norm( exact solution ) = ",l2_error/l2_norm_sol_exacte)
 219 print ("Maximum numerical solution = ", my_ResultField.max(), " Minimum numerical solution = ", my_ResultField.min())
 220 print ("Maximum exact solution = ", my_ExactSol.max(), " Minimum exact solution = ", my_ExactSol.min())
 221
 222 assert l2_error/l2_norm_sol_exacte <1.