+++ /dev/null
-# -*-coding:utf-8 -*
-#===============================================================================================================================
-# Name : Résolution EF de l'équation de Laplace-Beltrami -\triangle u = f sur la frontière d'un cube
-# Author : Michael Ndjinga
-# Copyright : CEA Saclay 2021
-# 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
-# Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
-# Résolution d'un système linéaire à matrice singulière : les vecteurs constants sont dans le noyau
-# Comparaison de la solution numérique avec la solution exacte définie face par face : u(x,y,z)= cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
-#================================================================================================================================
-
-import cdmath
-from math import cos, pi
-import numpy as np
-import PV_routines
-import VTK_routines
-import paraview.simple as pvs
-
-#Chargement du maillage triangulaire de la frontière du cube unité [0,1]x[0,1]x[0,1]
-#=======================================================================================
-my_mesh = cdmath.Mesh("meshCubeSkin.med")
-if(not my_mesh.isTriangular()) :
- raise ValueError("Wrong cell types : mesh is not made of triangles")
-if(my_mesh.getMeshDimension()!=2) :
- raise ValueError("Wrong mesh dimension : expected a surface of dimension 2")
-if(my_mesh.getSpaceDimension()!=3) :
- raise ValueError("Wrong space dimension : expected a space of dimension 3")
-
-nbNodes = my_mesh.getNumberOfNodes()
-nbCells = my_mesh.getNumberOfCells()
-
-print("Mesh building/loading done")
-print("nb of nodes=", nbNodes)
-print("nb of cells=", nbCells)
-
-#Discrétisation du second membre et détermination des noeuds intérieurs
-#======================================================================
-my_RHSfield = cdmath.Field("RHS field", cdmath.NODES, my_mesh, 1)
-maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
-
-eps=1e-6
-#parcours des noeuds pour discrétisation du second membre et extraction du nb max voisins d'un noeud
-for i in range(nbNodes):
- Ni=my_mesh.getNode(i)
- x = Ni.x()
- y = Ni.y()
- z = Ni.z()
-
- my_RHSfield[i]= 8*pi*pi*cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
-
- if my_mesh.isBorderNode(i): # Détection des noeuds frontière
- raise ValueError("Mesh should not contain borders")
- else:
- maxNbNeighbours = max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only for planar cells, otherwise use function Ni.getNumberOfEdges()
-
-print("Right hand side discretisation done")
-print("Max nb of neighbours=", maxNbNeighbours)
-print("Integral of the RHS", my_RHSfield.integral(0))
-
-# Construction de la matrice de rigidité et du vecteur second membre du système linéaire
-#=======================================================================================
-Rigidite=cdmath.SparseMatrixPetsc(nbNodes,nbNodes,maxNbNeighbours)# warning : third argument is number of non zero coefficients per line
-RHS=cdmath.Vector(nbNodes)
-
-# Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle
-GradShapeFunc0=cdmath.Vector(3)
-GradShapeFunc1=cdmath.Vector(3)
-GradShapeFunc2=cdmath.Vector(3)
-
-normalFace0=cdmath.Vector(3)
-normalFace1=cdmath.Vector(3)
-
-#On parcourt les triangles du domaine
-for i in range(nbCells):
-
- Ci=my_mesh.getCell(i)
-
- #Contribution à la matrice de rigidité
- nodeId0=Ci.getNodeId(0)
- nodeId1=Ci.getNodeId(1)
- nodeId2=Ci.getNodeId(2)
- N0=my_mesh.getNode(nodeId0)
- N1=my_mesh.getNode(nodeId1)
- N2=my_mesh.getNode(nodeId2)
-
- #Build normal to cell Ci
- normalFace0[0]=Ci.getNormalVector(0,0)
- normalFace0[1]=Ci.getNormalVector(0,1)
- normalFace0[2]=Ci.getNormalVector(0,2)
- normalFace1[0]=Ci.getNormalVector(1,0)
- normalFace1[1]=Ci.getNormalVector(1,1)
- normalFace1[2]=Ci.getNormalVector(1,2)
-
- normalCell = normalFace0.crossProduct(normalFace1)
- normalCell = normalCell*(1/normalCell.norm())
-
- cellMat=cdmath.Matrix(4)
- cellMat[0,0]=N0.x()
- cellMat[0,1]=N0.y()
- cellMat[0,2]=N0.z()
- cellMat[1,0]=N1.x()
- cellMat[1,1]=N1.y()
- cellMat[1,2]=N1.z()
- cellMat[2,0]=N2.x()
- cellMat[2,1]=N2.y()
- cellMat[2,2]=N2.z()
- cellMat[3,0]=normalCell[0]
- cellMat[3,1]=normalCell[1]
- cellMat[3,2]=normalCell[2]
- cellMat[0,3]=1
- cellMat[1,3]=1
- cellMat[2,3]=1
- cellMat[3,3]=0
-
- #Formule des gradients voir EF P1 -> calcul déterminants
- GradShapeFunc0[0]= cellMat.partMatrix(0,0).determinant()*0.5
- GradShapeFunc0[1]=-cellMat.partMatrix(0,1).determinant()*0.5
- GradShapeFunc0[2]= cellMat.partMatrix(0,2).determinant()*0.5
- GradShapeFunc1[0]=-cellMat.partMatrix(1,0).determinant()*0.5
- GradShapeFunc1[1]= cellMat.partMatrix(1,1).determinant()*0.5
- GradShapeFunc1[2]=-cellMat.partMatrix(1,2).determinant()*0.5
- GradShapeFunc2[0]= cellMat.partMatrix(2,0).determinant()*0.5
- GradShapeFunc2[1]=-cellMat.partMatrix(2,1).determinant()*0.5
- GradShapeFunc2[2]= cellMat.partMatrix(2,2).determinant()*0.5
-
- #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
- GradShapeFuncs={nodeId0 : GradShapeFunc0}
- GradShapeFuncs[nodeId1]=GradShapeFunc1
- GradShapeFuncs[nodeId2]=GradShapeFunc2
-
- # Remplissage de la matrice de rigidité et du second membre
- for j in [nodeId0,nodeId1,nodeId2] :
- #Ajout de la contribution de la cellule triangulaire i au second membre du noeud j
- RHS[j]=Ci.getMeasure()/3*my_RHSfield[j]+RHS[j] # intégrale dans le triangle du produit f x fonction de base
- #Contribution de la cellule triangulaire i à la ligne j du système linéaire
- for k in [nodeId0,nodeId1,nodeId2] :
- Rigidite.addValue(j,k,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
-
-print("Linear system matrix building done")
-
-# Conditionnement de la matrice de rigidité
-#=================================
-cond = Rigidite.getConditionNumber(True)
-print("Condition number is ",cond)
-
-# Résolution du système linéaire
-#=================================
-LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"GMRES","ILU")
-LS.setMatrixIsSingular()#En raison de l'absence de bord
-SolSyst=LS.solve()
-print("Preconditioner used : ", LS.getNameOfPc() )
-print("Number of iterations used : ", LS.getNumberOfIter() )
-print("Final residual : ", LS.getResidu() )
-print("Linear system solved")
-
-# Création du champ résultat
-#===========================
-my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
-for j in range(nbNodes):
- my_ResultField[j]=SolSyst[j];#remplissage des valeurs issues du système linéaire dans le champs résultat
-#sauvegarde sur le disque dur du résultat dans un fichier paraview
-my_ResultField.writeVTK("FiniteElementsOnCubeSkinPoisson")
-my_RHSfield.writeVTK("RHS_CubeSkinPoisson")
-
-print("Integral of the numerical solution", my_ResultField.integral(0))
-print("Numerical solution of Poisson equation on a cube skin using finite elements done")
-
-#Calcul de l'erreur commise par rapport à la solution exacte
-#===========================================================
-#The following formulas use the fact that the exact solution is equal the right hand side divided by 8*pi*pi
-max_abs_sol_exacte=0
-erreur_abs=0
-max_sol_num=0
-min_sol_num=0
-for i in range(nbNodes) :
- if max_abs_sol_exacte < abs(my_RHSfield[i]) :
- max_abs_sol_exacte = abs(my_RHSfield[i])
- if erreur_abs < abs(my_RHSfield[i]/(8*pi*pi) - my_ResultField[i]) :
- erreur_abs = abs(my_RHSfield[i]/(8*pi*pi) - my_ResultField[i])
- if max_sol_num < my_ResultField[i] :
- max_sol_num = my_ResultField[i]
- if min_sol_num > my_ResultField[i] :
- min_sol_num = my_ResultField[i]
-max_abs_sol_exacte = max_abs_sol_exacte/(8*pi*pi)
-
-print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
-print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
-print("Maximum exact solution = ", my_RHSfield.max()/(8*pi*pi), " Minimum exact solution = ", my_RHSfield.min()/(8*pi*pi) )
-
-#Postprocessing :
-#================
-# save 3D picture
-PV_routines.Save_PV_data_to_picture_file("FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"ResultField",'NODES',"FiniteElementsOnCubeSkinPoisson")
-resolution=100
-VTK_routines.Clip_VTK_data_to_VTK("FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"Clip_VTK_data_to_VTK_"+ "FiniteElementsOnCubeSkinPoisson"+'_0.vtu',[0.75,0.75,0.75], [0.,0.5,-0.5],resolution )
-PV_routines.Save_PV_data_to_picture_file("Clip_VTK_data_to_VTK_"+"FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"ResultField",'NODES',"Clip_VTK_data_to_VTK_"+"FiniteElementsOnCubeSkinPoisson")
-
-# Plot over slice circle
-finiteElementsOnCubeSkin_0vtu = pvs.XMLUnstructuredGridReader(FileName=["FiniteElementsOnCubeSkinPoisson"+'_0.vtu'])
-slice1 = pvs.Slice(Input=finiteElementsOnCubeSkin_0vtu)
-slice1.SliceType.Normal = [0, 1, 0]
-renderView1 = pvs.GetActiveViewOrCreate('RenderView')
-finiteElementsOnCubeSkin_0vtuDisplay = pvs.Show(finiteElementsOnCubeSkin_0vtu, renderView1)
-pvs.ColorBy(finiteElementsOnCubeSkin_0vtuDisplay, ('POINTS', 'ResultField'))
-slice1Display = pvs.Show(slice1, renderView1)
-pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_Slice"+'.png', magnification=1, quality=100, view=renderView1)
-plotOnSortedLines1 = pvs.PlotOnSortedLines(Input=slice1)
-lineChartView2 = pvs.CreateView('XYChartView')
-plotOnSortedLines1Display = pvs.Show(plotOnSortedLines1, lineChartView2)
-plotOnSortedLines1Display.UseIndexForXAxis = 0
-plotOnSortedLines1Display.XArrayName = 'arc_length'
-plotOnSortedLines1Display.SeriesVisibility = ['ResultField (1)']
-pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_PlotOnSortedLine_"+'.png', magnification=1, quality=100, view=lineChartView2)
-pvs.Delete(lineChartView2)
-
-assert erreur_abs/max_abs_sol_exacte <1.
+++ /dev/null
-import FiniteElements3DPoissonCubeSkin
-import matplotlib
-matplotlib.use("Agg")
-import matplotlib.pyplot as plt
-import numpy as np
-from math import log10, sqrt
-import time, json
-
-convergence_synthesis=dict(FiniteElements3DPoissonCubeSkin.test_desc)
-def test_validation3DCubeSkinEF():
- start = time.time()
- #### 3D cube skin FE triangle mesh
- meshList=['CubeSkin_1','CubeSkin_2','CubeSkin_3','CubeSkin_4','CubeSkin_5','CubeSkin_6','CubeSkin_7','CubeSkin_8']
- meshType="Unstructured_3D_triangles"
- testColor="Green"
- nbMeshes=len(meshList)
- error_tab=[0]*nbMeshes
- mesh_size_tab=[0]*nbMeshes
- time_tab=[0]*nbMeshes
- mesh_path='../../../ressources/3DCubeSkin/'
- mesh_name='CubeSkinWithTriangles'
- diag_data=[0]*nbMeshes
- resolution=100
- plt.close('all')
- i=0
- # Storing of numerical errors and mesh sizes
- for filename in meshList:
- error_tab[i], mesh_size_tab[i], min_sol_num, max_sol_num, time_tab[i] =FiniteElements3DPoissonCubeSkin.solve(mesh_path+filename, resolution,meshType,testColor)
- assert min_sol_num>-5.
- assert max_sol_num<6.
- error_tab[i]=log10(error_tab[i])
- time_tab[i]=log10(time_tab[i])
- # with open('./FiniteElementsOnCubeSkinPoisson_PlotOnSortedLines'+meshType+str(mesh_size_tab[i])+'.csv') as f:
- # lines = f.readlines()
- # y = [float(line.split(",")[0]) for line in lines[1:]]
- # x = [float(line.split(",")[1]) for line in lines[1:]]
-
- # plt.plot(x, y, label= str(mesh_size_tab[i]) + ' nodes')
- mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
- i=i+1
-
- end = time.time()
-
- # Plot over diagonal line
- plt.legend()
- plt.xlabel('Position on slice circle')
- plt.ylabel('Value on slice circle')
- plt.title('Plot over slice circle for finite elements \n for Laplace operator on 3D cube skin meshes')
- plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_Slice.png")
-
- # Least square linear regression
- # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
- # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
- a1=np.dot(mesh_size_tab,mesh_size_tab)
- a2=np.sum(mesh_size_tab)
- a3=nbMeshes
- b1=np.dot(error_tab,mesh_size_tab)
- b2=np.sum(error_tab)
-
- det=a1*a3-a2*a2
- assert det!=0, 'test_validation3DCubeSkinEF() : Make sure you use distinct meshes and at least two meshes'
- a=( a3*b1-a2*b2)/det
- b=(-a2*b1+a1*b2)/det
-
- print( "FE on 3D cube skin triangle mesh : scheme order is ", -a)
- assert abs(a+1.915)<0.001
-
- # Plot of convergence curves
- plt.close()
- plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
- plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
- plt.legend()
- plt.xlabel('log(sqrt(number of nodes))')
- plt.ylabel('log(error)')
- plt.title('Convergence of finite elements for \n Laplace operator on 3D cube skin triangular meshes')
- plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_ConvergenceCurve.png")
-
- # Plot of computational time
- plt.close()
- plt.plot(mesh_size_tab, time_tab, label='log(cpu time))')
- plt.legend()
- plt.xlabel('log(sqrt(number of nodes)')
- plt.ylabel('log(cpu time)')
- plt.title('Computational time of finite elements \n for Laplace operator on 3D cube skin triangular meshes')
- plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_ComputationalTime.png")
-
- plt.close('all')
-
- convergence_synthesis["Mesh_names"]=meshList
- convergence_synthesis["Mesh_type"]=meshType
- convergence_synthesis["Mesh_path"]=mesh_path
- convergence_synthesis["Mesh_description"]=mesh_name
- convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
- convergence_synthesis["Space_dimension"]=3
- convergence_synthesis["Mesh_dimension"]=2
- convergence_synthesis["Mesh_cell_type"]="3DTriangles"
- convergence_synthesis["Errors"]=[10**x for x in error_tab]
- convergence_synthesis["Scheme_order"]=-a
- convergence_synthesis["Test_color"]=testColor
- convergence_synthesis["Computational_time"]=end-start
-
- with open('Convergence_Poisson_32DFV_'+mesh_name+'.json', 'w') as outfile:
- json.dump(convergence_synthesis, outfile)
-
-if __name__ == """__main__""":
- test_validation3DCubeSkinEF()
--- /dev/null
+# -*-coding:utf-8 -*
+#===============================================================================================================================
+# Name : Résolution EF de l'équation de Laplace-Beltrami -\triangle u = f sur la frontière d'un cube
+# Author : Michael Ndjinga
+# Copyright : CEA Saclay 2021
+# 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
+# Création et sauvegarde du champ résultant ainsi que du champ second membre en utilisant la librairie CDMATH
+# Résolution d'un système linéaire à matrice singulière : les vecteurs constants sont dans le noyau
+# Comparaison de la solution numérique avec la solution exacte définie face par face : u(x,y,z)= cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
+#================================================================================================================================
+
+import cdmath
+from math import cos, pi
+import numpy as np
+import PV_routines
+import VTK_routines
+import paraview.simple as pvs
+
+#Chargement du maillage triangulaire de la frontière du cube unité [0,1]x[0,1]x[0,1]
+#=======================================================================================
+my_mesh = cdmath.Mesh("meshCubeSkin.med")
+if(not my_mesh.isTriangular()) :
+ raise ValueError("Wrong cell types : mesh is not made of triangles")
+if(my_mesh.getMeshDimension()!=2) :
+ raise ValueError("Wrong mesh dimension : expected a surface of dimension 2")
+if(my_mesh.getSpaceDimension()!=3) :
+ raise ValueError("Wrong space dimension : expected a space of dimension 3")
+
+nbNodes = my_mesh.getNumberOfNodes()
+nbCells = my_mesh.getNumberOfCells()
+
+print("Mesh building/loading done")
+print("nb of nodes=", nbNodes)
+print("nb of cells=", nbCells)
+
+#Discrétisation du second membre et détermination des noeuds intérieurs
+#======================================================================
+my_RHSfield = cdmath.Field("RHS field", cdmath.NODES, my_mesh, 1)
+maxNbNeighbours = 0#This is to determine the number of non zero coefficients in the sparse finite element rigidity matrix
+
+eps=1e-6
+#parcours des noeuds pour discrétisation du second membre et extraction du nb max voisins d'un noeud
+for i in range(nbNodes):
+ Ni=my_mesh.getNode(i)
+ x = Ni.x()
+ y = Ni.y()
+ z = Ni.z()
+
+ my_RHSfield[i]= 8*pi*pi*cos(2*pi*x)*cos(2*pi*y)*cos(2*pi*z)
+
+ if my_mesh.isBorderNode(i): # Détection des noeuds frontière
+ raise ValueError("Mesh should not contain borders")
+ else:
+ maxNbNeighbours = max(1+Ni.getNumberOfCells(),maxNbNeighbours) #true only for planar cells, otherwise use function Ni.getNumberOfEdges()
+
+print("Right hand side discretisation done")
+print("Max nb of neighbours=", maxNbNeighbours)
+print("Integral of the RHS", my_RHSfield.integral(0))
+
+# Construction de la matrice de rigidité et du vecteur second membre du système linéaire
+#=======================================================================================
+Rigidite=cdmath.SparseMatrixPetsc(nbNodes,nbNodes,maxNbNeighbours)# warning : third argument is number of non zero coefficients per line
+RHS=cdmath.Vector(nbNodes)
+
+# Vecteurs gradient de la fonction de forme associée à chaque noeud d'un triangle
+GradShapeFunc0=cdmath.Vector(3)
+GradShapeFunc1=cdmath.Vector(3)
+GradShapeFunc2=cdmath.Vector(3)
+
+normalFace0=cdmath.Vector(3)
+normalFace1=cdmath.Vector(3)
+
+#On parcourt les triangles du domaine
+for i in range(nbCells):
+
+ Ci=my_mesh.getCell(i)
+
+ #Contribution à la matrice de rigidité
+ nodeId0=Ci.getNodeId(0)
+ nodeId1=Ci.getNodeId(1)
+ nodeId2=Ci.getNodeId(2)
+ N0=my_mesh.getNode(nodeId0)
+ N1=my_mesh.getNode(nodeId1)
+ N2=my_mesh.getNode(nodeId2)
+
+ #Build normal to cell Ci
+ normalFace0[0]=Ci.getNormalVector(0,0)
+ normalFace0[1]=Ci.getNormalVector(0,1)
+ normalFace0[2]=Ci.getNormalVector(0,2)
+ normalFace1[0]=Ci.getNormalVector(1,0)
+ normalFace1[1]=Ci.getNormalVector(1,1)
+ normalFace1[2]=Ci.getNormalVector(1,2)
+
+ normalCell = normalFace0.crossProduct(normalFace1)
+ normalCell = normalCell*(1/normalCell.norm())
+
+ cellMat=cdmath.Matrix(4)
+ cellMat[0,0]=N0.x()
+ cellMat[0,1]=N0.y()
+ cellMat[0,2]=N0.z()
+ cellMat[1,0]=N1.x()
+ cellMat[1,1]=N1.y()
+ cellMat[1,2]=N1.z()
+ cellMat[2,0]=N2.x()
+ cellMat[2,1]=N2.y()
+ cellMat[2,2]=N2.z()
+ cellMat[3,0]=normalCell[0]
+ cellMat[3,1]=normalCell[1]
+ cellMat[3,2]=normalCell[2]
+ cellMat[0,3]=1
+ cellMat[1,3]=1
+ cellMat[2,3]=1
+ cellMat[3,3]=0
+
+ #Formule des gradients voir EF P1 -> calcul déterminants
+ GradShapeFunc0[0]= cellMat.partMatrix(0,0).determinant()*0.5
+ GradShapeFunc0[1]=-cellMat.partMatrix(0,1).determinant()*0.5
+ GradShapeFunc0[2]= cellMat.partMatrix(0,2).determinant()*0.5
+ GradShapeFunc1[0]=-cellMat.partMatrix(1,0).determinant()*0.5
+ GradShapeFunc1[1]= cellMat.partMatrix(1,1).determinant()*0.5
+ GradShapeFunc1[2]=-cellMat.partMatrix(1,2).determinant()*0.5
+ GradShapeFunc2[0]= cellMat.partMatrix(2,0).determinant()*0.5
+ GradShapeFunc2[1]=-cellMat.partMatrix(2,1).determinant()*0.5
+ GradShapeFunc2[2]= cellMat.partMatrix(2,2).determinant()*0.5
+
+ #Création d'un tableau (numéro du noeud, gradient de la fonction de forme
+ GradShapeFuncs={nodeId0 : GradShapeFunc0}
+ GradShapeFuncs[nodeId1]=GradShapeFunc1
+ GradShapeFuncs[nodeId2]=GradShapeFunc2
+
+ # Remplissage de la matrice de rigidité et du second membre
+ for j in [nodeId0,nodeId1,nodeId2] :
+ #Ajout de la contribution de la cellule triangulaire i au second membre du noeud j
+ RHS[j]=Ci.getMeasure()/3*my_RHSfield[j]+RHS[j] # intégrale dans le triangle du produit f x fonction de base
+ #Contribution de la cellule triangulaire i à la ligne j du système linéaire
+ for k in [nodeId0,nodeId1,nodeId2] :
+ Rigidite.addValue(j,k,GradShapeFuncs[j]*GradShapeFuncs[k]/Ci.getMeasure())
+
+print("Linear system matrix building done")
+
+# Conditionnement de la matrice de rigidité
+#=================================
+cond = Rigidite.getConditionNumber(True)
+print("Condition number is ",cond)
+
+# Résolution du système linéaire
+#=================================
+LS=cdmath.LinearSolver(Rigidite,RHS,100,1.E-6,"GMRES","ILU")
+LS.setMatrixIsSingular()#En raison de l'absence de bord
+SolSyst=LS.solve()
+print("Preconditioner used : ", LS.getNameOfPc() )
+print("Number of iterations used : ", LS.getNumberOfIter() )
+print("Final residual : ", LS.getResidu() )
+print("Linear system solved")
+
+# Création du champ résultat
+#===========================
+my_ResultField = cdmath.Field("ResultField", cdmath.NODES, my_mesh, 1)
+for j in range(nbNodes):
+ my_ResultField[j]=SolSyst[j];#remplissage des valeurs issues du système linéaire dans le champs résultat
+#sauvegarde sur le disque dur du résultat dans un fichier paraview
+my_ResultField.writeVTK("FiniteElementsOnCubeSkinPoisson")
+my_RHSfield.writeVTK("RHS_CubeSkinPoisson")
+
+print("Integral of the numerical solution", my_ResultField.integral(0))
+print("Numerical solution of Poisson equation on a cube skin using finite elements done")
+
+#Calcul de l'erreur commise par rapport à la solution exacte
+#===========================================================
+#The following formulas use the fact that the exact solution is equal the right hand side divided by 8*pi*pi
+max_abs_sol_exacte=0
+erreur_abs=0
+max_sol_num=0
+min_sol_num=0
+for i in range(nbNodes) :
+ if max_abs_sol_exacte < abs(my_RHSfield[i]) :
+ max_abs_sol_exacte = abs(my_RHSfield[i])
+ if erreur_abs < abs(my_RHSfield[i]/(8*pi*pi) - my_ResultField[i]) :
+ erreur_abs = abs(my_RHSfield[i]/(8*pi*pi) - my_ResultField[i])
+ if max_sol_num < my_ResultField[i] :
+ max_sol_num = my_ResultField[i]
+ if min_sol_num > my_ResultField[i] :
+ min_sol_num = my_ResultField[i]
+max_abs_sol_exacte = max_abs_sol_exacte/(8*pi*pi)
+
+print("Relative error = max(| exact solution - numerical solution |)/max(| exact solution |) = ",erreur_abs/max_abs_sol_exacte)
+print("Maximum numerical solution = ", max_sol_num, " Minimum numerical solution = ", min_sol_num)
+print("Maximum exact solution = ", my_RHSfield.max()/(8*pi*pi), " Minimum exact solution = ", my_RHSfield.min()/(8*pi*pi) )
+
+#Postprocessing :
+#================
+# save 3D picture
+PV_routines.Save_PV_data_to_picture_file("FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"ResultField",'NODES',"FiniteElementsOnCubeSkinPoisson")
+resolution=100
+VTK_routines.Clip_VTK_data_to_VTK("FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"Clip_VTK_data_to_VTK_"+ "FiniteElementsOnCubeSkinPoisson"+'_0.vtu',[0.75,0.75,0.75], [0.,0.5,-0.5],resolution )
+PV_routines.Save_PV_data_to_picture_file("Clip_VTK_data_to_VTK_"+"FiniteElementsOnCubeSkinPoisson"+'_0.vtu',"ResultField",'NODES',"Clip_VTK_data_to_VTK_"+"FiniteElementsOnCubeSkinPoisson")
+
+# Plot over slice circle
+finiteElementsOnCubeSkin_0vtu = pvs.XMLUnstructuredGridReader(FileName=["FiniteElementsOnCubeSkinPoisson"+'_0.vtu'])
+slice1 = pvs.Slice(Input=finiteElementsOnCubeSkin_0vtu)
+slice1.SliceType.Normal = [0, 1, 0]
+renderView1 = pvs.GetActiveViewOrCreate('RenderView')
+finiteElementsOnCubeSkin_0vtuDisplay = pvs.Show(finiteElementsOnCubeSkin_0vtu, renderView1)
+pvs.ColorBy(finiteElementsOnCubeSkin_0vtuDisplay, ('POINTS', 'ResultField'))
+slice1Display = pvs.Show(slice1, renderView1)
+pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_Slice"+'.png', magnification=1, quality=100, view=renderView1)
+plotOnSortedLines1 = pvs.PlotOnSortedLines(Input=slice1)
+lineChartView2 = pvs.CreateView('XYChartView')
+plotOnSortedLines1Display = pvs.Show(plotOnSortedLines1, lineChartView2)
+plotOnSortedLines1Display.UseIndexForXAxis = 0
+plotOnSortedLines1Display.XArrayName = 'arc_length'
+plotOnSortedLines1Display.SeriesVisibility = ['ResultField (1)']
+pvs.SaveScreenshot("./FiniteElementsOnCubeSkinPoisson"+"_PlotOnSortedLine_"+'.png', magnification=1, quality=100, view=lineChartView2)
+pvs.Delete(lineChartView2)
+
+assert erreur_abs/max_abs_sol_exacte <1.
--- /dev/null
+import FiniteElements3DPoissonCubeSkin
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from math import log10, sqrt
+import time, json
+
+convergence_synthesis=dict(FiniteElements3DPoissonCubeSkin.test_desc)
+def test_validation3DCubeSkinEF():
+ start = time.time()
+ #### 3D cube skin FE triangle mesh
+ meshList=['CubeSkin_1','CubeSkin_2','CubeSkin_3','CubeSkin_4','CubeSkin_5','CubeSkin_6','CubeSkin_7','CubeSkin_8']
+ meshType="Unstructured_3D_triangles"
+ testColor="Green"
+ nbMeshes=len(meshList)
+ error_tab=[0]*nbMeshes
+ mesh_size_tab=[0]*nbMeshes
+ time_tab=[0]*nbMeshes
+ mesh_path='../../../ressources/3DCubeSkin/'
+ mesh_name='CubeSkinWithTriangles'
+ diag_data=[0]*nbMeshes
+ resolution=100
+ plt.close('all')
+ i=0
+ # Storing of numerical errors and mesh sizes
+ for filename in meshList:
+ error_tab[i], mesh_size_tab[i], min_sol_num, max_sol_num, time_tab[i] =FiniteElements3DPoissonCubeSkin.solve(mesh_path+filename, resolution,meshType,testColor)
+ assert min_sol_num>-5.
+ assert max_sol_num<6.
+ error_tab[i]=log10(error_tab[i])
+ time_tab[i]=log10(time_tab[i])
+ # with open('./FiniteElementsOnCubeSkinPoisson_PlotOnSortedLines'+meshType+str(mesh_size_tab[i])+'.csv') as f:
+ # lines = f.readlines()
+ # y = [float(line.split(",")[0]) for line in lines[1:]]
+ # x = [float(line.split(",")[1]) for line in lines[1:]]
+
+ # plt.plot(x, y, label= str(mesh_size_tab[i]) + ' nodes')
+ mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
+ i=i+1
+
+ end = time.time()
+
+ # Plot over diagonal line
+ plt.legend()
+ plt.xlabel('Position on slice circle')
+ plt.ylabel('Value on slice circle')
+ plt.title('Plot over slice circle for finite elements \n for Laplace operator on 3D cube skin meshes')
+ plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_Slice.png")
+
+ # Least square linear regression
+ # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
+ # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
+ a1=np.dot(mesh_size_tab,mesh_size_tab)
+ a2=np.sum(mesh_size_tab)
+ a3=nbMeshes
+ b1=np.dot(error_tab,mesh_size_tab)
+ b2=np.sum(error_tab)
+
+ det=a1*a3-a2*a2
+ assert det!=0, 'test_validation3DCubeSkinEF() : Make sure you use distinct meshes and at least two meshes'
+ a=( a3*b1-a2*b2)/det
+ b=(-a2*b1+a1*b2)/det
+
+ print( "FE on 3D cube skin triangle mesh : scheme order is ", -a)
+ assert abs(a+1.915)<0.001
+
+ # Plot of convergence curves
+ plt.close()
+ plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
+ plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
+ plt.legend()
+ plt.xlabel('log(sqrt(number of nodes))')
+ plt.ylabel('log(error)')
+ plt.title('Convergence of finite elements for \n Laplace operator on 3D cube skin triangular meshes')
+ plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_ConvergenceCurve.png")
+
+ # Plot of computational time
+ plt.close()
+ plt.plot(mesh_size_tab, time_tab, label='log(cpu time))')
+ plt.legend()
+ plt.xlabel('log(sqrt(number of nodes)')
+ plt.ylabel('log(cpu time)')
+ plt.title('Computational time of finite elements \n for Laplace operator on 3D cube skin triangular meshes')
+ plt.savefig(mesh_name+"_3DCubeSkinPoissonFE_ComputationalTime.png")
+
+ plt.close('all')
+
+ convergence_synthesis["Mesh_names"]=meshList
+ convergence_synthesis["Mesh_type"]=meshType
+ convergence_synthesis["Mesh_path"]=mesh_path
+ convergence_synthesis["Mesh_description"]=mesh_name
+ convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
+ convergence_synthesis["Space_dimension"]=3
+ convergence_synthesis["Mesh_dimension"]=2
+ convergence_synthesis["Mesh_cell_type"]="3DTriangles"
+ convergence_synthesis["Errors"]=[10**x for x in error_tab]
+ convergence_synthesis["Scheme_order"]=-a
+ convergence_synthesis["Test_color"]=testColor
+ convergence_synthesis["Computational_time"]=end-start
+
+ with open('Convergence_Poisson_32DFV_'+mesh_name+'.json', 'w') as outfile:
+ json.dump(convergence_synthesis, outfile)
+
+if __name__ == """__main__""":
+ test_validation3DCubeSkinEF()
