CDMATH/tests/examples/Poisson1DEF/FiniteElements1DPoisson.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Poisson 1D -\triangle u = f avec conditions aux limites de Dirichlet u=0
   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 1D quelconque
   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=-sin(pi*x)
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import sin, pi, sqrt
  13 from numpy import linspace
  14 import matplotlib.pyplot as plt
  15 import PV_routines
  16 import VTK_routines
  17
  18 #Création d'un maillage uniforme du segment [0,1], définition des bords
  19 #======================================================================
  20 nx=100
  21 my_mesh = cdmath.Mesh(0,1,nx)
  22 if( my_mesh.getSpaceDimension()!=1 or my_mesh.getMeshDimension()!=1) :
  23     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 1")
  24
  25 eps=1e-6
  26 my_mesh.setGroupAtPlan(0.,0,eps,"DirichletBorder")#Bord GAUCHE
  27 my_mesh.setGroupAtPlan(1.,0,eps,"DirichletBorder")#Bord DROIT
  28
  29 nbNodes = my_mesh.getNumberOfNodes()
  30 nbCells = my_mesh.getNumberOfCells()
  31
  32 print("Mesh loading/building done")
  33 print("Number of nodes=", nbNodes)
  34 print("Number of cells=", nbCells)
  35
  36 #Discrétisation du second membre et détermination des noeuds intérieurs
  37 #======================================================================
  38 my_RHSfield = cdmath.Field("RHS_field", cdmath.NODES, my_mesh, 1)
  39 nbInteriorNodes = 0
  40 nbBoundaryNodes = 0
  41 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  42 interiorNodes=[]
  43 boundaryNodes=[]
  44
  45 #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
  46 for i in range(nbNodes):
  47         Ni=my_mesh.getNode(i)
  48         x = Ni.x()
  49
  50         my_RHSfield[i]=pi*pi*sin(pi*x)#mettre la fonction definie au second membre de l'edp
  51         if my_mesh.isBorderNode(i): # Détection des noeuds frontière
  52                 boundaryNodes.append(i)
  53                 nbBoundaryNodes=nbBoundaryNodes+1
  54         else: # Détection des noeuds intérieurs
  55                 interiorNodes.append(i)
  56                 nbInteriorNodes=nbInteriorNodes+1
  57                 maxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours)
  58
  59 print("Right hand side discretisation done")
  60 print("nb of interior nodes=", nbInteriorNodes)
  61 print("nb of boundary nodes=", nbBoundaryNodes)
  62 print("Max nb of neighbours=", maxNbNeighbours)
  63
  64 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
  65 #=======================================================================================
  66 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours)# warning : third argument is max number of non zero coefficients per line of the matrix
  67 RHS=cdmath.Vector(nbInteriorNodes)
  68
  69 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un segment (hypothèse 1D)
  70 GradShapeFunc0=cdmath.Vector(1)
  71 GradShapeFunc1=cdmath.Vector(1)
  72
  73 #On parcourt les segments du domaine
  74 for i in range(nbCells):
  75
  76         Ci=my_mesh.getCell(i)
  77
  78         #Contribution à la matrice de rigidité
  79         nodeId0=Ci.getNodeId(0)
  80         nodeId1=Ci.getNodeId(1)
  81
  82         N0=my_mesh.getNode(nodeId0)
  83         N1=my_mesh.getNode(nodeId1)
  84
  85         #Formule des gradients voir EF P1 -> calcul déterminants
  86         GradShapeFunc0[0]= 1
  87         GradShapeFunc1[0]=-1
  88
  89         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
  90         GradShapeFuncs={nodeId0 : GradShapeFunc0}
  91         GradShapeFuncs[nodeId1]=GradShapeFunc1
  92
  93         # Remplissage de  la matrice de rigidité et du second membre
  94         for j in [nodeId0,nodeId1] :
  95                 if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
  96                         j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
  97                         #Ajout de la contribution de la cellule i au second membre du noeud j
  98                         RHS[j_int]=Ci.getMeasure()/2*my_RHSfield[j]+RHS[j_int] # intégrale sur le segment du produit f x fonction de base
  99                         #Contribution de la cellule i à la ligne j_int du système linéaire
 100                         for k in [nodeId0,nodeId1] :
 101                                 if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
 102                                         k_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
 103                                         Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 104                                 #else: si condition limite non nulle au bord, ajouter la contribution du bord au second membre de la cellule j
 105
 106 print("Linear system matrix building done")
 107
 108 # Résolution du système linéaire
 109 #=================================
 110 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","ILU")#Remplacer CG par CHOLESKY pour solveur direct
 111 SolSyst=LS.solve()
 112
 113 print( "Preconditioner used : ", LS.getNameOfPc() )
 114 print( "Number of iterations used : ", LS.getNumberOfIter() )
 115 print( "Final residual : ", LS.getResidu() )
 116 print("Linear system solved")
 117
 118 # Création du champ résultat
 119 #===========================
 120 my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
 121 for j in range(nbInteriorNodes):
 122     my_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
 123 for j in range(nbBoundaryNodes):
 124     my_ResultField[boundaryNodes[j]]=0;#remplissage des valeurs pour les noeuds frontière (condition limite)
 125 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 126 my_ResultField.writeVTK("FiniteElements1DPoisson_ResultField")
 127
 128 # Postprocessing :
 129 #=================
 130 # save 1D picture
 131 PV_routines.Save_PV_data_to_picture_file("FiniteElements1DPoisson_ResultField"+'_0.vtu',"ResultField",'NODES',"FiniteElements1DPoisson_ResultField")
 132
 133 # extract and plot diagonal values
 134 resolution=100
 135 curv_abs=linspace(0,1,resolution+1)
 136 diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,0,0],[1,0,0], resolution)
 137 plt.plot(curv_abs, diag_data, label= '1D mesh with '+str(nbNodes) + ' nodes')
 138 plt.legend()
 139 plt.xlabel('Position')
 140 plt.ylabel('Value')
 141 plt.title('1D finite elements \n for Laplace operator')
 142 plt.savefig('FiniteElements1DPoisson_ResultField_'+str(nbNodes) + '_nodes'+'.png')
 143
 144 print("Numerical solution of the 1D Poisson equation using finite elements done")
 145
 146 #Calcul de l'erreur commise par rapport à la solution exacte
 147 #===========================================================
 148 #The following formulas use the fact that the exact solution is equal the right hand side divided by pi*pi
 149 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(pi*pi)
 150 max_sol_num=my_ResultField.max()
 151 min_sol_num=my_ResultField.min()
 152 erreur_abs=0
 153 for i in range(nbNodes) :
 154     if erreur_abs < abs(my_RHSfield[i]/(pi*pi) - my_ResultField[i]) :
 155         erreur_abs = abs(my_RHSfield[i]/(pi*pi) - my_ResultField[i])
 156
 157 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
 158 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
 159 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
 160 print("Maximum exact solution = ", my_RHSfield.max()/(pi*pi), " Minimum exact solution = ", my_RHSfield.min()/(pi*pi))
 161
 162 assert erreur_abs/max_abs_sol_exacte <1.