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

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