SALOME platform Git repositories - tools/solverlab.git/blob

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Poisson 2D -\triangle u = 0 sur un carré avec conditions aux limites de Dirichlet discontinues
   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 triangulaire
   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=atan((x*x-y*y)/(2*x*y))
   9 #================================================================================================================================
  10
  11 import cdmath
  12 from math import atan, pi
  13 import time, json
  14 import PV_routines
  15 import VTK_routines
  16
  17 test_desc={}
  18 test_desc["Initial_data"]="None"
  19 test_desc["Boundary_conditions"]="Dirichlet"
  20 test_desc["Global_name"]="FE simulation of the 2D Poisson equation"
  21 test_desc["Global_comment"]="Condition limite discontinue, Maillage triangulaire"
  22 test_desc["PDE_model"]="Poisson"
  23 test_desc["PDE_is_stationary"]=True
  24 test_desc["PDE_search_for_stationary_solution"]=False
  25 test_desc["Numerical_method_name"]="P1 FE"
  26 test_desc["Numerical_method_space_discretization"]="Finite elements"
  27 test_desc["Numerical_method_time_discretization"]="None"
  28 test_desc["Mesh_is_unstructured"]=True
  29 test_desc["Geometry"]="Square"
  30 test_desc["Part_of_mesh_convergence_analysis"]=True
  31
  32 # Fonction qui remplace successivement les colonnes d'une matrices par un vecteur donné et retourne la liste des déterminants
  33 def gradientNodal(M, values):
  34     matrices=[0]*(len(values)-1)
  35     for i in range(len(values)-1):
  36         matrices[i] = M.deepCopy()
  37         for j in range(len(values)):
  38             matrices[i][j,i] = values[j]
  39
  40     result = cdmath.Vector(len(values)-1)
  41     for i in range(len(values)-1):
  42         result[i] = matrices[i].determinant()
  43
  44     return result
  45
  46 def solve(filename,resolution,meshName, testColor):
  47     start = time.time()
  48     test_desc["Mesh_name"]=meshName
  49     test_desc["Test_color"]=testColor
  50
  51     #Chargement du maillage triangulaire du carré unité [0,1]x[0,1]
  52     #================================================================
  53     my_mesh = cdmath.Mesh(filename+".med")
  54     if( my_mesh.getSpaceDimension()!=2 or my_mesh.getMeshDimension()!=2) :
  55         raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 2")
  56     if(not my_mesh.isTriangular()) :
  57         raise ValueError("Wrong cell types : mesh is not made of triangles")
  58
  59     nbNodes = my_mesh.getNumberOfNodes()
  60     nbCells = my_mesh.getNumberOfCells()
  61
  62     test_desc["Space_dimension"]=my_mesh.getSpaceDimension()
  63     test_desc["Mesh_dimension"]=my_mesh.getMeshDimension()
  64     test_desc["Mesh_number_of_elements"]=my_mesh.getNumberOfNodes()
  65     test_desc["Mesh_cell_type"]=my_mesh.getElementTypesNames()
  66
  67     print("Mesh loading done")
  68     print("Number of nodes=", nbNodes)
  69     print("Number of cells=", nbCells)
  70
  71     #Détermination des noeuds intérieurs
  72     #===================================
  73     my_ExactSol = cdmath.Field("Exact_field", cdmath.NODES, my_mesh, 1)
  74     nbInteriorNodes = 0
  75     nbBoundaryNodes = 0
  76     maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  77     interiorNodes=[]
  78     boundaryNodes=[]
  79     eps=1e-6
  80
  81     #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
  82     for i in range(nbNodes):
  83         Ni=my_mesh.getNode(i)
  84         x = Ni.x()
  85         y = Ni.y()
  86
  87         #Robust calculation of atan((x*x-y*y)/(2*x*y))
  88         if x<-eps or x>1+eps or y<-eps or y>1+eps :
  89             print("!!! Warning Mesh ", meshName," !!! Node is not in the unit square.",", eps=",eps, ", x= ",x, ", y= ",y)
  90             #raise ValueError("!!! Domain should be the unit square.")
  91         if abs(x*y) > eps :
  92             my_ExactSol[i] = atan((x*x-y*y)/(2*x*y))
  93         elif x**2-y**2>0 :
  94             my_ExactSol[i] = pi/2
  95         elif x**2-y**2<0 :
  96             my_ExactSol[i] = -pi/2
  97         else : #x=0
  98             my_ExactSol[i] = 0
  99
 100         if my_mesh.isBorderNode(i): # Détection des noeuds frontière
 101             boundaryNodes.append(i)
 102             nbBoundaryNodes=nbBoundaryNodes+1
 103         else: # Détection des noeuds intérieurs
 104             interiorNodes.append(i)
 105             nbInteriorNodes=nbInteriorNodes+1
 106             maxNbNeighbours= max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only in 2D, need a function Ni.getNumberOfNeighbourNodes()
 107
 108     test_desc["Mesh_max_number_of_neighbours"]=maxNbNeighbours
 109
 110     print("Right hand side discretisation done")
 111     print("Number of interior nodes=", nbInteriorNodes)
 112     print("Number of boundary nodes=", nbBoundaryNodes)
 113     print("Max number of neighbours=", maxNbNeighbours)
 114
 115     # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
 116     #=======================================================================================
 117     Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours) # warning : third argument is maximum number of non zero coefficients per line of the matrix
 118     RHS=cdmath.Vector(nbInteriorNodes)
 119
 120     M=cdmath.Matrix(3,3)
 121     # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle (hypothèse 2D)
 122     GradShapeFunc0=cdmath.Vector(2)
 123     GradShapeFunc1=cdmath.Vector(2)
 124     GradShapeFunc2=cdmath.Vector(2)
 125
 126     #On parcourt les triangles du domaine
 127     for i in range(nbCells):
 128
 129         Ci=my_mesh.getCell(i)
 130
 131         #Contribution à la matrice de rigidité
 132         nodeId0=Ci.getNodeId(0)
 133         nodeId1=Ci.getNodeId(1)
 134         nodeId2=Ci.getNodeId(2)
 135
 136         N0=my_mesh.getNode(nodeId0)
 137         N1=my_mesh.getNode(nodeId1)
 138         N2=my_mesh.getNode(nodeId2)
 139
 140         M[0,0]=N0.x()
 141         M[0,1]=N0.y()
 142         M[0,2]=1
 143         M[1,0]=N1.x()
 144         M[1,1]=N1.y()
 145         M[1,2]=1
 146         M[2,0]=N2.x()
 147         M[2,1]=N2.y()
 148         M[2,2]=1
 149
 150         #Values of each shape function at each node
 151         values0=[1,0,0]
 152         values1=[0,1,0]
 153         values2=[0,0,1]
 154
 155         GradShapeFunc0 = gradientNodal(M,values0)/2
 156         GradShapeFunc1 = gradientNodal(M,values1)/2
 157         GradShapeFunc2 = gradientNodal(M,values2)/2
 158
 159         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme)
 160         GradShapeFuncs={nodeId0 : GradShapeFunc0}
 161         GradShapeFuncs[nodeId1]=GradShapeFunc1
 162         GradShapeFuncs[nodeId2]=GradShapeFunc2
 163
 164
 165         # Remplissage de  la matrice de rigidité et du second membre
 166         for j in [nodeId0,nodeId1,nodeId2] :
 167             if boundaryNodes.count(j)==0 : #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
 168                 j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
 169                 #Pas de contribution au second membre car pas de terme source
 170                 boundaryContributionAdded=False#Needed in case j is a border cell
 171                 #Contribution de la cellule triangulaire i à la ligne j_int du système linéaire
 172                 for k in [nodeId0,nodeId1,nodeId2] :
 173                     if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
 174                         k_int=interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
 175                         Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 176                     elif boundaryContributionAdded == False: #si condition limite non nulle au bord (ou maillage non recouvrant), ajouter la contribution du bord au second membre de la cellule j
 177                         if boundaryNodes.count(nodeId0)!=0 :
 178                             u0=my_ExactSol[nodeId0]
 179                         else:
 180                             u0=0
 181                         if boundaryNodes.count(nodeId1)!=0 :
 182                             u1=my_ExactSol[nodeId1]
 183                         else:
 184                             u1=0
 185                         if boundaryNodes.count(nodeId2)!=0 :
 186                             u2=my_ExactSol[nodeId2]
 187                         else:
 188                             u2=0
 189                         boundaryContributionAdded=True#Contribution from the boundary to matrix line j is done in one step
 190                         GradGh = gradientNodal(M,[u0,u1,u2])/2
 191                         RHS[j_int] += -(GradGh*GradShapeFuncs[j])/Ci.getMeasure()
 192
 193     print("Linear system matrix building done")
 194
 195     # Résolution du système linéaire
 196     #=================================
 197     LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","LU")#Remplacer CG par CHOLESKY pour solveur direct
 198     LS.setComputeConditionNumber()
 199     SolSyst=LS.solve()
 200
 201     print( "Preconditioner used : ", LS.getNameOfPc() )
 202     print( "Number of iterations used : ", LS.getNumberOfIter() )
 203     print("Linear system solved")
 204
 205     test_desc["Linear_solver_algorithm"]=LS.getNameOfMethod()
 206     test_desc["Linear_solver_preconditioner"]=LS.getNameOfPc()
 207     test_desc["Linear_solver_precision"]=LS.getTolerance()
 208     test_desc["Linear_solver_maximum_iterations"]=LS.getNumberMaxOfIter()
 209     test_desc["Linear_system_max_actual_iterations_number"]=LS.getNumberOfIter()
 210     test_desc["Linear_system_max_actual_error"]=LS.getResidu()
 211     test_desc["Linear_system_max_actual_condition number"]=LS.getConditionNumber()
 212
 213
 214     # Création du champ résultat
 215     #===========================
 216     my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
 217     for j in range(nbInteriorNodes):
 218         my_ResultField[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs
 219     for j in range(nbBoundaryNodes):
 220         my_ResultField[boundaryNodes[j]]=my_ExactSol[boundaryNodes[j]];#remplissage des valeurs pour les noeuds frontière (condition limite)
 221     #sauvegarde sur le disque dur du résultat dans un fichier paraview
 222     my_ResultField.writeVTK("FiniteElements2DPoissonStiffBC_SQUARE_"+meshName+str(nbNodes))
 223
 224     print("Numerical solution of 2D Poisson equation on a square using finite elements done")
 225
 226     end = time.time()
 227
 228     #Calcul de l'erreur commise par rapport à la solution exacte
 229     #===========================================================
 230     l2_norm_sol_exacte=my_ExactSol.normL2()[0]
 231     l2_error = (my_ExactSol - my_ResultField).normL2()[0]
 232
 233     print("L2 relative error = norm( exact solution - numerical solution )/norm( exact solution ) = ", l2_error/l2_norm_sol_exacte)
 234     print("Maximum numerical solution = ", my_ResultField.max(), " Minimum numerical solution = ", my_ResultField.min())
 235     print("Maximum exact solution = ", my_ExactSol.max(), " Minimum exact solution = ", my_ExactSol.min())
 236
 237     #Postprocessing :
 238     #================
 239     # Extraction of the diagonal data
 240     diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,0,0],[1,1,0], resolution)
 241     # save 2D picture
 242     PV_routines.Save_PV_data_to_picture_file("FiniteElements2DPoissonStiffBC_SQUARE_"+meshName+str(nbNodes)+'_0.vtu',"ResultField",'NODES',"FiniteElements2DPoissonStiffBC_SQUARE_"+meshName+str(nbNodes))
 243
 244     test_desc["Computational_time_taken_by_run"]=end-start
 245     test_desc["Absolute_error"]=l2_error
 246     test_desc["Relative_error"]=l2_error/l2_norm_sol_exacte
 247
 248     with open('test_PoissonStiffBC'+str(my_mesh.getMeshDimension())+'D_EF_'+"SQUARE_"+meshName+str(nbCells)+ "Cells.json", 'w') as outfile:
 249         json.dump(test_desc, outfile)
 250
 251     return l2_error/l2_norm_sol_exacte, nbNodes, diag_data, my_ResultField.min(), my_ResultField.max(), end - start
 252
 253 if __name__ == """__main__""":
 254     solve("squareWithTriangles",100,"Unstructured_triangles","Green")