Coverage for src/meshpy/four_c/solid_shell_thickness_direction.py: 87%

1# The MIT License (MIT) 

2# 

3# Copyright (c) 2018-2025 MeshPy Authors 

4# 

5# Permission is hereby granted, free of charge, to any person obtaining a copy 

6# of this software and associated documentation files (the "Software"), to deal 

7# in the Software without restriction, including without limitation the rights 

8# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 

9# copies of the Software, and to permit persons to whom the Software is 

10# furnished to do so, subject to the following conditions: 

11# 

12# The above copyright notice and this permission notice shall be included in 

13# all copies or substantial portions of the Software. 

14# 

15# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 

16# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

17# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 

18# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 

19# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 

20# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 

21# THE SOFTWARE. 

22"""This function allows to reorder the connectivity of solid shell elements 

23such that the solid shell direction is correctly represented.""" 

24 

25from typing import List as _List 

26 

27import numpy as _np 

28import pyvista as _pv 

29 

30from meshpy.core.element import Element as _Element 

31from meshpy.core.element_volume import VolumeElement as _VolumeElement 

32from meshpy.core.element_volume import VolumeHEX8 as _VolumeHEX8 

33from meshpy.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates 

34 

35 

36def shape_functions_hex8(xi1, xi2, xi3): 

37 """Return the shape functions for a hex8 element.""" 

38 shape_functions = _np.zeros(8) 

39 one_over_eight = 0.125 

40 shape_functions[0] = one_over_eight * (1 - xi1) * (1 - xi2) * (1 - xi3) 

41 shape_functions[1] = one_over_eight * (1 + xi1) * (1 - xi2) * (1 - xi3) 

42 shape_functions[2] = one_over_eight * (1 + xi1) * (1 + xi2) * (1 - xi3) 

43 shape_functions[3] = one_over_eight * (1 - xi1) * (1 + xi2) * (1 - xi3) 

44 shape_functions[4] = one_over_eight * (1 - xi1) * (1 - xi2) * (1 + xi3) 

45 shape_functions[5] = one_over_eight * (1 + xi1) * (1 - xi2) * (1 + xi3) 

46 shape_functions[6] = one_over_eight * (1 + xi1) * (1 + xi2) * (1 + xi3) 

47 shape_functions[7] = one_over_eight * (1 - xi1) * (1 + xi2) * (1 + xi3) 

48 return shape_functions 

49 

50 

51def shape_functions_derivative_hex8(xi1, xi2, xi3): 

52 """Return the derivative of the shape functions for a hex8 element.""" 

53 derivatives = _np.zeros((3, 8)) 

54 one_over_eight = 0.125 

55 

56 # Derivatives with respect to xi1 

57 derivatives[0, 0] = -one_over_eight * (1 - xi2) * (1 - xi3) 

58 derivatives[0, 1] = one_over_eight * (1 - xi2) * (1 - xi3) 

59 derivatives[0, 2] = one_over_eight * (1 + xi2) * (1 - xi3) 

60 derivatives[0, 3] = -one_over_eight * (1 + xi2) * (1 - xi3) 

61 derivatives[0, 4] = -one_over_eight * (1 - xi2) * (1 + xi3) 

62 derivatives[0, 5] = one_over_eight * (1 - xi2) * (1 + xi3) 

63 derivatives[0, 6] = one_over_eight * (1 + xi2) * (1 + xi3) 

64 derivatives[0, 7] = -one_over_eight * (1 + xi2) * (1 + xi3) 

65 

66 # Derivatives with respect to xi2 

67 derivatives[1, 0] = -one_over_eight * (1 - xi1) * (1 - xi3) 

68 derivatives[1, 1] = -one_over_eight * (1 + xi1) * (1 - xi3) 

69 derivatives[1, 2] = one_over_eight * (1 + xi1) * (1 - xi3) 

70 derivatives[1, 3] = one_over_eight * (1 - xi1) * (1 - xi3) 

71 derivatives[1, 4] = -one_over_eight * (1 - xi1) * (1 + xi3) 

72 derivatives[1, 5] = -one_over_eight * (1 + xi1) * (1 + xi3) 

73 derivatives[1, 6] = one_over_eight * (1 + xi1) * (1 + xi3) 

74 derivatives[1, 7] = one_over_eight * (1 - xi1) * (1 + xi3) 

75 

76 # Derivatives with respect to xi3 

77 derivatives[2, 0] = -one_over_eight * (1 - xi1) * (1 - xi2) 

78 derivatives[2, 1] = -one_over_eight * (1 + xi1) * (1 - xi2) 

79 derivatives[2, 2] = -one_over_eight * (1 + xi1) * (1 + xi2) 

80 derivatives[2, 3] = -one_over_eight * (1 - xi1) * (1 + xi2) 

81 derivatives[2, 4] = one_over_eight * (1 - xi1) * (1 - xi2) 

82 derivatives[2, 5] = one_over_eight * (1 + xi1) * (1 - xi2) 

83 derivatives[2, 6] = one_over_eight * (1 + xi1) * (1 + xi2) 

84 derivatives[2, 7] = one_over_eight * (1 - xi1) * (1 + xi2) 

85 

86 return derivatives 

87 

88 

89def get_hex8_element_center_and_jacobian_mapping(element): 

90 """Return the center of a hex8 element and the Jacobian mapping for that 

91 point.""" 

92 

93 nodal_coordinates = _get_nodal_coordinates(element.nodes) 

94 if not len(nodal_coordinates) == 8: 

95 raise ValueError(f"Expected 8 nodes, got {len(nodal_coordinates)}") 

96 

97 N = shape_functions_hex8(0, 0, 0) 

98 dN = shape_functions_derivative_hex8(0, 0, 0) 

99 

100 reference_position_center = _np.dot(N, nodal_coordinates) 

101 jacobian_center = _np.dot(dN, nodal_coordinates) 

102 

103 return reference_position_center, jacobian_center 

104 

105 

106def get_reordering_index_thickness(jacobian, *, identify_threshold=None): 

107 """Return the reordering index from the Jacobian such that the thinnest 

108 direction is the 3rd parameter direction. 

109 

110 Additionally it is checked, that the thinnest direction is at least 

111 identify_threshold times smaller than the next thinnest, to avoid 

112 wrongly detected directions. 

113 """ 

114 

115 # The direction with the smallest parameter derivative is the thickness direction 

116 parameter_derivative_norms = [ 

117 _np.linalg.norm(parameter_direction) for parameter_direction in jacobian 

118 ] 

119 thickness_direction = _np.argmin(parameter_derivative_norms) 

120 

121 if identify_threshold is not None: 

122 # Check that the minimal parameter direction is at least a given factor 

123 # smaller than the other directions. This helps to identify cases where 

124 # it is unlikely that a unique direction is found. 

125 min_norm = parameter_derivative_norms[thickness_direction] 

126 relative_difference = [ 

127 parameter_derivative_norms[i] / min_norm 

128 for i in range(3) 

129 if not i == thickness_direction 

130 ] 

131 if _np.min(relative_difference) < 1.5: 

132 raise ValueError("Could not uniquely identify the thickness direction.") 

133 

134 return thickness_direction 

135 

136 

137def get_reordering_index_director_projection( 

138 jacobian, director, *, identify_threshold=None 

139): 

140 """Return the reordering index from the Jacobian such that the thickness 

141 direction is the one that has the largest dot product with the given 

142 director.""" 

143 

144 projections = [] 

145 for parameter_director in jacobian: 

146 parameter_director = parameter_director / _np.linalg.norm(parameter_director) 

147 projections.append(_np.abs(_np.dot(director, parameter_director))) 

148 thickness_direction = _np.argmax(projections) 

149 

150 if identify_threshold is not None: 

151 # Check that the maximal dot product is at least a given factor larger than 

152 # the other projections. This helps to identify cases where it is unlikely 

153 # that a unique direction is found. 

154 max_norm = projections[thickness_direction] 

155 relative_difference = [ 

156 projections[i] / max_norm for i in range(3) if not i == thickness_direction 

157 ] 

158 if _np.max(relative_difference) > 1.0 / identify_threshold: 

159 raise ValueError("Could not uniquely identify the thickness direction.") 

160 

161 return thickness_direction 

162 

163 

164def set_solid_shell_thickness_direction( 

165 elements: _List[_Element], 

166 *, 

167 selection_type="thickness", 

168 director=None, 

169 director_function=None, 

170 identify_threshold=2.0, 

171): 

172 """Set the solid shell directions for all solid shell elements in the 

173 element list. 

174 

175 Args: 

176 ---- 

177 elements: 

178 A list containing all elements that should be checked 

179 selection_type: 

180 The type of algorithm that shall be used to select the thickness direction 

181 "thickness": 

182 The "smallest" dimension of the element will be set to the 

183 thickness direction 

184 "projection_director": 

185 The parameter director that aligns most with a given director 

186 will be set as the thickness direction 

187 "projection_director_function": 

188 The parameter director that aligns most with a director obtained 

189 by a given function of the element centroid coordinates will be 

190 set as the thickness direction 

191 identify_threshold: float/None 

192 To ensure that the found directions are well-defined, i.e., that not multiple 

193 directions are almost equally suited to the thickness direction. 

194 """ 

195 

196 if len(elements) == 0: 

197 raise ValueError("Expected a non empty element list") 

198 

199 for element in elements: 

200 is_hex8 = isinstance(element, _VolumeHEX8) 

201 if is_hex8: 

202 is_solid_shell = "SOLIDSH8" in element.string_pre_nodes # type: ignore[attr-defined] 

203 

204 if is_solid_shell: 

205 # Get the element center and the Jacobian at the center 

206 ( 

207 reference_position_center, 

208 jacobian_center, 

209 ) = get_hex8_element_center_and_jacobian_mapping(element) 

210 

211 # Depending on the chosen method, get the thickness direction 

212 if selection_type == "thickness": 

213 thickness_direction = get_reordering_index_thickness( 

214 jacobian_center, identify_threshold=identify_threshold 

215 ) 

216 elif selection_type == "projection_director": 

217 thickness_direction = get_reordering_index_director_projection( 

218 jacobian_center, director, identify_threshold=identify_threshold 

219 ) 

220 elif selection_type == "projection_director_function": 

221 director = director_function(reference_position_center) 

222 thickness_direction = get_reordering_index_director_projection( 

223 jacobian_center, director, identify_threshold=identify_threshold 

224 ) 

225 else: 

226 raise ValueError( 

227 f'Got unexpected selection_type of value "{selection_type}"' 

228 ) 

229 

230 n_apply_mapping = 0 

231 if thickness_direction == 2: 

232 # We already have the orientation we want 

233 continue 

234 elif thickness_direction == 1: 

235 # We need to apply the connectivity mapping once, i.e., the 2nd parameter 

236 # direction has to become the 3rd 

237 n_apply_mapping = 1 

238 elif thickness_direction == 0: 

239 # We need to apply the connectivity mapping twice, i.e., the 1nd parameter 

240 # direction has to become the 3rd 

241 n_apply_mapping = 2 

242 

243 # This permutes the parameter coordinate the following way: 

244 # [xi,eta,zeta]->[zeta,xi,eta] 

245 mapping = [2, 6, 7, 3, 1, 5, 4, 0] 

246 for _ in range(n_apply_mapping): 

247 element.nodes = [ 

248 element.nodes[local_index] for local_index in mapping 

249 ] 

250 

251 

252def get_visualization_third_parameter_direction_hex8(mesh): 

253 """Return a pyvista mesh with cell data for the third parameter direction 

254 for hex8 elements.""" 

255 

256 vtk_solid = mesh.get_vtk_representation()[1].grid 

257 pv_solid = _pv.UnstructuredGrid(vtk_solid) 

258 

259 cell_thickness_direction = [] 

260 for element in mesh.elements: 

261 if isinstance(element, _VolumeHEX8): 

262 _, jacobian_center = get_hex8_element_center_and_jacobian_mapping(element) 

263 cell_thickness_direction.append(jacobian_center[2]) 

264 elif isinstance(element, _VolumeElement): 

265 cell_thickness_direction.append([0, 0, 0]) 

266 

267 if not len(cell_thickness_direction) == pv_solid.number_of_cells: 

268 raise ValueError( 

269 "Expected the same number of cells from the mesh and from the " 

270 f"pyvista object. Got {len(cell_thickness_direction)} form the " 

271 f"mesh and {pv_solid.number_of_cells} form pyvista" 

272 ) 

273 

274 pv_solid.cell_data["thickness_direction"] = cell_thickness_direction 

275 return pv_solid 

276 

277 

278def visualize_third_parameter_direction_hex8(mesh): 

279 """Visualize the third parameter direction for hex8 elements. 

280 

281 This can be used to check the correct definition of the shell 

282 thickness for solid shell elements. 

283 """ 

284 

285 grid = get_visualization_third_parameter_direction_hex8(mesh) 

286 grid = grid.clean() 

287 cell_centers = grid.cell_centers() 

288 thickness_direction = cell_centers.glyph( 

289 orient="thickness_direction", scale="thickness_direction", factor=5 

290 ) 

291 

292 plotter = _pv.Plotter() 

293 plotter.renderer.add_axes() 

294 plotter.add_mesh(grid, color="white", show_edges=True, opacity=0.5) 

295 plotter.add_mesh(thickness_direction, color="red") 

296 plotter.show()