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

   1 import FiniteElements2DPoisson_SQUARE
   2 import matplotlib
   3 matplotlib.use("Agg")
   4 import matplotlib.pyplot as plt
   5 import numpy as np
   6 from math import log10, sqrt
   7 import time, json
   8
   9 convergence_synthesis=dict(FiniteElements2DPoisson_SQUARE.test_desc)
  10
  11 def test_validation2DEF_skinny_triangles():
  12     start = time.time()
  13     #### 2D FE skinny triangle mesh
  14     #meshList=[5,9,15,21,31]
  15     meshList=['squareWithSkinnyTriangles_0','squareWithSkinnyTriangles_1','squareWithSkinnyTriangles_2','squareWithSkinnyTriangles_3']#,'squareWithSkinnyTriangles_4'
  16     mesh_path='../../../ressources/2DSkinnyTriangles/'
  17     meshType="Regular_skinny_triangles"
  18     mesh_name='squareWithSkinnyTriangles'
  19     testColor="Green"
  20     nbMeshes=len(meshList)
  21     error_tab=[0]*nbMeshes
  22     mesh_size_tab=[0]*nbMeshes
  23     diag_data=[0]*nbMeshes
  24     time_tab=[0]*nbMeshes
  25     max_tab=[0]*nbMeshes
  26     min_tab=[0]*nbMeshes
  27     resolution=100
  28     curv_abs=np.linspace(0,sqrt(2),resolution+1)
  29     plt.close('all')
  30     i=0
  31     # Storing of numerical errors, mesh sizes and diagonal values
  32     for filename in meshList:
  33     #for nx in meshList:
  34         #my_mesh=cdmath.Mesh(0,1,nx,0,1,nx*nx)#Use script provided by Adrien to split each quadrangle in 2 triangles
  35         error_tab[i], mesh_size_tab[i], diag_data[i], min_tab[i], max_tab[i], time_tab[i] =FiniteElements2DPoisson_SQUARE.solve(mesh_path+filename, resolution, meshType, testColor)
  36         assert min_tab[i]>-0.01
  37         assert max_tab[i]<1.2
  38         plt.plot(curv_abs, diag_data[i], label= str(mesh_size_tab[i]) + ' nodes')
  39         error_tab[i]=log10(error_tab[i])
  40         time_tab[i]=log10(time_tab[i])
  41         i=i+1
  42
  43     end = time.time()
  44
  45     # Plot over diagonal line
  46     plt.legend()
  47     plt.xlabel('Position on diagonal line')
  48     plt.ylabel('Value on diagonal line')
  49     plt.title('Plot over diagonal line for finite elements \n for Laplace operator on 2D skinny triangle meshes')
  50     plt.savefig(mesh_name+"_2DPoissonEF_PlotOverDiagonalLine.png")
  51
  52
  53     # Plot of min and max curves
  54     plt.close()
  55     plt.plot(mesh_size_tab, min_tab, label='Minimum value')
  56     plt.plot(mesh_size_tab, max_tab, label='Maximum value')
  57     plt.legend()
  58     plt.xlabel('Number of nodes')
  59     plt.ylabel('Value')
  60     plt.title('Min/Max of Finite elements \n for Laplace operator on 2D skinny triangle meshes')
  61     plt.savefig(mesh_name+"_2DPoissonEF_MinMax.png")
  62
  63     for i in range(nbMeshes):
  64         mesh_size_tab[i] = 0.5*log10(mesh_size_tab[i])
  65
  66     # Least square linear regression
  67     # Find the best a,b such that f(x)=ax+b best approximates the convergence curve
  68     # The vector X=(a,b) solves a symmetric linear system AX=B with A=(a1,a2\\a2,a3), B=(b1,b2)
  69     a1=np.dot(mesh_size_tab,mesh_size_tab)
  70     a2=np.sum(mesh_size_tab)
  71     a3=nbMeshes
  72     b1=np.dot(error_tab,mesh_size_tab)
  73     b2=np.sum(error_tab)
  74
  75     det=a1*a3-a2*a2
  76     assert det!=0, 'test_validation2DEF_skinny_triangles() : Make sure you use distinct meshes and at least two meshes'
  77     a=( a3*b1-a2*b2)/det
  78     b=(-a2*b1+a1*b2)/det
  79
  80     print( "FE on 2D skinny triangle mesh : scheme order is ", -a )
  81     assert abs(a+1.397)<0.1
  82
  83     # Plot of convergence curves
  84     plt.close()
  85     plt.plot(mesh_size_tab, error_tab, label='log(|numerical-exact|)')
  86     plt.plot(mesh_size_tab, a*np.array(mesh_size_tab)+b,label='least square slope : '+'%.3f' % a)
  87     plt.legend()
  88     plt.xlabel('log(sqrt(number of nodes))')
  89     plt.ylabel('log(error)')
  90     plt.title('Convergence of finite elements \n for Laplace operator on 2D skinny triangle meshes')
  91     plt.savefig(mesh_name+"_2DPoissonEF_ConvergenceCurve.png")
  92
  93     # Plot of computational time
  94     plt.close()
  95     plt.plot(mesh_size_tab, time_tab, label='log(cpu time)')
  96     plt.legend()
  97     plt.xlabel('log(sqrt(number of nodes))')
  98     plt.ylabel('log(cpu time)')
  99     plt.title('Computational time of finite elements \n for Laplace operator on 2D skinny triangle meshes')
 100     plt.savefig(mesh_name+"_2DPoissonEF_ComputationalTime.png")
 101
 102     plt.close('all')
 103
 104     convergence_synthesis["Mesh_names"]=meshList
 105     convergence_synthesis["Mesh_type"]=meshType
 106     convergence_synthesis["Mesh_path"]=mesh_path
 107     convergence_synthesis["Mesh_description"]=mesh_name
 108     convergence_synthesis["Mesh_sizes"]=[10**x for x in mesh_size_tab]
 109     convergence_synthesis["Space_dimension"]=2
 110     convergence_synthesis["Mesh_dimension"]=2
 111     convergence_synthesis["Mesh_cell_type"]="Triangles"
 112     convergence_synthesis["Errors"]=[10**x for x in error_tab]
 113     convergence_synthesis["Scheme_order"]=-a
 114     convergence_synthesis["Test_color"]=testColor
 115     convergence_synthesis["Computational_time"]=end-start
 116
 117     with open('Convergence_Poisson_2DFE_'+mesh_name+'.json', 'w') as outfile:
 118         json.dump(convergence_synthesis, outfile)
 119
 120 if __name__ == """__main__""":
 121     test_validation2DEF_skinny_triangles()