Coverage for src/meshpy/core/nurbs_patch.py: 93%

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 module implements NURBS patches for the mesh.""" 

23 

24import numpy as _np 

25 

26from meshpy.core.conf import mpy as _mpy 

27from meshpy.core.element import Element as _Element 

28from meshpy.four_c.material import ( 

29 MaterialStVenantKirchhoff as _MaterialStVenantKirchhoff, 

30) 

31from meshpy.utils.environment import fourcipp_is_available as _fourcipp_is_available 

32 

33 

34class NURBSPatch(_Element): 

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

36 

37 # A list of valid material types for this element 

38 valid_material = [_MaterialStVenantKirchhoff] 

39 

40 def __init__( 

41 self, 

42 knot_vectors, 

43 polynomial_orders, 

44 element_string, 

45 material=None, 

46 nodes=None, 

47 element_description=None, 

48 ): 

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

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 # Set the element definitions 

61 self.element_string = element_string 

62 self.element_description = element_description 

63 

64 def get_nurbs_dimension(self): 

65 """Return the number of dimensions of the NURBS structure.""" 

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 dump_element_specific_section(self, yaml_dict): 

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

77 

78 if _fourcipp_is_available(): 

79 raise ValueError("Port this functionality to not use the legacy format.") 

80 

81 knotvectors_section = "STRUCTURE KNOTVECTORS" 

82 

83 if knotvectors_section not in yaml_dict.keys(): 

84 yaml_dict[knotvectors_section] = [] 

85 

86 section = yaml_dict[knotvectors_section] 

87 section.append(f"NURBS_DIMENSION {self.get_nurbs_dimension()}") 

88 section.append("BEGIN NURBSPATCH") 

89 section.append(f"ID {self.n_nurbs_patch}") 

90 

91 for dir_manifold in range(len(self.knot_vectors)): 

92 num_knotvectors = len(self.knot_vectors[dir_manifold]) 

93 section.append(f"NUMKNOTS {num_knotvectors}") 

94 section.append(f"DEGREE {self.polynomial_orders[dir_manifold]}") 

95 

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

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

98 # is an open (interpolated) knot vector. 

99 knotvector_type = "Interpolated" 

100 

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

102 if ( 

103 abs( 

104 self.knot_vectors[dir_manifold][i] 

105 - self.knot_vectors[dir_manifold][i + 1] 

106 ) 

107 > _mpy.eps_knot_vector 

108 ) or ( 

109 abs( 

110 self.knot_vectors[dir_manifold][num_knotvectors - 2 - i] 

111 - self.knot_vectors[dir_manifold][num_knotvectors - 1 - i] 

112 ) 

113 > _mpy.eps_knot_vector 

114 ): 

115 knotvector_type = "Periodic" 

116 break 

117 

118 section.append(f"TYPE {knotvector_type}") 

119 

120 for knot_vector_val in self.knot_vectors[dir_manifold]: 

121 section.append(knot_vector_val) 

122 

123 section.append("END NURBSPATCH") 

124 

125 def get_number_elements(self): 

126 """Get the number of elements in this patch by checking the amount of 

127 nonzero knot spans in the knot vector.""" 

128 

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

130 

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

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

133 if ( 

134 abs( 

135 self.knot_vectors[i_dir][i_knot] 

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

137 ) 

138 > _mpy.eps_knot_vector 

139 ): 

140 num_elements_dir[i_dir] += 1 

141 

142 total_num_elements = _np.prod(num_elements_dir) 

143 

144 return total_num_elements 

145 

146 def _check_material(self): 

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

148 element.""" 

149 for material_type in self.valid_material: 

150 if isinstance(self.material, material_type): 

151 return 

152 raise TypeError( 

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

154 f"ype {type(self.material)}!" 

155 ) 

156 

157 

158class NURBSSurface(NURBSPatch): 

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

160 

161 def __init__(self, *args, element_string=None, **kwargs): 

162 if element_string is None: 

163 element_string = "WALLNURBS" 

164 super().__init__(*args, element_string, **kwargs) 

165 

166 def dump_to_list(self): 

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

168 patch.""" 

169 

170 if _fourcipp_is_available(): 

171 raise ValueError("Port this functionality to not use the legacy format.") 

172 

173 # Check the material 

174 self._check_material() 

175 

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

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

178 

179 def get_ids_ctrlpts_surface(knot_span_u, knot_span_v): 

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

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

182 knot_span_v)""" 

183 

184 id_u = knot_span_u - self.polynomial_orders[0] 

185 id_v = knot_span_v - self.polynomial_orders[1] 

186 

187 element_ctrlpts_ids = [] 

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

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

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

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

192 # Cottrell, p. 314. 

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

194 element_ctrlpts_ids.append(index_global) 

195 

196 return element_ctrlpts_ids 

197 

198 patch_elements = [] 

199 

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

201 j = 0 

202 

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

204 for knot_span_v in range( 

205 self.polynomial_orders[1], 

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

207 ): 

208 for knot_span_u in range( 

209 self.polynomial_orders[0], 

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

211 ): 

212 element_cps_ids = get_ids_ctrlpts_surface(knot_span_u, knot_span_v) 

213 

214 string_cps = " ".join( 

215 [str(self.nodes[i].i_global) for i in element_cps_ids] 

216 ) 

217 

218 patch_elements.append( 

219 "{} {} NURBS{} {} MAT {} {}".format( 

220 self.i_global + j, 

221 self.element_string, 

222 (self.polynomial_orders[0] + 1) 

223 * (self.polynomial_orders[1] + 1), 

224 string_cps, 

225 self.material.i_global, 

226 self.element_description, 

227 ) 

228 ) 

229 j += 1 

230 

231 return patch_elements 

232 

233 

234class NURBSVolume(NURBSPatch): 

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

236 

237 def __init__(self, *args, element_string=None, **kwargs): 

238 if element_string is not None: 

239 raise ValueError("element_string is not yet implemented for NURBS volumes") 

240 super().__init__(*args, element_string, **kwargs) 

241 

242 def dump_to_list(self): 

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

244 patch.""" 

245 

246 if _fourcipp_is_available(): 

247 raise ValueError("Port this functionality to not use the legacy format.") 

248 

249 # Check the material 

250 self._check_material() 

251 

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

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

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

255 

256 def get_ids_ctrlpts_volume(knot_span_u, knot_span_v, knot_span_w): 

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

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

259 knot_span_v, knot_span_w)""" 

260 

261 id_u = knot_span_u - self.polynomial_orders[0] 

262 id_v = knot_span_v - self.polynomial_orders[1] 

263 id_w = knot_span_w - self.polynomial_orders[2] 

264 

265 element_ctrlpts_ids = [] 

266 

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

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

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

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

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

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

273 index_global = ( 

274 ctrlpts_size_u * ctrlpts_size_v * (id_w + k) 

275 + ctrlpts_size_u * (id_v + j) 

276 + id_u 

277 + i 

278 ) 

279 element_ctrlpts_ids.append(index_global) 

280 

281 return element_ctrlpts_ids 

282 

283 patch_elements = [] 

284 

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

286 increment_ele = 0 

287 

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

289 for knot_span_w in range( 

290 self.polynomial_orders[2], 

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

292 ): 

293 for knot_span_v in range( 

294 self.polynomial_orders[1], 

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

296 ): 

297 for knot_span_u in range( 

298 self.polynomial_orders[0], 

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

300 ): 

301 element_cps_ids = get_ids_ctrlpts_volume( 

302 knot_span_u, knot_span_v, knot_span_w 

303 ) 

304 

305 string_cps = " ".join( 

306 [str(self.nodes[i].i_global) for i in element_cps_ids] 

307 ) 

308 

309 num_cp_in_element = ( 

310 (self.polynomial_orders[0] + 1) 

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

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

313 ) 

314 

315 patch_elements.append( 

316 "{} SOLID NURBS{} {} MAT {} {}".format( 

317 self.i_global + increment_ele, 

318 num_cp_in_element, 

319 string_cps, 

320 self.material.i_global, 

321 self.element_description, 

322 ) 

323 ) 

324 increment_ele += 1 

325 

326 return patch_elements