Coverage for /usr/local/lib/python3.8/site-packages/spam/mesh/structured.py: 100%

1# Library of SPAM functions for dealing with structured meshes 

2# Copyright (C) 2020 SPAM Contributors 

3# 

4# This program is free software: you can redistribute it and/or modify it 

5# under the terms of the GNU General Public License as published by the Free 

6# Software Foundation, either version 3 of the License, or (at your option) 

7# any later version. 

8# 

9# This program is distributed in the hope that it will be useful, but WITHOUT 

10# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

11# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

12# more details. 

13# 

14# You should have received a copy of the GNU General Public License along with 

15# this program. If not, see <http://www.gnu.org/licenses/>. 

16 

17""" 

18This module offers a set of tools in order to manipulate structured meshes. 

19 

20>>> # import module 

21>>> import spam.mesh 

22 

23 

24The strucutred VTK files used to save the data have the form: 

25 

26.. code-block:: text 

27 

28 # vtk DataFile Version 2.0 

29 VTK file from spam: spam.vtk 

30 ASCII 

31 

32 DATASET STRUCTURED_POINTS 

33 DIMENSIONS nx ny nz 

34 ASPECT_RATIO lx ly lz 

35 ORIGIN ox oy oz 

36 

37 POINT_DATA nx x ny x nz 

38 

39 SCALARS myNodalField1 float 

40 LOOKUP_TABLE default 

41 nodalValue_1 

42 nodalValue_2 

43 nodalValue_3 

44 ... 

45 

46 VECTORS myNodalField2 float 

47 LOOKUP_TABLE default 

48 nodalValue_1_X nodalValue_1_Y nodalValue_1_Z 

49 nodalValue_2_X nodalValue_2_Y nodalValue_2_Z 

50 nodalValue_3_X nodalValue_3_Y nodalValue_3_Z 

51 ... 

52 

53 CELL_DATA (nx-1) x (ny-1) x (nz-1) 

54 

55 SCALARS myCellField1 float 

56 LOOKUP_TABLE default 

57 cellValue_1 

58 cellValue_2 

59 cellValue_3 

60 ... 

61 

62where nx, ny and nz are the number of nodes in each axis, lx, ly, lz, the mesh length in each axis and ox, oy, oz the spatial postiion of the origin. 

63""" 

64 

65 

66def createCylindricalMask(shape, radius, voxSize=1.0, centre=None): 

67 """ 

68 Create a image mask of a cylinder in the z direction. 

69 

70 Parameters 

71 ---------- 

72 shape: array, int 

73 The shape of the array the where the cylinder is saved 

74 

75 radius: float 

76 The radius of the cylinder 

77 

78 voxSize: float (default=1.0) 

79 The physical size of a voxel 

80 

81 centre: array of floats of size 2, (default None) 

82 The center [y,x] of the axis of rotation of the cylinder. 

83 If None it is taken to be the centre of the array. 

84 

85 Returns 

86 ------- 

87 cyl: array, bool 

88 The cylinder 

89 

90 """ 

91 import numpy 

92 

93 cyl = numpy.zeros(shape).astype(bool) 

94 

95 if centre is None: 

96 centre = [float(shape[1]) / 2.0, float(shape[2]) / 2.0] 

97 

98 for iy in range(cyl.shape[1]): 

99 y = (float(iy) + 0.5) * float(voxSize) 

100 for ix in range(cyl.shape[2]): 

101 x = (float(ix) + 0.5) * float(voxSize) 

102 dist = numpy.sqrt((x - centre[1]) ** 2 + (y - centre[0]) ** 2) 

103 if dist < radius: 

104 cyl[:, iy, ix] = True 

105 

106 return cyl 

107 

108 

109def structuringElement(radius=1, order=2, margin=0, dim=3): 

110 """ 

111 This function construct a structural element. 

112 

113 Parameters 

114 ----------- 

115 radius : int, default=1 

116 The `radius` of the structural element 

117 

118 .. code-block:: text 

119 

120 radius = 1 gives 3x3x3 arrays 

121 radius = 2 gives 5x5x5 arrays 

122 ... 

123 radius = n gives (2n+1)x(2n+1)x(2n+1) arrays 

124 

125 order : int, default=2 

126 Defines the shape of the structuring element by setting the order of the norm 

127 used to compute the distance between the centre and the border. 

128 

129 A representation for the slices of a 5x5x5 element (size=2) from the center to on corner (1/8 of the cube) 

130 

131 .. code-block:: text 

132 

133 order=numpy.inf: the cube 

134 1 1 1 1 1 1 1 1 1 

135 1 1 1 1 1 1 1 1 1 

136 1 1 1 1 1 1 1 1 1 

137 

138 order=2: the sphere 

139 1 0 0 0 0 0 0 0 0 

140 1 1 0 1 1 0 0 0 0 

141 1 1 1 1 1 0 1 0 0 

142 

143 order=1: the diamond 

144 1 0 0 0 0 0 0 0 0 

145 1 1 0 1 0 0 0 0 0 

146 1 1 1 1 1 0 1 0 0 

147 

148 margin : int, default=0 

149 Gives a 0 valued margin of size margin. 

150 

151 dim : int, default=3 

152 Spatial dimension (2 or 3). 

153 

154 Returns 

155 -------- 

156 array 

157 The structural element 

158 """ 

159 import numpy 

160 

161 tb = tuple([2 * radius + 2 * margin + 1 for _ in range(dim)]) 

162 ts = tuple([2 * radius + 1 for _ in range(dim)]) 

163 c = numpy.abs(numpy.indices(ts) - radius) 

164 d = numpy.zeros(tb) 

165 s = tuple([slice(margin, margin + 2 * radius + 1) for _ in range(dim)]) 

166 d[s] = numpy.power(numpy.sum(numpy.power(c, order), axis=0), 1.0 / float(order)) <= radius 

167 return d.astype("<u1") 

168 

169 

170def createLexicoCoordinates(lengths, nNodes, origin=(0, 0, 0)): 

171 """ 

172 Create a list of coordinates following the lexicographical order. 

173 

174 Parameters 

175 ---------- 

176 lengths : array of floats 

177 The length of the cuboids in every directions. 

178 nNodes : array of int 

179 The number of nodes of the mesh in every directions. 

180 origin : array of floats 

181 The coordinates of the origin of the mesh. 

182 

183 Returns 

184 ------- 

185 array 

186 The list of coordinates. ``shape=(nx*ny*nz, 3)`` 

187 

188 """ 

189 import numpy 

190 

191 x = numpy.linspace(origin[0], lengths[0] + origin[0], nNodes[0]) 

192 y = numpy.linspace(origin[1], lengths[1] + origin[1], nNodes[1]) 

193 z = numpy.linspace(origin[2], lengths[2] + origin[2], nNodes[2]) 

194 cx = numpy.tile(x, (1, nNodes[1] * nNodes[2])) 

195 cy = numpy.tile(numpy.sort(numpy.tile(y, (1, nNodes[0]))), (1, nNodes[2])) 

196 cz = numpy.sort(numpy.tile(z, (1, nNodes[0] * nNodes[1]))) 

197 return numpy.transpose([cx[0], cy[0], cz[0]])