Coverage for /opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/hypervehicle/components/component.py: 48%

1import pymeshfix

2import numpy as np

3from stl import mesh

4from copy import deepcopy

5from abc import abstractmethod

6from collections import Counter

7from multiprocess.pool import Pool

8from typing import Callable, Union, Optional

9from hypervehicle.geometry.vector import Vector3

10from hypervehicle.geometry.surface import StructuredGrid, ParametricSurface

11from hypervehicle.utilities import surfce_to_stl

12from hypervehicle.geometry.geometry import (

13 CurvedPatch,

14 RotatedPatch,

15 MirroredPatch,

16 OffsetPatchFunction,

17)

20class AbstractComponent:

21 componenttype = None

23 @abstractmethod

24 def __init__(self, params: dict, verbosity: int = 1) -> None:

25 pass

27 @abstractmethod

28 def __repr__(self):

29 pass

31 @abstractmethod

32 def __str__(self):

33 pass

35 @property

36 @abstractmethod

37 def componenttype(self):

38 # This is a placeholder for a class variable defining the component type

39 pass

41 @abstractmethod

42 def generate_patches(self):

43 """Generates the parametric patches from the parameter dictionary."""

44 pass

46 @abstractmethod

47 def curve(self):

48 """Applies a curvature function to the parametric patches."""

49 pass

51 @abstractmethod

52 def rotate(self, angle: float = 0, axis: str = "y"):

53 """Rotates the parametric patches."""

54 pass

56 @abstractmethod

57 def reflect(self):

58 """Reflects the parametric patches."""

59 pass

61 @abstractmethod

62 def grid(self):

63 """Creates a discrete grid from the parametric patches."""

64 pass

66 @abstractmethod

67 def surface(self):

68 """Creates the discretised surface data from the

69 parametric patches."""

70 pass

72 @abstractmethod

73 def to_vtk(self):

74 """Writes the component to VTK file format."""

75 pass

77 @abstractmethod

78 def to_stl(self):

79 """Writes the component to STL file format."""

80 pass

82 @abstractmethod

83 def analyse(self):

84 """Evaluates properties of the STL mesh."""

85 pass

88class Component(AbstractComponent):

89 def __init__(

90 self,

91 params: Optional[dict] = None,

92 edges: Optional[list] = None,

93 stl_resolution: Optional[int] = 2,

94 verbosity: Optional[int] = 1,

95 name: Optional[str] = None,

96 ) -> None:

97 # Set verbosity

98 self.verbosity = verbosity

100 # Save parameters

101 self.params = params

102

103 # Save edges

104 self.edges = edges if edges is not None else []

105

106 # Processed objects

107 self.patches: dict[str, ParametricSurface] = (

108 {}

109 ) # Parametric patches (continuous)

110 self.patch_res_r = {} # corresponding stl resolution in 'r' direction

111 self.patch_res_s = {} # corresponding stl resolution in 's' direction

112

113 # VTK Attributes

114 self.grids = None # Structured grids

115

116 # STL Attributes

117 self.surfaces = None # STL surfaces for each patch

118 self.stl_resolution = stl_resolution # STL cells per edge

119 self._mesh = None # STL mesh for entire component

120

121 # Curvature functions

122 self._curvatures = None

123

124 # Clustering

125 self._clustering = {}

126

127 # Transformations

128 self._transformations = []

129

130 # Modifier function

131 self._modifier_function = None

132

133 # Component reflection

134 self._reflection_axis = None

135 self._append_reflection = True

136

137 # Ghost component

138 self._ghost = False

139

140 # Component name

141 self.name = name

142

143 # Multiprocessing flag for STL generation

144 # TODO - make this configurable by user

145 self._multiprocess = True

146

147 def __repr__(self):

148 s = f"{self.componenttype} component"

149 if self.name:

150 s += f" (tagged '{self.name}')"

151

152 return s

153

154 def __str__(self):

155 if self.name:

156 return self.name

157 else:

158 return f"{self.componenttype} component"

159

160 def __copy__(self):

161 cls = self.__class__

162 result = cls.__new__(cls)

163 result.__dict__.update(self.__dict__)

164 return result

165

166 def __deepcopy__(self, memo):

167 cls = self.__class__

168 result = cls.__new__(cls)

169 memo[id(self)] = result

170 for k, v in self.__dict__.items():

171 setattr(result, k, deepcopy(v, memo))

172 return result

173

174 def curve(self):

175 if self._curvatures is not None:

176 for curvature in self._curvatures:

177 for key, patch in self.patches.items():

178 self.patches[key] = CurvedPatch(

179 underlying_surf=patch,

180 direction=curvature[0],

181 fun=curvature[1],

182 fun_dash=curvature[2],

183 )

184

185 @property

186 def mesh(self):

187 if not self._mesh:

188 # Check for processed surfaces

189 if self.surfaces is None:

190 if self.verbosity > 1:

191 print(" Generating surfaces for component.")

192

193 # Generate surfaces

194 self.surface()

195

196 # Combine all surface data

197 surface_data = np.concatenate([s[1].data for s in self.surfaces.items()])

198

199 # Create nominal STL mesh

200 self._mesh = mesh.Mesh(surface_data)

201

202 return self._mesh

203

204 @mesh.setter

205 def mesh(self, value):

206 self._mesh = value

207

208 def rotate(self, angle: float = 0, axis: str = "y"):

209 for key, patch in self.patches.items():

210 self.patches[key] = RotatedPatch(patch, np.deg2rad(angle), axis=axis)

211

212 def translate(self, offset: Union[Callable, Vector3]):

213 if isinstance(offset, Vector3):

214 # Translate vector into lambda function

215 offset_function = lambda x, y, z: offset

216 else:

217 offset_function = offset

218

219 # Could wrap it in a lambda if provided

220 for key, patch in self.patches.items():

221 self.patches[key] = OffsetPatchFunction(patch, offset_function)

222

223 def transform(self):

224 for transform in self._transformations:

225 func = getattr(self, transform[0])

226 func(*transform[1:])

227

228 def apply_modifier(self):

229 if self._modifier_function:

230 for key, patch in self.patches.items():

231 self.patches[key] = OffsetPatchFunction(patch, self._modifier_function)

232

233 def reflect(self, axis: str = None):

234 axis = self._reflection_axis if self._reflection_axis is not None else axis

235 if axis is not None:

236 # Create mirrored patches

237 mirrored_patches = {}

238 for key, patch in self.patches.items():

239 mirrored_patches[f"{key}_mirrored"] = MirroredPatch(patch, axis=axis)

240

241 # Now either add the reflected patches to the same component, or

242 # replace the patches with their reflected versions

243 if self._append_reflection:

244 # Append mirrored patches to original patches

245 for key, patch in mirrored_patches.items():

246 self.patches[key] = patch

247 else:

248 # Overwrite existing patches

249 self.patches = mirrored_patches

250

251 def grid(self):

252 for key in self.patches:

253 self.grids[key] = StructuredGrid(

254 surface=self.patches[key],

255 niv=self.vtk_resolution,

256 njv=self.vtk_resolution,

257 )

258

259 def surface(self, resolution: int = None):

260 stl_resolution = self.stl_resolution if resolution is None else resolution

261

262 # Check for patches

263 if len(self.patches) == 0:

264 raise Exception(

265 "No patches have been generated. Please call .generate_patches()."

266 )

267

268 # Create case list

269 if isinstance(stl_resolution, int):

270 case_list = [

271 [k, patch, stl_resolution, stl_resolution]

272 for k, patch in self.patches.items()

273 ]

274 else:

275 case_list = [

276 [k, patch, self.patch_res_r[k], self.patch_res_r[k]]

277 for k, patch in self.patches.items()

278 ]

279

280 # Prepare multiprocessing arguments iterable

281 def wrapper(key: str, patch: ParametricSurface, res_r: int, res_s: int):

282 surface = surfce_to_stl(

283 parametric_surface=patch,

284 triangles_per_edge_r=res_r,

285 triangles_per_edge_s=res_s,

286 **self._clustering,

287 )

288 return (key, surface)

289

290 self.surfaces = {}

291 if self._multiprocess is True:

292 # Initialise surfaces and pool

293 pool = Pool()

294

295 # Submit tasks

296 for result in pool.starmap(wrapper, case_list):

297 self.surfaces[result[0]] = result[1]

298

299 else:

300 for case in case_list:

301 k = case[0]

302 pat = case[1]

303 print(f"START: Creating stl for '{k}'.")

304 result = wrapper(k, pat, case[2], case[3])

305 self.surfaces[result[0]] = result[1]

306 print(f" DONE: Creating stl for '{k}'.")

307

308 def to_vtk(self):

309 raise NotImplementedError("This method has not been implemented yet.")

310 # TODO - check for processed grids

311 for key, grid in self.grids.items():

312 grid.write_to_vtk_file(f"{self.vtk_filename}-wing_{key}.vtk")

313

314 def to_stl(self, outfile: str = None):

315 if not self._ghost:

316 if self.verbosity > 1:

317 print("Writing patches to STL format. ")

318 if outfile is not None:

319 print(f"Output file = {outfile}.")

320

321 # Get mesh

322 stl_mesh = self.mesh

323

324 if outfile is not None:

325 # Write STL to file

326 stl_mesh.save(outfile)

327

328 # Clean it

329 # NOTE - disabling this for now as it can cause severe mesh deformation

330 # print("CLEANING STL FILE WITH PYMESHFIX")

331 # pymeshfix.clean_from_file(outfile, outfile)

332

333 def analyse(self):

334 # Get mass properties

335 volume, cog, inertia = self.mesh.get_mass_properties()

336

337 # Print results

338 print(f"Volume: {volume} m^3")

339 print(f"COG location: {cog}")

340 print("Moment of intertia metrix at COG:")

341 print(inertia)

342

343 def add_clustering_options(

344 self,

345 i_clustering_func: Optional[Callable] = None,

346 j_clustering_func: Optional[Callable] = None,

347 ):

348 """Add a clustering option to this component.

349

350 Parameters

351 -----------

352 i_clustering_func : Callable, optional

353 The clustering function in the i direction. The default is None.

354

355 j_clustering_func : Callable, optional

356 The clustering function in the j direction. The default is None.

357 """

358 if i_clustering_func:

359 self._clustering.update({"i_clustering_func": i_clustering_func})

360

361 if j_clustering_func:

362 self._clustering.update({"j_clustering_func": j_clustering_func})

363

364 def stl_check(

365 self, project_area: Optional[bool] = True, matching_lines: Optional[bool] = True

366 ):

367 """Check the STL mesh."""

368 # TODO - add comment annotations

369 pass_flag = True

370 mesh = self.mesh

371 print(f" MESH CHECK")

372 print(f" mesh: {mesh}")

373 if project_area:

374 print(f" Project Area Check")

375 normals = np.asarray(mesh.normals, dtype=np.float64)

376 allowed_max_errors = np.abs(normals).sum(axis=0) * np.finfo(np.float32).eps

377 sum_normals = normals.sum(axis=0)

378 print(f" allowed error: {allowed_max_errors}")

379 print(f" normals sum: {sum_normals}")

380 if ((np.abs(sum_normals)) <= allowed_max_errors).all():

381 print(f" PASS")

382 else:

383 print(f" FAIL")

384 pass_flag = False

385

386 if matching_lines:

387 print(f" Matching Lines Check")

388 reversed_triangles = (

389 np.cross(mesh.v1 - mesh.v0, mesh.v2 - mesh.v0) * mesh.normals

390 ).sum(axis=1) < 0

391 import itertools

392

393 directed_edges = {

394 tuple(edge.ravel() if not rev else edge[::-1, :].ravel())

395 for rev, edge in zip(

396 itertools.cycle(reversed_triangles),

397 itertools.chain(

398 mesh.vectors[:, (0, 1), :],

399 mesh.vectors[:, (1, 2), :],

400 mesh.vectors[:, (2, 0), :],

401 ),

402 )

403 }

404 undirected_edges = {frozenset((edge[:3], edge[3:])) for edge in directed_edges}

405 edge_check = len(directed_edges) == 3 * mesh.data.size

406 print(f" len(directed_edges) == 3 * mesh.data.size")

407 if edge_check:

408 print(f" {len(directed_edges)} == {3*mesh.data.size}")

409 print(f" PASS")

410 else:

411 print(f" {len(directed_edges)} != {3*mesh.data.size}")

412 print(f" FAIL")

413 print(f" The number of edges should be N_cells*3")

414 print(f" len(directed_edges)={len(directed_edges)}")

415 print(f" mesh.data.size={mesh.data.size}")

416 pass_flag = False

417 pair_check = len(directed_edges) == 2 * len(undirected_edges)

418 print(f" len(directed_edges) == 2 * len(undirected_edges)")

419 if pair_check:

420 print(f" {len(directed_edges)} == {2 * len(undirected_edges)}")

421 print(f" PASS")

422 return pass_flag

423 else:

424 print(f" {len(directed_edges)} != {2 * len(undirected_edges)}")

425 print(f" FAIL")

426 print(f" All edges should be in pairs")

427 print(f" len_directed:{len(directed_edges)}")

428 print(f" len_undirect (pairs+leftover):{len(undirected_edges)}")

429 pass_flag = False

430

431 edge_list = []

432 for edge in directed_edges:

433 edge_list.append(frozenset((edge[:3], edge[3:])))

434

435 counter_out = Counter(edge_list)

436

437 k_list = []

438 for k, v in counter_out.items():

439 if v == 2:

440 k_list.append(k)

441 for k in k_list:

442 del counter_out[k]

443

444 return pass_flag

445

446 def stl_repair(self, small_distance: float = 1e-6):

447 """Attempts to repair stl mesh issues.

448

449 Parameters

450 -----------

451 small_distance : Float, optional

452 Vectors less than this distance apart will be set to the same value.

453 """

454 mesh = self.mesh

455 N_faces = np.shape(mesh)[0]

456 vectors = np.vstack((mesh.v0, mesh.v1, mesh.v2))

457

458 def find_groups(vector, small_distance):

459 """find indices for groups of points that within small_distance from one another."""

460 sort_i = np.argsort(vector)

461 groups = []

462 li = []

463 count = 0

464 for i in range(len(vector) - 1):

465 if count == 0:

466 x0 = vector[sort_i[i]]

467 x1 = vector[sort_i[i + 1]]

468 if abs(x1 - x0) < small_distance:

469 if count == 0:

470 li.append(sort_i[i])

471 li.append(sort_i[i + 1])

472 count = 1

473 else:

474 li.append(sort_i[i + 1])

475 else:

476 groups.append(li)

477 li = []

478 count = 0

479 if i == len(vector) - 2 and count > 0:

480 groups.append(li)

481 return groups

482

483 iix = find_groups(vectors[:, 0], small_distance)

484 iiy = find_groups(vectors[:, 1], small_distance)

485 iiz = find_groups(vectors[:, 2], small_distance)

486

487 # find intersecting sets and take average

488 for ix in iix:

489 for iy in iiy:

490 common0 = set(ix) & set(iy)

491 if len(common0) == 0:

492 continue # jumpt to next loop iteration if only one entry

493 for iz in iiz:

494 common = list(common0 & set(iz))

495 if common:

496 vectors[common] = np.mean(vectors[common], axis=0)

497 mesh.v0 = vectors[0:N_faces, :]

498 mesh.v1 = vectors[N_faces : 2 * N_faces, :]

499 mesh.v2 = vectors[2 * N_faces : 3 * N_faces, :]