Coverage for /usr/local/lib/python3.10/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""" 

64import numpy 

65 

66 

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

68 """ 

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

70 

71 Parameters 

72 ---------- 

73 shape: array, int 

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

75 

76 radius: float 

77 The radius of the cylinder 

78 

79 voxSize: float (default=1.0) 

80 The physical size of a voxel 

81 

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

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

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

85 

86 Returns 

87 ------- 

88 cyl: array, bool 

89 The cylinder 

90 

91 """ 

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

93 

94 if centre is None: 

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

96 

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

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

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

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

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

102 if dist < radius: 

103 cyl[:, iy, ix] = True 

104 

105 return cyl 

106 

107 

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

109 """ 

110 This function construct a structural element. 

111 

112 Parameters 

113 ----------- 

114 radius : int, default=1 

115 The `radius` of the structural element 

116 

117 .. code-block:: text 

118 

119 radius = 1 gives 3x3x3 arrays 

120 radius = 2 gives 5x5x5 arrays 

121 ... 

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

123 

124 order : int, default=2 

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

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

127 

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

129 

130 .. code-block:: text 

131 

132 order=numpy.inf: the cube 

133 1 1 1 1 1 1 1 1 1 

134 1 1 1 1 1 1 1 1 1 

135 1 1 1 1 1 1 1 1 1 

136 

137 order=2: the sphere 

138 1 0 0 0 0 0 0 0 0 

139 1 1 0 1 1 0 0 0 0 

140 1 1 1 1 1 0 1 0 0 

141 

142 order=1: the diamond 

143 1 0 0 0 0 0 0 0 0 

144 1 1 0 1 0 0 0 0 0 

145 1 1 1 1 1 0 1 0 0 

146 

147 margin : int, default=0 

148 Gives a 0 valued margin of size margin. 

149 

150 dim : int, default=3 

151 Spatial dimension (2 or 3). 

152 

153 Returns 

154 -------- 

155 array 

156 The structural element 

157 """ 

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

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

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

161 d = numpy.zeros(tb) 

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

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

164 return d.astype("<u1") 

165 

166 

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

168 """ 

169 Create a list of coordinates following the lexicographical order. 

170 

171 Parameters 

172 ---------- 

173 lengths : array of floats 

174 The length of the cuboids in every directions. 

175 nNodes : array of int 

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

177 origin : array of floats 

178 The coordinates of the origin of the mesh. 

179 

180 Returns 

181 ------- 

182 array 

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

184 

185 """ 

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

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

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

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

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

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

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