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

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