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

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