2 #===============================================================================================================================
3 # Name : Résolution EF de l'équation de Poisson 3D -\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 éléménts finis P1 avec champs u et f discrétisés aux noeuds d'un maillage tétraédrique
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)**2*cos(phi)
9 #================================================================================================================================
12 from math import sin, cos, atan2, sqrt
14 import matplotlib.pyplot as plt
18 #Chargement du maillage tétraédrique du domaine cubique [0,1]x[0,1]x[0,1], définition des bords
19 #==============================================================================================
20 my_mesh = cdmath.Mesh("ballWithTetrahedra.med")
21 if( my_mesh.getSpaceDimension()!=3 or my_mesh.getMeshDimension()!=3) :
22 raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 3")
23 if(not my_mesh.isTetrahedral()) :
24 raise ValueError("Wrong cell types : mesh is not made of tetrahedra")
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_SOL", 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 theta=atan2(y,z*sin(phi))
54 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
55 my_ExactSol[i]=-(r-1)*r**2*sin(theta)**2*cos(phi)
56 if my_mesh.isBorderNode(i): # Détection des noeuds frontière
57 boundaryNodes.append(i)
58 nbBoundaryNodes=nbBoundaryNodes+1
59 else: # Détection des noeuds intérieurs
60 interiorNodes.append(i)
61 nbInteriorNodes=nbInteriorNodes+1
62 maxNbNeighbours= max(1+Ni.getNumberOfEdges(),maxNbNeighbours)
64 print("Right hand side discretisation done")
65 print("Number of interior nodes=", nbInteriorNodes)
66 print("Number of boundary nodes=", nbBoundaryNodes)
67 print("Maximum number of neighbours per node=", maxNbNeighbours)
69 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
70 #=======================================================================================
71 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours) # warning : third argument is max number of non zero coefficients per line of the matrix
72 RHS=cdmath.Vector(nbInteriorNodes)
74 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un hexaèdre
75 GradShapeFunc0=cdmath.Vector(3)
76 GradShapeFunc1=cdmath.Vector(3)
77 GradShapeFunc2=cdmath.Vector(3)
78 GradShapeFunc3=cdmath.Vector(3)
80 #On parcourt les triangles du domaine
81 for i in range(nbCells):
85 #Contribution à la matrice de rigidité
86 nodeId0=Ci.getNodeId(0)
87 nodeId1=Ci.getNodeId(1)
88 nodeId2=Ci.getNodeId(2)
89 nodeId3=Ci.getNodeId(3)
90 N0=my_mesh.getNode(nodeId0)
91 N1=my_mesh.getNode(nodeId1)
92 N2=my_mesh.getNode(nodeId2)
93 N3=my_mesh.getNode(nodeId3)
95 #Formule des gradients voir EF P1 -> calcul déterminants
96 GradShapeFunc0[0]= (N2.y()*N3.z()-N2.z()*N3.y()-N1.y()*N3.z()+N3.y()*N1.z()+N1.y()*N2.z()-N2.y()*N1.z())/6
97 GradShapeFunc0[1]=-(N2.x()*N3.z()-N2.z()*N3.x()-N1.x()*N3.z()+N3.x()*N1.z()+N1.x()*N2.z()-N2.x()*N1.z())/6
98 GradShapeFunc0[2]=(N2.x()*N3.y()-N2.y()*N3.x()-N1.x()*N3.y()+N3.x()*N1.y()+N1.x()*N2.y()-N2.x()*N1.y())/6
99 GradShapeFunc1[0]=- (N2.y()*N3.z()-N2.z()*N3.y()-N0.y()*N3.z()+N3.y()*N0.z()+N0.y()*N2.z()-N2.y()*N0.z())/6
100 GradShapeFunc1[1]=(N2.x()*N3.z()-N2.z()*N3.x()-N0.x()*N3.z()+N3.x()*N0.z()+N0.x()*N2.z()-N2.x()*N0.z())/6
101 GradShapeFunc1[2]=-(N2.x()*N3.y()-N2.y()*N3.x()-N0.x()*N3.y()+N3.x()*N0.y()+N0.x()*N2.y()-N2.x()*N0.y())/6
102 GradShapeFunc2[0]= -(N0.y()*N3.z()-N0.z()*N3.y()-N1.y()*N3.z()+N3.y()*N1.z()+N1.y()*N0.z()-N0.y()*N1.z())/6
103 GradShapeFunc2[1]=(N0.x()*N3.z()-N0.z()*N3.x()-N1.x()*N3.z()+N3.x()*N1.z()+N1.x()*N0.z()-N0.x()*N1.z())/6
104 GradShapeFunc2[2]= -(N0.x()*N3.y()-N0.y()*N3.x()-N1.x()*N3.y()+N3.x()*N1.y()+N1.x()*N0.y()-N0.x()*N1.y())/6
105 GradShapeFunc3[0]=-(N2.y()*N0.z()-N2.z()*N0.y()-N1.y()*N0.z()+N0.y()*N1.z()+N1.y()*N2.z()-N2.y()*N1.z())/6
106 GradShapeFunc3[1]=(N2.x()*N0.z()-N2.z()*N0.x()-N1.x()*N0.z()+N0.x()*N1.z()+N1.x()*N2.z()-N2.x()*N1.z())/6
107 GradShapeFunc3[2]=-(N2.x()*N0.y()-N2.y()*N0.x()-N1.x()*N0.y()+N0.x()*N1.y()+N1.x()*N2.y()-N2.x()*N1.y())/6
109 #Création d'un tableau (numéro du noeud, gradient de la fonction de forme)
110 GradShapeFuncs={nodeId0 : GradShapeFunc0}
111 GradShapeFuncs[nodeId1]=GradShapeFunc1
112 GradShapeFuncs[nodeId2]=GradShapeFunc2
113 GradShapeFuncs[nodeId3]=GradShapeFunc3
115 # Remplissage de la matrice de rigidité et du second membre
116 for j in [nodeId0,nodeId1,nodeId2,nodeId3] :
117 if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
118 j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
119 #Ajout de la contribution de la cellule triangulaire i au second membre du noeud j
120 RHS[j_int]=Ci.getMeasure()/4*my_RHSfield[j]+RHS[j_int] # intégrale dans le triangle du produit f x fonction de base
121 #Contribution de la cellule triangulaire i à la ligne j_int du système linéaire
122 for k in [nodeId0,nodeId1,nodeId2,nodeId3] :
123 if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
124 k_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
125 Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
126 #else: si condition limite non nulle au bord, ajouter la contribution du bord au second membre de la cellule j
128 print("Linear system matrix building done")
130 # Résolution du système linéaire
131 #=================================
132 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","ILU")#,"ILU" Remplacer CG par CHOLESKY pour solveur direct
134 print "Preconditioner used : ", LS.getNameOfPc()
135 print "Number of iterations used : ", LS.getNumberOfIter()
136 print "Final residual : ", LS.getResidu()
137 print("Linear system solved")
139 # Création du champ résultat
140 #===========================
141 my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
142 for j in range(nbInteriorNodes):
143 my_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
145 for j in range(nbBoundaryNodes):
146 my_ResultField[boundaryNodes[j]]=0;#remplissage des valeurs pour les noeuds frontière (condition limite)
147 #sauvegarde sur le disque dur du résultat dans un fichier paraview
148 my_ResultField.writeVTK("FiniteElements3DPoisson_BALL_ResultField")
154 VTK_routines.Clip_VTK_data_to_VTK("FiniteElements3DPoisson_BALL_ResultField"+'_0.vtu',"Clip_VTK_data_to_VTK_"+ "FiniteElements3DPoisson_BALL_ResultField"+'_0.vtu',[0.5,0.5,0.5], [-0.5,-0.5,-0.5],resolution )
155 PV_routines.Save_PV_data_to_picture_file("Clip_VTK_data_to_VTK_"+"FiniteElements3DPoisson_BALL_ResultField"+'_0.vtu',"ResultField",'NODES',"Clip_VTK_data_to_VTK_"+"FiniteElements3DPoisson_BALL_ResultField")
157 # extract and plot diagonal values
158 curv_abs=np.linspace(0,sqrt(3),resolution+1)
159 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[-0.577,-0.577,-0.577],[0.577,0.577,0.577], resolution)
160 plt.plot(curv_abs, diag_data, label= str(nbNodes) + ' nodes 3D mesh')
162 plt.xlabel('Position on diagonal line')
163 plt.ylabel('Value on diagonal line')
164 plt.title('Plot over diagonal line for finite elements \n for Laplace operator on a 3D ball tetrahedral mesh')
165 plt.savefig("FiniteElements3DPoisson_BALL_ResultField_"+str(nbNodes) + '_nodes'+"_PlotOverDiagonalLine.png")
167 print("Numerical solution of 3D Poisson equation on a 3D ball using finite elements done")
170 #Calcul de l'erreur commise par rapport à la solution exacte
171 #===========================================================
172 max_abs_sol_exacte=max(my_ExactSol.max(),-my_ExactSol.min())
173 max_sol_num=my_ResultField.max()
174 min_sol_num=my_ResultField.min()
176 for i in range(nbNodes) :
177 if erreur_abs < abs(my_ExactSol[i] - my_ResultField[i]) :
178 erreur_abs = abs(my_ExactSol[i] - my_ResultField[i])
180 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
181 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
182 print ("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
184 assert erreur_abs/max_abs_sol_exacte <1.