CDMATH/tests/examples/Poisson3DEF_BALL/FiniteElements3DPoisson_BALL.py

   1 # -*-coding:utf-8 -*
   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 #================================================================================================================================
  10
  11 import cdmath
  12 from math import sin, cos, atan2, sqrt
  13 import numpy as np
  14 import matplotlib.pyplot as plt
  15 import PV_routines
  16 import VTK_routines
  17
  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")
  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_SOL", 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         z = Ni.z()
  49
  50         r=sqrt(x*x+y*y+z*z)
  51         phi=atan2(y,x)
  52         theta=atan2(y,z*sin(phi))
  53
  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)
  63
  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)
  68
  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)
  73
  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)
  79
  80 #On parcourt les triangles du domaine
  81 for i in range(nbCells):
  82
  83         Ci=my_mesh.getCell(i)
  84
  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)
  94
  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
 108
 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
 114
 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
 127
 128 print("Linear system matrix building done")
 129
 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
 133 SolSyst=LS.solve()
 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")
 138
 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
 144
 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")
 149
 150 #Postprocessing :
 151 #================
 152 # save 3D picture
 153 resolution=100
 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")
 156
 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')
 161 plt.legend()
 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")
 166
 167 print("Numerical solution of 3D Poisson equation on a 3D ball using finite elements done")
 168
 169
 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()
 175 erreur_abs=0
 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])
 179
 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)
 183
 184 assert erreur_abs/max_abs_sol_exacte <1.