From: michael Date: Mon, 7 Jun 2021 18:14:08 +0000 (+0200) Subject: Added mini report on Poisson problem on cube skin X-Git-Tag: V9_8_0~98 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=f01b94a584db5ec34484812d71f63be4cea10098;p=tools%2Fsolverlab.git Added mini report on Poisson problem on cube skin --- diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ComputationalTime.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ComputationalTime.png new file mode 100644 index 0000000..7c837de Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ComputationalTime.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ConvergenceCurve.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ConvergenceCurve.png new file mode 100644 index 0000000..1e551e0 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkinWithTriangles_3DCubeSkinPoissonFE_ConvergenceCurve.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_3.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_3.png new file mode 100755 index 0000000..91a6609 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_3.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_5.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_5.png new file mode 100755 index 0000000..f7d1217 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_5.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_7.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_7.png new file mode 100755 index 0000000..e8d5240 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_7.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_8.png b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_8.png new file mode 100755 index 0000000..3a66fb7 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/CubeSkin_8.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElements3DPoissonCubeSkin.py b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElements3DPoissonCubeSkin.py new file mode 100644 index 0000000..e167456 --- /dev/null +++ b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElements3DPoissonCubeSkin.py @@ -0,0 +1,216 @@ +# -*-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. diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_13225.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_13225.png new file mode 100644 index 0000000..4db4f9a Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_13225.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_1923.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_1923.png new file mode 100644 index 0000000..8e28cfe Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_1923.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_412.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_412.png new file mode 100644 index 0000000..aa662ce Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_412.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_7042.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_7042.png new file mode 100644 index 0000000..02ea4e9 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_7042.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_13225.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_13225.png new file mode 100644 index 0000000..4a32d91 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_13225.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_1923.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_1923.png new file mode 100644 index 0000000..59865fe Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_1923.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_412.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_412.png new file mode 100644 index 0000000..d18b5d2 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_412.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_7042.png b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_7042.png new file mode 100644 index 0000000..ae28204 Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/FiniteElementsOnCubeSkinPoisson_Clip_7042.png differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/PoissonCubeSkin.pdf b/CDMATH/tests/doc/3DCubeSkinEF/PoissonCubeSkin.pdf new file mode 100644 index 0000000..748c6ca Binary files /dev/null and b/CDMATH/tests/doc/3DCubeSkinEF/PoissonCubeSkin.pdf differ diff --git a/CDMATH/tests/doc/3DCubeSkinEF/test_validation3DCubeSkinPoissonEF.py b/CDMATH/tests/doc/3DCubeSkinEF/test_validation3DCubeSkinPoissonEF.py new file mode 100755 index 0000000..7d2f71b --- /dev/null +++ b/CDMATH/tests/doc/3DCubeSkinEF/test_validation3DCubeSkinPoissonEF.py @@ -0,0 +1,106 @@ +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()