3 #===============================================================================================================================
4 # Name : Résolution EF de l'équation de Poisson 3D -\triangle T = f avec conditions aux limites de Dirichlet T=0
5 # Author : Michaël Ndjinga
6 # Copyright : CEA Saclay 2019
7 # Description : Utilisation de la méthode des éléménts finis P1 avec champs T et f discrétisés aux noeuds d'un maillage tétraédrique
8 # Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
9 # Comparaison de la solution numérique avec la solution exacte T=-sin(pi*x)*sin(pi*y)*sin(pi*z)
10 #================================================================================================================================
12 import CoreFlows as cf
14 from math import sin, pi
16 def StationaryDiffusionEquation_3DEF_StructuredTriangles():
18 # Prepare for the mesh
19 print("Building mesh " );
29 M=cm.Mesh(xinf,xsup,nx,yinf,ysup,ny,zinf,zsup,nz,1)#Regular tetrahedral mesh
30 # set the limit field for each boundary
32 M.setGroupAtPlan(xsup,0,eps,"Bord1")
33 M.setGroupAtPlan(xinf,0,eps,"Bord2")
34 M.setGroupAtPlan(ysup,1,eps,"Bord3")
35 M.setGroupAtPlan(yinf,1,eps,"Bord4")
36 M.setGroupAtPlan(zsup,2,eps,"Bord5")
37 M.setGroupAtPlan(zinf,2,eps,"Bord6")
39 print( "Built a regular tetrahedra 3D mesh from a cube mesh with ", nx,"x" ,ny,"x" ,nz, " cells")
42 myProblem = cf.StationaryDiffusionEquation(spaceDim,FEComputation);
45 # set the limit value for each boundary
53 myProblem.setDirichletBoundaryCondition("Bord1",T1)
54 myProblem.setDirichletBoundaryCondition("Bord2",T2)
55 myProblem.setDirichletBoundaryCondition("Bord3",T3)
56 myProblem.setDirichletBoundaryCondition("Bord4",T4)
57 myProblem.setDirichletBoundaryCondition("Bord5",T5)
58 myProblem.setDirichletBoundaryCondition("Bord6",T6)
60 #Set the right hand side function
61 my_RHSfield = cm.Field("RHS_field", cm.NODES, M, 1)
62 for i in range(M.getNumberOfNodes()):
68 my_RHSfield[i]=2*pi*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)#mettre la fonction definie au second membre de l'edp
70 myProblem.setHeatPowerField(my_RHSfield)
71 myProblem.setLinearSolver(cf.GMRES,cf.ILU);
74 fileName = "StationnaryDiffusion_3DEF_StructuredTriangles";
76 # computation parameters
77 myProblem.setFileName(fileName);
80 myProblem.initialize();
81 print("Running python "+ fileName );
83 ok = myProblem.solveStationaryProblem();
85 print( "Python simulation of " + fileName + " failed ! " );
88 ####################### Postprocessing #########################
89 my_ResultField = myProblem.getOutputTemperatureField()
90 #The following formulas use the fact that the exact solution is equal the right hand side divided by 3*pi*pi
91 max_abs_sol_exacte=max(my_RHSfield.max(),-my_RHSfield.min())/(3*pi*pi)
92 max_sol_num=my_ResultField.max()
93 min_sol_num=my_ResultField.min()
95 for i in range(M.getNumberOfNodes()) :
96 if erreur_abs < abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i]) :
97 erreur_abs = abs(my_RHSfield[i]/(3*pi*pi) - my_ResultField[i])
99 print("Absolute error = max(| exact solution - numerical solution |) = ",erreur_abs )
100 print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
101 print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
103 assert erreur_abs/max_abs_sol_exacte <1.
106 print( "------------ !!! End of calculation !!! -----------" );
108 myProblem.terminate();
111 if __name__ == """__main__""":
112 StationaryDiffusionEquation_3DEF_StructuredTriangles()