Coverage for src/beamme/core/nurbs_patch.py: 96%

1# The MIT License (MIT) 

2# 

3# Copyright (c) 2018-2025 BeamMe 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 module implements NURBS patches for the mesh.""" 

23 

24from typing import Any as _Any 

25 

26import numpy as _np 

27 

28from beamme.core.conf import mpy as _mpy 

29from beamme.core.element import Element as _Element 

30from beamme.core.material import ( 

31 MaterialSolidBase as _MaterialSolidBase, 

32) 

33 

34 

35class NURBSPatch(_Element): 

36 """A base class for a NURBS patch.""" 

37 

38 # A list of valid material types for this element 

39 valid_material = [_MaterialSolidBase] 

40 

41 def __init__( 

42 self, 

43 knot_vectors, 

44 polynomial_orders, 

45 material=None, 

46 nodes=None, 

47 data=None, 

48 ): 

49 super().__init__(nodes=nodes, material=material, data=data) 

50 

51 # Knot vectors 

52 self.knot_vectors = knot_vectors 

53 

54 # Polynomial degrees 

55 self.polynomial_orders = polynomial_orders 

56 

57 # Set numbers for elements 

58 self.n_nurbs_patch = None 

59 

60 def get_nurbs_dimension(self) -> int: 

61 """Determine the number of dimensions of the NURBS structure. 

62 

63 Returns: 

64 Number of dimensions of the NURBS object. 

65 """ 

66 n_knots = len(self.knot_vectors) 

67 n_polynomial = len(self.polynomial_orders) 

68 if not n_knots == n_polynomial: 

69 raise ValueError( 

70 "The variables n_knots and polynomial_orders should have " 

71 f"the same length. Got {n_knots} and {n_polynomial}" 

72 ) 

73 return n_knots 

74 

75 def get_number_of_control_points_per_dir(self) -> list[int]: 

76 """Determine the number of control points in each parameter direction 

77 of the patch. 

78 

79 Returns: 

80 List of control points per direction. 

81 """ 

82 n_dim = len(self.knot_vectors) 

83 n_cp_per_dim = [] 

84 for i_dim in range(n_dim): 

85 knot_vector_size = len(self.knot_vectors[i_dim]) 

86 polynomial_order = self.polynomial_orders[i_dim] 

87 n_cp_per_dim.append(knot_vector_size - polynomial_order - 1) 

88 return n_cp_per_dim 

89 

90 def dump_element_specific_section(self, input_file) -> None: 

91 """Set the knot vectors of the NURBS patch in the input file.""" 

92 

93 patch_data: dict[str, _Any] = { 

94 "knot_vectors": [], 

95 } 

96 

97 for dir_manifold in range(self.get_nurbs_dimension()): 

98 knotvector = self.knot_vectors[dir_manifold] 

99 num_knots = len(knotvector) 

100 

101 # Check the type of knot vector, in case that the multiplicity of the first and last 

102 # knot vectors is not p + 1, then it is a closed (periodic) knot vector, otherwise it 

103 # is an open (interpolated) knot vector. 

104 knotvector_type = "Interpolated" 

105 

106 for i in range(self.polynomial_orders[dir_manifold] - 1): 

107 if (abs(knotvector[i] - knotvector[i + 1]) > _mpy.eps_knot_vector) or ( 

108 abs(knotvector[num_knots - 2 - i] - knotvector[num_knots - 1 - i]) 

109 > _mpy.eps_knot_vector 

110 ): 

111 knotvector_type = "Periodic" 

112 break 

113 

114 patch_data["knot_vectors"].append( 

115 { 

116 "DEGREE": self.polynomial_orders[dir_manifold], 

117 "TYPE": knotvector_type, 

118 "knots": [ 

119 knot_vector_val 

120 for knot_vector_val in self.knot_vectors[dir_manifold] 

121 ], 

122 } 

123 ) 

124 

125 if "STRUCTURE KNOTVECTORS" not in input_file: 

126 input_file.add({"STRUCTURE KNOTVECTORS": []}) 

127 input_file["STRUCTURE KNOTVECTORS"] = [] 

128 

129 patches = input_file["STRUCTURE KNOTVECTORS"] 

130 patch_data["ID"] = len(patches) + 1 

131 patches.append(patch_data) 

132 

133 def get_number_elements(self) -> int: 

134 """Determine the number of elements in this patch by checking the 

135 amount of nonzero knot spans in the knot vector. 

136 

137 Returns: 

138 Number of elements for this patch. 

139 """ 

140 

141 num_elements_dir = _np.zeros(len(self.knot_vectors), dtype=int) 

142 

143 for i_dir in range(len(self.knot_vectors)): 

144 for i_knot in range(len(self.knot_vectors[i_dir]) - 1): 

145 if ( 

146 abs( 

147 self.knot_vectors[i_dir][i_knot] 

148 - self.knot_vectors[i_dir][i_knot + 1] 

149 ) 

150 > _mpy.eps_knot_vector 

151 ): 

152 num_elements_dir[i_dir] += 1 

153 

154 total_num_elements = _np.prod(num_elements_dir) 

155 

156 return total_num_elements 

157 

158 def _check_material(self) -> None: 

159 """Check if the linked material is valid for this type of NURBS solid 

160 element.""" 

161 for material_type in type(self).valid_material: 

162 if isinstance(self.material, material_type): 

163 return 

164 raise TypeError( 

165 f"NURBS solid of type {type(self)} can not have a material of" 

166 f" type {type(self.material)}!" 

167 ) 

168 

169 

170class NURBSSurface(NURBSPatch): 

171 """A patch of a NURBS surface.""" 

172 

173 def __init__(self, *args, **kwargs): 

174 super().__init__(*args, **kwargs) 

175 

176 def dump_to_list(self): 

177 """Return a list with all the element definitions contained in this 

178 patch.""" 

179 

180 # Check the material 

181 self._check_material() 

182 

183 # Calculate how many control points are on the u direction 

184 ctrlpts_size_u = len(self.knot_vectors[0]) - self.polynomial_orders[0] - 1 

185 

186 def get_ids_ctrlpts_surface(knot_span_u, knot_span_v): 

187 """For an interpolated patch, calculate control points involved in 

188 evaluation of the surface point at the knot span (knot_span_u, 

189 knot_span_v)""" 

190 

191 id_u = knot_span_u - self.polynomial_orders[0] 

192 id_v = knot_span_v - self.polynomial_orders[1] 

193 

194 element_ctrlpts_ids = [] 

195 for j in range(self.polynomial_orders[1] + 1): 

196 for i in range(self.polynomial_orders[0] + 1): 

197 # Calculating the global index of the control point, based on the book 

198 # "Isogeometric Analysis: toward Integration of CAD and FEA" of J. Austin 

199 # Cottrell, p. 314. 

200 index_global = ctrlpts_size_u * (id_v + j) + id_u + i 

201 element_ctrlpts_ids.append(index_global) 

202 

203 return element_ctrlpts_ids 

204 

205 patch_elements = [] 

206 

207 # Adding an increment j to self.global to obtain the ID of an element in the patch 

208 j = 0 

209 

210 # Loop over the knot spans to obtain the elements inside the patch 

211 for knot_span_v in range( 

212 self.polynomial_orders[1], 

213 len(self.knot_vectors[1]) - self.polynomial_orders[1] - 1, 

214 ): 

215 for knot_span_u in range( 

216 self.polynomial_orders[0], 

217 len(self.knot_vectors[0]) - self.polynomial_orders[0] - 1, 

218 ): 

219 element_cps_ids = get_ids_ctrlpts_surface(knot_span_u, knot_span_v) 

220 

221 connectivity = [self.nodes[i] for i in element_cps_ids] 

222 

223 num_cp_in_element = (self.polynomial_orders[0] + 1) * ( 

224 self.polynomial_orders[1] + 1 

225 ) 

226 

227 patch_elements.append( 

228 { 

229 "id": self.i_global + j, 

230 "cell": { 

231 "type": f"NURBS{num_cp_in_element}", 

232 "connectivity": connectivity, 

233 }, 

234 "data": { 

235 "type": "WALLNURBS", 

236 "MAT": self.material, 

237 **(self.data if self.data else {}), 

238 }, 

239 } 

240 ) 

241 j += 1 

242 

243 return patch_elements 

244 

245 

246class NURBSVolume(NURBSPatch): 

247 """A patch of a NURBS volume.""" 

248 

249 def __init__(self, *args, **kwargs): 

250 super().__init__(*args, **kwargs) 

251 

252 def dump_to_list(self): 

253 """Return a list with all the element definitions contained in this 

254 patch.""" 

255 

256 # Check the material 

257 self._check_material() 

258 

259 # Calculate how many control points are on the u and v directions 

260 ctrlpts_size_u = len(self.knot_vectors[0]) - self.polynomial_orders[0] - 1 

261 ctrlpts_size_v = len(self.knot_vectors[1]) - self.polynomial_orders[1] - 1 

262 

263 def get_ids_ctrlpts_volume(knot_span_u, knot_span_v, knot_span_w): 

264 """For an interpolated patch, calculate control points involved in 

265 evaluation of the surface point at the knot span (knot_span_u, 

266 knot_span_v, knot_span_w)""" 

267 

268 id_u = knot_span_u - self.polynomial_orders[0] 

269 id_v = knot_span_v - self.polynomial_orders[1] 

270 id_w = knot_span_w - self.polynomial_orders[2] 

271 

272 element_ctrlpts_ids = [] 

273 

274 for k in range(self.polynomial_orders[2] + 1): 

275 for j in range(self.polynomial_orders[1] + 1): 

276 for i in range(self.polynomial_orders[0] + 1): 

277 # Calculating the global index of the control point, based on the paper 

278 # "Isogeometric analysis: an overview and computer implementation aspects" 

279 # of Vinh-Phu Nguyen, Mathematics and Computers in Simulation, Jun-2015. 

280 index_global = ( 

281 ctrlpts_size_u * ctrlpts_size_v * (id_w + k) 

282 + ctrlpts_size_u * (id_v + j) 

283 + id_u 

284 + i 

285 ) 

286 element_ctrlpts_ids.append(index_global) 

287 

288 return element_ctrlpts_ids 

289 

290 patch_elements = [] 

291 

292 # Adding an increment to self.global to obtain the ID of an element in the patch 

293 increment_ele = 0 

294 

295 # Loop over the knot spans to obtain the elements inside the patch 

296 for knot_span_w in range( 

297 self.polynomial_orders[2], 

298 len(self.knot_vectors[2]) - self.polynomial_orders[2] - 1, 

299 ): 

300 for knot_span_v in range( 

301 self.polynomial_orders[1], 

302 len(self.knot_vectors[1]) - self.polynomial_orders[1] - 1, 

303 ): 

304 for knot_span_u in range( 

305 self.polynomial_orders[0], 

306 len(self.knot_vectors[0]) - self.polynomial_orders[0] - 1, 

307 ): 

308 element_cps_ids = get_ids_ctrlpts_volume( 

309 knot_span_u, knot_span_v, knot_span_w 

310 ) 

311 

312 connectivity = [self.nodes[i] for i in element_cps_ids] 

313 

314 num_cp_in_element = ( 

315 (self.polynomial_orders[0] + 1) 

316 * (self.polynomial_orders[1] + 1) 

317 * (self.polynomial_orders[2] + 1) 

318 ) 

319 

320 patch_elements.append( 

321 { 

322 "id": self.i_global + increment_ele, 

323 "cell": { 

324 "type": f"NURBS{num_cp_in_element}", 

325 "connectivity": connectivity, 

326 }, 

327 "data": { 

328 "type": "SOLID", 

329 "MAT": self.material, 

330 **(self.data if self.data else {}), 

331 }, 

332 } 

333 ) 

334 increment_ele += 1 

335 

336 return patch_elements