Coverage for pyxfem/postpro/interpolator.py: 17%

1''' define interpolator from mesh field ''' 

2from warnings import warn 

3import numpy as np 

4from matplotlib.tri import Triangulation 

5from pyxfem.core.config import npfloat 

6from pyxfem.core.mapping import MappingT3 

7 

8def get_triangulation(xy, t2v): 

9 ''' return a matplotlib.tri.Triangulation from 

10 - a set of point xy and triangle connectivity t2xy ''' 

11 return Triangulation(xy[:,0], xy[:,1], t2v) 

12 

13def get_trifinder( xy, t2v): 

14 ''' return a matplotlib "trifinder " : a function that return the triangle number in 

15 which each point it is called on is in ''' 

16 tri = get_triangulation(xy, t2v) 

17 trif = tri.get_trifinder() 

18 def finder(xyp): 

19 return trif(xyp[...,0], xyp[...,1]) 

20 return finder 

21 

22class ConstantPerTriInterpolator(): 

23 ''' Constant per Element field interpolator : find in which element is a point and return 

24 the corresponding element consant value''' 

25 def __init__(self, xy, t2v, values, trifinder = None): 

26 ''' values are the per element values, mus also pass either a trifinder objet or a sMesh''' 

27 if trifinder is None: 

28 self.trifinder = get_trifinder(xy, t2v) 

29 else : 

30 self.trifinder = trifinder 

31 self.values = values 

32 def __call__(self, xy): 

33 elems = self.trifinder(xy) 

34 return self.values[elems] 

35 

36class LinearPerTriInterpolator(): 

37 ''' linear per Element field interpolator field is supposed C0 on the input mesh, 

38 and linear per element input : value at node and mesh. 

39 ''' 

40 def __init__(self, xy, t2v, values, trifinder = None): 

41 ''' values are the per node values, must also pass either a trifinder objet or a sMesh''' 

42 if trifinder is None: 

43 self.trifinder = get_trifinder(xy, t2v) 

44 else : 

45 self.trifinder = trifinder 

46 self.t2v = t2v 

47 t3xy = xy[t2v] 

48 self.map = MappingT3(t3xy) 

49 self.values = values 

50 assert(values.shape[0] == xy.shape[0]) 

51 self.fieldshape = values.shape[1:] 

52 

53 def __call__(self, xy): 

54 xy = np.atleast_2d(xy) 

55 elems = np.asarray(self.trifinder(xy), dtype =int) 

56 xy0 = self.map.xy[elems, 0, :] 

57 xy0s = xy - xy0 

58 uv = np.einsum('...jl,...l->...j', self.map.inv_F[elems], xy0s) 

59 valelems = self.values[self.t2v[elems]] 

60 u = uv[...,0] 

61 v = uv[...,1] 

62 shape = u.shape + tuple(1 for i in range(len(self.fieldshape))) 

63 u = u.reshape(shape) 

64 v = v.reshape(shape) 

65 nodeaxis = -1-len(self.fieldshape) 

66 valelemsv0 = valelems.take(0, axis = nodeaxis ) 

67 valelemsv1 = valelems.take(1, axis = nodeaxis ) 

68 valelemsv2 = valelems.take(2, axis = nodeaxis ) 

69 val = valelemsv0 + (valelemsv1-valelemsv0)*u + (valelemsv2-valelemsv0)*v 

70 if np.any(elems < 0) : 

71 warn("some point are outside of the mesh linearPerTriInterpolator," \ 

72 " corresponding values are set to np.inf") 

73 val[elems < 0] = np.inf 

74 return np.asarray(val, dtype = npfloat) 

75 

76if __name__ == '__main__': 

77 import gmsh 

78 import pyxfem.core.gmsh_mesh_generator as mg 

79 from pyxfem.core.mesh_simplex_one_level import gmshmodel2smesh 

80 import pyxfem.postpro.gmsh_post_pro_interface as gp 

81 gmsh.initialize() 

82 gmsh.option.setNumber("PostProcessing.AnimationCycle", 1) 

83 gmsh.option.setNumber('General.Terminal', 0) 

84 #model_name1, ____ = mg.square(length=1., relh=0.1, transfinite = True) 

85 model_name1, ____ = mg.crack_in_square_contour(relh = 0.1) 

86 gmsh.fltk.run() 

87 mesh1 = gmshmodel2smesh(model_name1) 

88 model_name2, ____ = mg.l_shape(relh=0.05) 

89 mesh2 = gmshmodel2smesh(model_name2) 

90 def fun_scal(xy): 

91 ''' a simple scalar function of coordinate ''' 

92 return xy[...,0]**2 

93 def fun_vec(xy): 

94 ''' a simple vector function of coordinate ''' 

95 return xy 

96 fun = fun_scal 

97 mesh1_trixy = mesh1.t2v_coord() 

98 gp.listPlotFieldTri(mesh1_trixy, fun(mesh1_trixy), 

99 viewname='base scalar field on mesh1', Range =[0.,0.25]) 

100 interp = LinearPerTriInterpolator(mesh1.xy, mesh1.t2v, fun(mesh1.xy)) 

101 val2 = interp(mesh2.xy) 

102 gp.listPlotFieldTri(mesh2.t2v_coord(), val2[np.asarray(mesh2.t2v, dtype=int)], 

103 viewname='interp scalar field on mesh2', Range =[0.,0.25]) 

104 

105 fun = fun_vec 

106 mesh1_trixy = mesh1.t2v_coord() 

107 gp.listPlotFieldTri(mesh1_trixy, fun(mesh1_trixy), 

108 viewname='base vector field on mesh1', Range = [0., np.sqrt(2)]) 

109 interp = LinearPerTriInterpolator(mesh1.xy, mesh1.t2v, fun(mesh1.xy), trifinder=interp.trifinder) 

110 val2 = interp(mesh2.xy) 

111 gp.listPlotFieldTri(mesh2.t2v_coord(), val2[np.asarray(mesh2.t2v, dtype=int)], 

112 viewname='interp vector field on mesh2', Range = [0., np.sqrt(2)]) 

113 gmsh.fltk.run()