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 #================================================================================================================================
12 from math import sin, sqrt, atan2
14 import matplotlib.pyplot as plt
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")
26 nbNodes = my_mesh.getNumberOfNodes()
27 nbCells = my_mesh.getNumberOfCells()
29 print("Mesh loading done")
30 print("Number of nodes=", nbNodes)
31 print("Number of cells=", nbCells)
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)
39 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
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):
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
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()
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)
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)
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)
78 #On parcourt les triangles du domaine
79 for i in range(nbCells):
83 #Contribution à la matrice de rigidité
84 nodeId0=Ci.getNodeId(0)
85 nodeId1=Ci.getNodeId(1)
86 nodeId2=Ci.getNodeId(2)
88 N0=my_mesh.getNode(nodeId0)
89 N1=my_mesh.getNode(nodeId1)
90 N2=my_mesh.getNode(nodeId2)
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
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
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
119 print("Linear system matrix building done")
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
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")
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")
144 PV_routines.Save_PV_data_to_picture_file("FiniteElements2DPoisson_DISK_ResultField"+'_0.vtu',"ResultField",'NODES',"FiniteElements2DPoisson_DISK_ResultField")
146 # extract and plot diagonal values
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')
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")
157 print("Numerical solution of 2D Poisson equation on a disk using finite elements done")
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()
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])
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() )
174 assert erreur_abs/max_abs_sol_exacte <1.