CDMATH/tests/examples/Poisson3DEF_RadiatorAndWindow/FiniteElements3DPoisson_CUBE_RadiatorAndWindow.py

   1 # -*-coding:utf-8 -*
   2 #===============================================================================================================================
   3 # Name        : Résolution EF de l'équation de Laplace 3D -\Delta T = 0 avec conditions aux limites de Dirichlet u non nulle
   4 # Authors     : Michaël Ndjinga, Sédrick Kameni Ngwamou
   5 # Copyright   : CEA Saclay 2019
   6 # Description : Utilisation de la méthode des éléménts finis P1 avec champs u discrétisés aux noeuds d'un maillage tétraédrique
   7 #               Condition limites correspondant au refroidissement dû à une fenêtre et au chauffage dû à un radiateur
   8 #                               Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
   9 #================================================================================================================================
  10
  11 import cdmath
  12 import numpy as np
  13
  14 # Fonction qui remplace successivement les colonnes d'une matrices par un vecteur donné et retourne la liste des déterminants
  15 def gradientNodal(M, values):
  16         matrices=[0]*(len(values)-1)
  17         for i in range(len(values)-1):
  18                 matrices[i] = M.deepCopy()
  19                 for j in range(len(values)):
  20                         matrices[i][j,i] = values[j]
  21
  22         result = cdmath.Vector(len(values)-1)
  23         for i in range(len(values)-1):
  24                 result[i] = matrices[i].determinant()
  25
  26         return result
  27
  28 # Fonction qui calcule la valeur de la condition limite en un noeud donné
  29 def boundaryValue(nodeId):
  30         Ni=my_mesh.getNode(nodeId)
  31
  32         # 4 groupes sont considérés sur le bord
  33         if boundaryNodes.count(nodeId)==0:
  34                 return 0
  35         elif Ni.getGroupName()=='Fenetre':
  36                 return Tfenetre;
  37         elif Ni.getGroupName()=='Radiateur_sous-fenetre':
  38                 return Tradiateur
  39         #elif Ni.getGroupName()=='Radiateur_Devant':
  40                 #return Tradiateur;
  41         #elif Ni.getGroupName()=='Radiateur_droit':
  42                 #return Tradiateur
  43         else:
  44                 return Tmur;
  45
  46 #Chargement du maillage tétraédrique du domaine
  47 #==============================================
  48 my_mesh = cdmath.Mesh("Mesh_RadiatorAndWindow.med")
  49 if( my_mesh.getSpaceDimension()!=3 or my_mesh.getMeshDimension()!=3) :
  50     raise ValueError("Wrong space or mesh dimension : space and mesh dimensions should be 3")
  51 if(not my_mesh.isTetrahedral()) :
  52         raise ValueError("Wrong cell types : mesh is not made of tetrahedra")
  53
  54 nbNodes = my_mesh.getNumberOfNodes()
  55 nbCells = my_mesh.getNumberOfCells()
  56
  57 print("Mesh loading done")
  58 print("nb of nodes=", nbNodes)
  59 print("nb of cells=", nbCells)
  60
  61 #Conditions limites
  62 Tmur=20
  63 Tfenetre=0
  64 Tradiateur=40
  65
  66 #Détermination des noeuds intérieurs
  67 #======================================================================
  68 nbInteriorNodes = 0
  69 nbBoundaryNodes = 0
  70 maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
  71 interiorNodes = []
  72 boundaryNodes = []
  73
  74 #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
  75 for i in range(nbNodes):
  76         Ni=my_mesh.getNode(i)
  77
  78         if my_mesh.isBorderNode(i): # Détection des noeuds frontière getGroupName my_mesh.getNode(i)
  79                 boundaryNodes.append(i)
  80         else: # Détection des noeuds intérieurs
  81                 interiorNodes.append(i)
  82                 maxNbNeighbours= max(1+Ni.getNumberOfEdges(),maxNbNeighbours)
  83
  84 nbInteriorNodes=len(interiorNodes)
  85 nbBoundaryNodes=len(boundaryNodes)
  86
  87
  88 print("nb of interior nodes=", nbInteriorNodes)
  89 print("nb of Boundary nodes=", nbBoundaryNodes)
  90 print("Max nb of neighbours=", maxNbNeighbours)
  91
  92
  93 # Construction de la matrice de rigidité et du vecteur second membre du système linéaire
  94 #=======================================================================================
  95 Rigidite=cdmath.SparseMatrixPetsc(nbInteriorNodes,nbInteriorNodes,maxNbNeighbours)
  96 RHS=cdmath.Vector(nbInteriorNodes)
  97
  98 # Vecteurs gradient de la fonction de forme associée à chaque noeud d'un tétrèdre (hypothèse 3D)
  99 M=cdmath.Matrix(4,4)
 100 GradShapeFunc0=cdmath.Vector(3)
 101 GradShapeFunc1=cdmath.Vector(3)
 102 GradShapeFunc2=cdmath.Vector(3)
 103 GradShapeFunc3=cdmath.Vector(3)
 104
 105 #On parcourt les tétraèdres du domaine pour remplir la matrice
 106 for i in range(nbCells):
 107
 108         Ci=my_mesh.getCell(i)
 109
 110         #Extraction des noeuds de la cellule
 111         nodeId0=Ci.getNodeId(0)
 112         nodeId1=Ci.getNodeId(1)
 113         nodeId2=Ci.getNodeId(2)
 114         nodeId3=Ci.getNodeId(3)
 115         N0=my_mesh.getNode(nodeId0)
 116         N1=my_mesh.getNode(nodeId1)
 117         N2=my_mesh.getNode(nodeId2)
 118         N3=my_mesh.getNode(nodeId3)
 119
 120         M[0,0]=N0.x()
 121         M[0,1]=N0.y()
 122         M[0,2]=N0.z()
 123         M[0,3]=1
 124         M[1,0]=N1.x()
 125         M[1,1]=N1.y()
 126         M[1,2]=N1.z()
 127         M[1,3]=1
 128         M[2,0]=N2.x()
 129         M[2,1]=N2.y()
 130         M[2,2]=N2.z()
 131         M[2,3]=1
 132         M[3,0]=N3.x()
 133         M[3,1]=N3.y()
 134         M[3,2]=N3.z()
 135         M[3,3]=1
 136
 137         #Values of each shape function at each node
 138         values0=[1,0,0,0]
 139         values1=[0,1,0,0]
 140         values2=[0,0,1,0]
 141         values3=[0,0,0,1]
 142
 143         GradShapeFunc0 = gradientNodal(M,values0)/6
 144         GradShapeFunc1 = gradientNodal(M,values1)/6
 145         GradShapeFunc2 = gradientNodal(M,values2)/6
 146         GradShapeFunc3 = gradientNodal(M,values3)/6
 147
 148         #Création d'un tableau (numéro du noeud, gradient de la fonction de forme)
 149         GradShapeFuncs={nodeId0 : GradShapeFunc0}
 150         GradShapeFuncs[nodeId1]=GradShapeFunc1
 151         GradShapeFuncs[nodeId2]=GradShapeFunc2
 152         GradShapeFuncs[nodeId3]=GradShapeFunc3
 153
 154         # Remplissage de  la matrice de rigidité et du second membre
 155         for j in [nodeId0,nodeId1,nodeId2,nodeId3] :
 156                 if boundaryNodes.count(j)==0: #seuls les noeuds intérieurs contribuent au système linéaire (matrice de rigidité et second membre)
 157                         j_int=interiorNodes.index(j)#indice du noeud j en tant que noeud intérieur
 158                         boundaryContributionAdded=False#Needed in case j is a border cell
 159                         #Ajout de la contribution de la cellule ttétraédrique i au second membre du noeud j
 160                         for k in [nodeId0,nodeId1,nodeId2,nodeId3] :
 161                                 if boundaryNodes.count(k)==0 : #seuls les noeuds intérieurs contribuent à la matrice du système linéaire
 162                                         k_int = interiorNodes.index(k)#indice du noeud k en tant que noeud intérieur
 163                                         Rigidite.addValue(j_int,k_int,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
 164                                 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
 165                                         # Valeurs de g_h aux noeuds du tétraèdre
 166                                         T0 = boundaryValue(nodeId0)
 167                                         T1 = boundaryValue(nodeId1)
 168                                         T2 = boundaryValue(nodeId2)
 169                                         T3 = boundaryValue(nodeId3)
 170                                         boundaryContributionAdded=True#Contribution from the boundary to matrix line j is done in one step
 171                                         GradGh = gradientNodal(M,[T0,T1,T2,T3])/6
 172                                         RHS[j_int] += -(GradGh*GradShapeFuncs[j])/Ci.getMeasure()
 173
 174
 175
 176 print("Linear system matrix building done")
 177
 178 LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"CG","ILU")#,"ILU" Remplacer CG par CHOLESKY pour solveur direct
 179
 180 SolSyst=LS.solve()
 181
 182 # Création du champ résultat
 183 #===========================
 184 my_Temperature = cdmath.Field("Temperature", cdmath.NODES, my_mesh, 1)
 185 for j in range(nbInteriorNodes):
 186         my_Temperature[interiorNodes[j]]=SolSyst[j];#remplissage des valeurs pour les noeuds intérieurs SolSyst[j]
 187 #Remplissage des valeurs pour les noeuds frontière (condition limite)
 188 for j in boundaryNodes:
 189     my_Temperature[j]=boundaryValue(j)
 190
 191
 192 #sauvegarde sur le disque dur du résultat dans un fichier paraview
 193 my_Temperature.writeVTK("FiniteElements3DTemperature")
 194
 195 print( "Minimum temperature= ", my_Temperature.min(), ", maximum temperature= ", my_Temperature.max() )
 196 assert my_Temperature.min()>= min(Tmur,Tfenetre,Tradiateur) and my_Temperature.max()<= max(Tmur,Tfenetre,Tradiateur)
 197
 198 print( "Numerical solution of 3D Laplace equation using finite elements done" )
 199