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 #================================================================================================================================
12 from math import sin, pi, sqrt
13 from numpy import linspace
14 import matplotlib.pyplot as plt
18 #Création d'un maillage uniforme du segment [0,1], définition des bords
19 #======================================================================
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")
26 my_mesh.setGroupAtPlan(0.,0,eps,"DirichletBorder")#Bord GAUCHE
27 my_mesh.setGroupAtPlan(1.,0,eps,"DirichletBorder")#Bord DROIT
29 nbNodes = my_mesh.getNumberOfNodes()
30 nbCells = my_mesh.getNumberOfCells()
32 print("Mesh loading/building done")
33 print("Number of nodes=", nbNodes)
34 print("Number of cells=", nbCells)
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)
41 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
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):
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)
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)
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)
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)
73 #On parcourt les segments du domaine
74 for i in range(nbCells):
78 #Contribution à la matrice de rigidité
79 nodeId0=Ci.getNodeId(0)
80 nodeId1=Ci.getNodeId(1)
82 N0=my_mesh.getNode(nodeId0)
83 N1=my_mesh.getNode(nodeId1)
85 #Formule des gradients voir EF P1 -> calcul déterminants
89 #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
90 GradShapeFuncs={nodeId0 : GradShapeFunc0}
91 GradShapeFuncs[nodeId1]=GradShapeFunc1
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
106 print("Linear system matrix building done")
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
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")
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")
131 PV_routines.Save_PV_data_to_picture_file("FiniteElements1DPoisson_ResultField"+'_0.vtu',"ResultField",'NODES',"FiniteElements1DPoisson_ResultField")
133 # extract and plot diagonal values
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')
139 plt.xlabel('Position')
141 plt.title('1D finite elements \n for Laplace operator')
142 plt.savefig("FiniteElements1DPoisson_ResultField_"+str(nbNodes) + '_nodes'+".png")
144 print("Numerical solution of the 1D Poisson equation using finite elements done")
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()
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])
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))
162 assert erreur_abs/max_abs_sol_exacte <1.