2 #===============================================================================================================================
3 # Name : Résolution VF de l'équation de Poisson -\triangle u = f sur la boule 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 volumes finis avec champs u et f discrétisés aux cellules d'un maillage quelconque
7 # Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant CDMATH
8 # Comparaison de la solution numérique avec la solution exacte u=-(r-1)*r**2*sin(theta)**2*cos(phi)
9 #================================================================================================================================
12 from math import sin, cos, pi, sqrt, atan2
14 import matplotlib.pyplot as plt
20 #Chargement du maillage tétraédrique du domaine cubique [0,1]x[0,1]x[0,1], définition des bords
21 #==============================================================================================
22 my_mesh = cdmath.Mesh("ballWithTetrahedra.med")
23 if( my_mesh.getSpaceDimension()!=3 or my_mesh.getMeshDimension()!=3) :
24 raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 3")
27 my_mesh.setGroupAtPlan(0,0,eps,"DirichletBorder")#Bord GAUCHE
28 my_mesh.setGroupAtPlan(1,0,eps,"DirichletBorder")#Bord DROIT
29 my_mesh.setGroupAtPlan(0,1,eps,"DirichletBorder")#Bord BAS
30 my_mesh.setGroupAtPlan(1,1,eps,"DirichletBorder")#Bord HAUT
31 my_mesh.setGroupAtPlan(0,2,eps,"DirichletBorder")#Bord AVANT
32 my_mesh.setGroupAtPlan(1,2,eps,"DirichletBorder")#Bord ARRIERE
34 nbCells = my_mesh.getNumberOfCells()
36 print("Mesh loading done")
37 print("Number of cells ", nbCells)
39 #Discrétisation du second membre et extraction du nb max de voisins d'une cellule
40 #================================================================================
41 my_RHSfield = cdmath.Field("RHS_field", cdmath.CELLS, my_mesh, 1)
42 my_ExactSol = cdmath.Field("EXACT_SOL", cdmath.CELLS, my_mesh, 1)
44 maxNbNeighbours=0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
45 #parcours des cellules pour discrétisation du second membre et extraction du nb max de voisins d'une cellule
46 for i in range(nbCells):
47 Ci = my_mesh.getCell(i)
54 theta=atan2(y,z*sin(phi))
56 my_RHSfield[i]=6*r*sin(theta)**2*cos(phi)+3*(r-1)*cos(phi)#mettre la fonction definie au second membre de l'edp
57 my_ExactSol[i]=-(r-1)*r**2*sin(theta)**2*cos(phi)
58 # compute maximum number of neighbours
59 maxNbNeighbours= max(1+Ci.getNumberOfFaces(),maxNbNeighbours)
61 print("Right hand side discretisation done")
62 print("Maximum number of neighbours=", maxNbNeighbours)
64 # Construction de la matrice et du vecteur second membre du système linéaire
65 #===========================================================================
66 Rigidite=cdmath.SparseMatrixPetsc(nbCells,nbCells,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
67 RHS=cdmath.Vector(nbCells)
68 #Parcours des cellules du domaine
69 for i in range(nbCells):
70 RHS[i]=my_RHSfield[i] #la valeur moyenne du second membre f dans la cellule i
72 for j in range(Ci.getNumberOfFaces()):# parcours des faces voisinnes
73 Fj=my_mesh.getFace(Ci.getFaceId(j))
79 distance=Ci.getBarryCenter().distance(Ck.getBarryCenter())
80 coeff=Fj.getMeasure()/Ci.getMeasure()/distance
81 Rigidite.setValue(i,k,-coeff) # terme extradiagonal
83 coeff=Fj.getMeasure()/Ci.getMeasure()/Ci.getBarryCenter().distance(Fj.getBarryCenter())
84 Rigidite.addValue(i,i,coeff) # terme diagonal
86 print("Linear system matrix building done")
88 # Résolution du système linéaire
89 #=================================
90 LS=cdmath.LinearSolver(Rigidite,RHS,500,1.E-6,"GMRES","ILU")
93 print( "Preconditioner used : ", LS.getNameOfPc() )
94 print( "Number of iterations used : ", LS.getNumberOfIter() )
95 print( "Final residual : ", LS.getResidu() )
96 print("Linear system solved")
98 # Création du champ résultat
99 #===========================
100 my_ResultField = cdmath.Field("ResultField", cdmath.CELLS, my_mesh, 1)
101 for i in range(nbCells):
102 my_ResultField[i]=SolSyst[i];
103 #sauvegarde sur le disque dur du résultat dans un fichier paraview
104 my_ResultField.writeVTK("FiniteVolumes3DPoisson_BALL_ResultField")
110 VTK_routines.Clip_VTK_data_to_VTK("FiniteVolumes3DPoisson_BALL_ResultField"+'_0.vtu',"Clip_VTK_data_to_VTK_"+ "FiniteVolumes3DPoisson_BALL_ResultField"+'_0.vtu',[0.5,0.5,0.5], [-0.5,-0.5,-0.5],resolution )
111 PV_routines.Save_PV_data_to_picture_file("Clip_VTK_data_to_VTK_"+"FiniteVolumes3DPoisson_BALL_ResultField"+'_0.vtu',"ResultField",'CELLS',"Clip_VTK_data_to_VTK_"+"FiniteVolumes3DPoisson_BALL_ResultField")
113 # extract and plot diagonal values
115 curv_abs=np.linspace(0,sqrt(3),resolution+1)
116 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,1,0],[1,0,0], resolution)
118 plt.xlabel('Position on diagonal line')
119 plt.ylabel('Value on diagonal line')
120 if len(sys.argv) >1 :
121 plt.title('Plot over diagonal line for finite volumes \n for Laplace operator on a 3D BALL with mesh '+my_mesh.getName())
122 plt.plot(curv_abs, diag_data, label= str(nbCells)+ ' cells mesh')
123 plt.savefig("FiniteVolumes3DPoisson_BALL_ResultField_"+str(nbCells)+ '_cells'+"_PlotOverDiagonalLine.png")
125 plt.title('Plot over diagonal line for finite volumes \n for Laplace operator on a 3D BALL with a rectangular grid')
126 plt.plot(curv_abs, diag_data, label= str(nbCells)+ ' cells mesh')
127 plt.savefig("FiniteVolumes3DPoisson_BALL_ResultField_"+str(nbCells)+ '_cells'+"_PlotOverDiagonalLine.png")
129 print("Numerical solution of 3D poisson equation on a 3D ball using finite volumes done")
131 #Calcul de l'erreur commise par rapport à la solution exacte
132 #===========================================================
133 max_abs_sol_exacte=max(my_ExactSol.max(),-my_ExactSol.min())
134 max_sol_num=my_ResultField.max()
135 min_sol_num=my_ResultField.min()
137 for i in range(nbCells) :
138 if erreur_abs < abs(my_ExactSol[i] - my_ResultField[i]) :
139 erreur_abs = abs(my_ExactSol[i] - my_ResultField[i])
141 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
142 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
143 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
145 assert erreur_abs/max_abs_sol_exacte <1.