Coverage for /opt/hostedtoolcache/Python/3.11.9/x64/lib/python3.11/site-packages/hypervehicle/geometry/geometry.py: 52%

514 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2024-08-25 22:58 +0000

1import numpy as np 

2from hypervehicle.geometry.vector import Vector3 

3from hypervehicle.geometry.surface import CoonsPatch, ParametricSurface 

4from hypervehicle.geometry.path import Line, Path, ArcLengthParameterizedPath 

5 

6 

7class SubRangedPath(Path): 

8 """ 

9 A Path reparameterized to a subset of the original between t0 and t1. 

10 If t1 < t0 the path will be traveresed in the reverse direction. 

11 """ 

12 

13 __slots__ = ["underlying_path", "t0", "t1"] 

14 

15 def __init__(self, underlying_path, t0, t1): 

16 if isinstance(underlying_path, Path): 

17 self.underlying_path = underlying_path 

18 self.t0 = t0 

19 self.t1 = t1 

20 else: 

21 raise NotImplementedError("underlying_path should be a type of Path") 

22 

23 def __repr__(self) -> str: 

24 return "ArcLengthParameterizedPath(underlying_path={}, t0={}, t1={})".format( 

25 self.underlying_path, self.t0, self.t1 

26 ) 

27 

28 def __str__(self) -> str: 

29 return "ArcLengthParameterizedPath" 

30 

31 def __call__(self, t): 

32 t_dash = self.t0 + t * (self.t1 - self.t0) 

33 return self.underlying_path(t_dash) 

34 

35 def length(self): 

36 return self.underlying_path.length() 

37 

38 

39class ReversedPath(SubRangedPath): 

40 def __init__(self, underlying_path: Path) -> None: 

41 self.underlying_path = underlying_path 

42 self.t0 = 1 

43 self.t1 = 0 

44 

45 

46class ElipsePath(Path): 

47 """ 

48 A path following a quarter elipse from a -> b, around c 

49 

50 b - _ 

51 -_ 

52 \ 

53 c a 

54 

55 """ 

56 

57 __slots__ = ["centre", "thickness", "LE_width", "side"] 

58 

59 def __init__(self, centre, thickness, LE_width, side): 

60 self.centre = centre 

61 self.thickness = thickness 

62 self.LE_width = LE_width 

63 self.side = side 

64 

65 def __repr__(self): 

66 return "ElipsePath" 

67 

68 def __call__(self, r): 

69 # establish local elipse shapes 

70 a = self.LE_width 

71 b = abs(self.thickness) 

72 # a = b*self.LE_ratio(t).y 

73 

74 # construct elipse 

75 if self.side == "bot": 

76 theta = 0 + r * np.pi / 2 

77 elif self.side == "top": 

78 theta = -np.pi / 2 + r * np.pi / 2 

79 else: 

80 raise Exception("Value of 'side' not supported") 

81 

82 x_elipse = a * b / np.sqrt(b * b + (a * np.tan(theta)) ** 2) 

83 y_elipse = x_elipse * np.tan(theta) 

84 # set elipse angle 

85 angle = np.pi / 2 

86 # add elipse to overall shape 

87 x = x_elipse * np.cos(angle) 

88 y = x_elipse * np.sin(angle) 

89 z = y_elipse 

90 elipse_base = self.centre 

91 return elipse_base + Vector3(x=x, y=y, z=z) 

92 

93 

94class OffsetPathFunction(Path): 

95 """ 

96 Creates a line with an offset applied 

97 """ 

98 

99 __slots__ = ["underlying_line", "offset_function"] 

100 

101 def __init__(self, underlying_line, offset_function): 

102 self.underlying_line = underlying_line 

103 self.offset_function = offset_function 

104 

105 def __repr__(self): 

106 return "Offset path function" 

107 

108 def __call__(self, t): 

109 # calculate postion 

110 pos = self.underlying_line(t) 

111 

112 # calculate local thickness 

113 offset = self.offset_function(x=pos.x, y=pos.y) 

114 

115 return pos + offset 

116 

117 

118class GeometricMeanPathFunction(Path): 

119 """ 

120 Creates a line which is the geometric mean of two underlying lines. 

121 """ 

122 

123 __slots__ = ["underlying_line_1", "underlying_line_2"] 

124 

125 def __init__(self, underlying_line_1, underlying_line_2): 

126 self.underlying_line_1 = underlying_line_1 

127 self.underlying_line_2 = underlying_line_2 

128 

129 def __repr__(self): 

130 return "Mean path function" 

131 

132 def __call__(self, t): 

133 # calculate postion 

134 pos_1 = self.underlying_line_1(t) 

135 pos_2 = self.underlying_line_2(t) 

136 

137 # calculate local thickness 

138 mean_point = 0.5 * (pos_1 + pos_2) 

139 

140 return mean_point 

141 

142 

143class OffsetPatchFunction(ParametricSurface): 

144 """ 

145 Creates a patch with an offset applied. 

146 """ 

147 

148 __slots__ = ["underlying_surf", "function"] 

149 

150 def __init__(self, underlying_surf, function): 

151 self.underlying_surf = underlying_surf 

152 self.function = function 

153 

154 def __repr__(self): 

155 str = " + offset function" 

156 str = self.underlying_surf.__repr__() + str 

157 return str 

158 

159 def __call__(self, r, s): 

160 pos = self.underlying_surf(r, s) 

161 offset = self.function(pos.x, pos.y, pos.z) 

162 return self.underlying_surf(r, s) + offset 

163 

164 

165class LeadingEdgePatchFunction(ParametricSurface): 

166 """Creates Leading Edge by pair of guiding lines and LE_width function.""" 

167 

168 __slots__ = [ 

169 "centralLine", 

170 "thickness_function", 

171 "LE_width_function", 

172 "t0", 

173 "t1", 

174 "side", 

175 ] 

176 

177 def __init__( 

178 self, centralLine, thickness_function, LE_width_function, t0, t1, side="top" 

179 ): 

180 self.centralLine = centralLine 

181 self.thickness_function = thickness_function 

182 self.t0 = t0 

183 self.t1 = t1 

184 self.LE_width_function = LE_width_function 

185 self.side = side 

186 

187 def __repr__(self): 

188 return "Leading edge surface" 

189 

190 def __call__(self, t_raw, r): 

191 # convert to global t 

192 # t = self.t0 + t_raw * (self.t1 - self.t0) 

193 t = self.t1 - t_raw * (self.t1 - self.t0) 

194 # calculate local thickness 

195 pos = self.centralLine(t) 

196 thickness = self.thickness_function( 

197 x=pos.x, y=pos.y 

198 ).z # we only need to get z-value 

199 LE_width = self.LE_width_function(t) 

200 

201 # establish local elipse shapes 

202 elipse = ElipsePath( 

203 centre=Vector3(x=0.0, y=0.0, z=0), 

204 thickness=thickness, 

205 LE_width=LE_width, 

206 side=self.side, 

207 ) 

208 elipse_norm = ArcLengthParameterizedPath(underlying_path=elipse) 

209 pos_elipse = elipse_norm(r) 

210 x_elipse = pos_elipse.y 

211 # y_elipse = pos_elipse.x 

212 z_elipse = pos_elipse.z 

213 

214 # set elipse angle 

215 if t == 0: 

216 angle = np.pi / 2 

217 elif t == 1: 

218 angle = 0 

219 else: 

220 dt = min([0.001, abs(1.0 - t), abs(t - 0.0)]) 

221 plus = self.centralLine(t + dt) 

222 minus = self.centralLine(t - dt) 

223 angle = np.arctan2((plus.y - minus.y), (plus.x - minus.x)) + np.pi / 2 

224 

225 # print(angle, x_elipse) 

226 

227 # add elipse to overall shape 

228 x = x_elipse * np.cos(angle) 

229 y = x_elipse * np.sin(angle) 

230 z = z_elipse 

231 elipse_base = self.centralLine(t) 

232 return elipse_base + Vector3(x=x, y=y, z=z) 

233 

234 

235class MeanLeadingEdgePatchFunction(ParametricSurface): 

236 """Creates Leading Edge by mean line and guiding line, and LE_width function.""" 

237 

238 __slots__ = ["mean_line", "guide_line", "LE_width_function", "t0", "t1", "side"] 

239 

240 def __init__(self, mean_line, guide_line, LE_width_function, t0, t1, side="top"): 

241 self.mean_line = mean_line 

242 self.guide_line = guide_line 

243 self.t0 = t0 

244 self.t1 = t1 

245 self.LE_width_function = LE_width_function 

246 self.side = side 

247 

248 def __repr__(self): 

249 return "Leading edge surface patch" 

250 

251 def __call__(self, t_raw, r): 

252 # Convert to global t 

253 t = self.t1 - t_raw * (self.t1 - self.t0) 

254 

255 # Calculate local thickness 

256 mean_point = self.mean_line(t) 

257 guide_point = self.guide_line(t) 

258 thickness = mean_point.z - guide_point.z 

259 LE_width = self.LE_width_function(t) 

260 

261 # Establish local elipse shapes 

262 elipse = ElipsePath( 

263 centre=Vector3(x=0.0, y=0.0, z=0), 

264 thickness=thickness, 

265 LE_width=LE_width, 

266 side=self.side, 

267 ) 

268 elipse_norm = ArcLengthParameterizedPath(underlying_path=elipse) 

269 pos_elipse = elipse_norm(r) 

270 x_elipse = pos_elipse.y 

271 z_elipse = pos_elipse.z 

272 

273 # set elipse angle 

274 if t == 0: 

275 angle = np.pi / 2 

276 elif t == 1: 

277 angle = 0 

278 else: 

279 dt = min([0.001, abs(1.0 - t), abs(t - 0.0)]) 

280 plus = self.mean_line(t + dt) 

281 minus = self.mean_line(t - dt) 

282 angle = np.arctan2((plus.y - minus.y), (plus.x - minus.x)) + np.pi / 2 

283 

284 # Add elipse to overall shape 

285 x = x_elipse * np.cos(angle) 

286 y = x_elipse * np.sin(angle) 

287 z = z_elipse 

288 

289 return mean_point + Vector3(x=x, y=y, z=z) 

290 

291 

292class FlatLeadingEdgePatchFunction(ParametricSurface): 

293 """ 

294 Creates a flat leading edge between two paths. 

295 

296 TODO - adapt this to allow sharp LE 

297 """ 

298 

299 __slots__ = ["path1", "path2", "t0", "t1"] 

300 

301 def __init__(self, path1, path2, t0, t1): 

302 self.path1 = path1 

303 self.path2 = path2 

304 

305 # Parametric section of mean and guide line 

306 self.t0 = t0 

307 self.t1 = t1 

308 

309 def __repr__(self): 

310 return "Flat leading edge surface patch" 

311 

312 def __call__(self, t_raw, r): 

313 # Convert to global t 

314 t = self.t1 - t_raw * (self.t1 - self.t0) 

315 

316 # Calculate local thickness 

317 point1 = self.path1(t) 

318 point2 = self.path2(t) 

319 

320 return Vector3(x=point1.x, y=point1.y, z=(1 - r) * point1.z + r * point2.z) 

321 

322 

323class TrailingEdgePath(Path): 

324 """ 

325 A path following rear of wing 

326 """ 

327 

328 __slots__ = ["A0", "B0", "thickness_function"] 

329 

330 def __init__(self, A0, B0, thickness_function): 

331 self.A0 = A0 

332 self.B0 = B0 

333 self.thickness_function = thickness_function 

334 self.Line = Line(p0=A0, p1=B0) 

335 

336 def __repr__(self): 

337 return "Trailing Egde Path" 

338 

339 def __call__(self, r): 

340 # calculate postion 

341 pos = self.Line(r) 

342 # calculate local thickness 

343 offset = self.thickness_function(x=pos.x, y=pos.y) 

344 return pos + offset 

345 

346 

347class TrailingEdgePatch(ParametricSurface): 

348 """ 

349 Fucntion to create a trailing edge patch 

350 """ 

351 

352 __slots__ = ["A0", "B0", "TE_path", "flap_length", "flap_angle", "side"] 

353 

354 def __init__(self, A0, B0, TE_path, flap_length, flap_angle=0.0, side="top"): 

355 self.A0 = A0 

356 self.B0 = B0 

357 self.TE_path = TE_path 

358 self.flap_length = flap_length 

359 self.flap_angle = flap_angle 

360 self.side = side 

361 self.preparation() 

362 

363 def __repr__(self): 

364 return "Trailing Edge for '{0}' with flap_angle={1} deg".format( 

365 self.side, np.rad2deg(self.flap_angle) 

366 ) 

367 

368 def __call__(self, r, s): 

369 return self.patch(r, s) 

370 

371 def preparation(self): 

372 if self.side == "top": 

373 north = self.TE_path 

374 south = Line( 

375 p0=Vector3( 

376 x=self.A0.x - self.flap_length, 

377 y=self.A0.y, 

378 z=0.0 + self.flap_length * np.sin(self.flap_angle), 

379 ), 

380 p1=Vector3( 

381 x=self.B0.x - self.flap_length, 

382 y=self.B0.y, 

383 z=0.0 + self.flap_length * np.sin(self.flap_angle), 

384 ), 

385 ) 

386 west = Line(p0=south(0.0), p1=north(0.0)) 

387 east = Line(p0=south(1.0), p1=north(1.0)) 

388 elif self.side == "bot": 

389 south = self.TE_path 

390 north = Line( 

391 p0=Vector3( 

392 x=self.A0.x - self.flap_length, 

393 y=self.A0.y, 

394 z=0.0 + self.flap_length * np.sin(self.flap_angle), 

395 ), 

396 p1=Vector3( 

397 x=self.B0.x - self.flap_length, 

398 y=self.B0.y, 

399 z=0.0 + self.flap_length * np.sin(self.flap_angle), 

400 ), 

401 ) 

402 west = Line(p0=south(0.0), p1=north(0.0)) 

403 east = Line(p0=south(1.0), p1=north(1.0)) 

404 else: 

405 raise Exception("Value of 'side' not supported") 

406 self.patch = CoonsPatch(south=south, north=north, west=west, east=east) 

407 

408 

409class MeanTrailingEdgePatch(ParametricSurface): 

410 """ 

411 Fucntion to create a trailing edge patch, with flap angle defined 

412 from the geometric mean of upper and lower surfaces. 

413 """ 

414 

415 __slots__ = ["mean_line", "TE_path", "flap_length", "flap_angle", "side"] 

416 

417 def __init__(self, mean_line, TE_path, flap_length, flap_angle=0, side="top"): 

418 self.mean_line = mean_line 

419 self.TE_path = TE_path 

420 self.flap_length = flap_length 

421 self.flap_angle = flap_angle 

422 self.side = side 

423 self.preparation() 

424 

425 def __repr__(self): 

426 return f"Trailing Edge for '{self.side}' with flap_angle={np.rad2deg(self.flap_angle)} deg" 

427 

428 def __call__(self, r, s): 

429 return self.patch(r, s) 

430 

431 def preparation(self): 

432 if self.side == "top": 

433 north = self.TE_path 

434 south = Line( 

435 p0=Vector3( 

436 x=self.mean_line(0).x - self.flap_length * np.cos(self.flap_angle), 

437 y=self.mean_line(0).y, 

438 z=self.mean_line(0).z + self.flap_length * np.sin(self.flap_angle), 

439 ), 

440 p1=Vector3( 

441 x=self.mean_line(1).x - self.flap_length * np.cos(self.flap_angle), 

442 y=self.mean_line(1).y, 

443 z=self.mean_line(1).z + self.flap_length * np.sin(self.flap_angle), 

444 ), 

445 ) 

446 

447 west = Line(p0=south(0.0), p1=north(0.0)) 

448 east = Line(p0=south(1.0), p1=north(1.0)) 

449 

450 elif self.side == "bot": 

451 south = self.TE_path 

452 north = Line( 

453 p0=Vector3( 

454 x=self.mean_line(0).x - self.flap_length * np.cos(self.flap_angle), 

455 y=self.mean_line(0).y, 

456 z=self.mean_line(0).z + self.flap_length * np.sin(self.flap_angle), 

457 ), 

458 p1=Vector3( 

459 x=self.mean_line(1).x - self.flap_length * np.cos(self.flap_angle), 

460 y=self.mean_line(1).y, 

461 z=self.mean_line(1).z + self.flap_length * np.sin(self.flap_angle), 

462 ), 

463 ) 

464 west = Line(p0=south(0.0), p1=north(0.0)) 

465 east = Line(p0=south(1.0), p1=north(1.0)) 

466 

467 else: 

468 raise Exception("Value of 'side' not supported") 

469 

470 self.patch = CoonsPatch(south=south, north=north, west=west, east=east) 

471 

472 

473class CurvedPatch(ParametricSurface): 

474 """Adds curvature in x or y direction to an existing patch""" 

475 

476 __slots__ = ["underlying_surf", "direction", "fun", "fun_dash"] 

477 

478 def __init__(self, underlying_surf, direction=None, fun=None, fun_dash=None): 

479 self.underlying_surf = underlying_surf 

480 self.direction = direction 

481 self.fun = fun 

482 self.fun_dash = fun_dash 

483 

484 def __repr__(self): 

485 return self.underlying_surf.__repr__() + " with added curvature" 

486 

487 def __call__(self, r, s): 

488 if self.fun == None or self.fun_dash == None: 

489 raise Exception("Both 'fun' and 'fun_dash' need to be specified.") 

490 pos = self.underlying_surf(r, s) 

491 offset = self.fun(pos.x, pos.y) 

492 slope = self.fun_dash(pos.x, pos.y) 

493 angle = np.arctan2(slope, 1) 

494 if self.direction == "x": 

495 pos_new = Vector3( 

496 x=pos.x - pos.z * np.sin(angle), 

497 y=pos.y, 

498 z=pos.z * np.cos(angle) + offset, 

499 ) 

500 return pos_new 

501 if self.direction == "y": 

502 pos_new = Vector3( 

503 x=pos.x, 

504 y=pos.y - pos.z * np.sin(angle), 

505 z=offset + pos.z * np.cos(angle), 

506 ) 

507 return pos_new 

508 

509 

510class ConePatch(ParametricSurface): 

511 """ 

512 Creates a patch describing a cone (or cylinder) between two rings. 

513 """ 

514 

515 __slots__ = ["x0", "x1", "r0", "r1", "angle0", "angle1"] 

516 

517 def __init__(self, x0, x1, r0, r1, angle0, angle1): 

518 self.x0 = x0 

519 self.x1 = x1 

520 self.r0 = r0 

521 self.r1 = r1 

522 self.angle0 = angle0 

523 self.angle1 = angle1 

524 

525 def __repr__(self): 

526 str = "Cone Patch" 

527 return str 

528 

529 def __call__(self, r, s): 

530 angle = self.angle0 * (1 - s) + self.angle1 * s 

531 x0 = self.x0 

532 y0 = self.r0 * np.cos(angle) 

533 z0 = self.r0 * np.sin(angle) 

534 x1 = self.x1 

535 y1 = self.r1 * np.cos(angle) 

536 z1 = self.r1 * np.sin(angle) 

537 x = x0 * (1 - r) + x1 * r 

538 y = y0 * (1 - r) + y1 * r 

539 z = z0 * (1 - r) + z1 * r 

540 return Vector3(x=x, y=y, z=z) 

541 

542 

543class RevolvedPatch(ParametricSurface): 

544 """Creates a path by revolving a line about a central 

545 axis. 

546 """ 

547 

548 def __init__(self, line, angle0=0, angle1=2 * np.pi): 

549 self.line = line 

550 self.angle0 = angle0 

551 self.angle1 = angle1 

552 

553 def __repr__(self): 

554 return "Revolved Patch" 

555 

556 def __call__(self, r, s): 

557 # Map s to angular coordinate 

558 angle = self.angle0 * (1 - s) + self.angle1 * s 

559 

560 # Map r to distance along line 

561 point = self.line(r) 

562 

563 # Calculate points 

564 x = point.x 

565 y = point.y * np.cos(angle) + point.z * np.sin(angle) 

566 z = point.z * np.cos(angle) - point.y * np.sin(angle) 

567 

568 return Vector3(x=x, y=y, z=z) 

569 

570 

571class BluntConePatch(ParametricSurface): 

572 """ 

573 Creates a patch describing a blunt cone. 

574 

575 Parameters: 

576 x0: x coordinate of base of cone 

577 y0: y coordinate of base of cone 

578 rn: spherical nose radius 

579 rb: radius of cone at base of cone 

580 L: length of nominal cone (base to tip before spherical blunting) 

581 angle0: start revolve angle 

582 angle1: end revolve angle 

583 

584 """ 

585 

586 __slots__ = ["x0", "y0", "rn", "rb", "L", "angle0", "angle1"] 

587 

588 def __init__(self, x0, y0, rn, rb, L, angle0, angle1): 

589 self.L = L 

590 self.x0 = x0 

591 self.y0 = y0 

592 self.rn = rn 

593 self.rb = rb 

594 self.angle0 = angle0 

595 self.angle1 = angle1 

596 

597 def __repr__(self): 

598 str = "Blunt Cone Patch" 

599 return str 

600 

601 def __call__(self, r, s): 

602 angle = self.angle0 * (1 - s) + self.angle1 * s 

603 x0 = self.x0 

604 y0 = self.r0 * np.cos(angle) 

605 z0 = self.r0 * np.sin(angle) 

606 x1 = self.x1 

607 y1 = self.r1 * np.cos(angle) 

608 z1 = self.r1 * np.sin(angle) 

609 x = x0 * (1 - r) + x1 * r 

610 y = y0 * (1 - r) + y1 * r 

611 z = z0 * (1 - r) + z1 * r 

612 return Vector3(x=x, y=y, z=z) 

613 

614 

615class SweptPatch(ParametricSurface): 

616 """Creates a swept patch from a series of cross sections.""" 

617 

618 __slots__ = ["cross_sections", "section_origins"] 

619 

620 def __init__(self, cross_sections: list, sweep_axis: str = "z") -> None: 

621 """Construct the SweptPatch object. 

622 

623 Parameters 

624 ----------- 

625 cross_sections : list 

626 A list containing the cross sections to be blended. These must 

627 be defined in the x-y plane. 

628 sweep_axis : str, optional 

629 The axis to sweep the cross sections along. Note that each cross 

630 section should vary along this axis. The default is "z". 

631 """ 

632 self.cross_sections = cross_sections 

633 self.section_origins = [getattr(cs(0, 0), sweep_axis) for cs in cross_sections] 

634 self._sweep_axis = sweep_axis 

635 self._other_axes = {"x", "y", "z"} - set(sweep_axis) 

636 

637 if min(self.section_origins) == max(self.section_origins): 

638 raise Exception( 

639 "There is no axial variation in the cross " + "sections provided!" 

640 ) 

641 

642 self.perimeters = [SurfacePerimeter(s) for s in cross_sections] 

643 self.min_origin = min(self.section_origins) 

644 self.origin_dist = max(self.section_origins) - self.min_origin 

645 

646 def __repr__(self): 

647 return "Swept Patch" 

648 

649 def __call__(self, r, s) -> Vector3: 

650 # Calculate physical axial distance 

651 dist = self.min_origin + r * self.origin_dist 

652 

653 # Find index of bounding cross sections 

654 for i, lv in enumerate(self.section_origins): 

655 if dist == self.section_origins[-1]: 

656 i = len(self.section_origins) - 2 

657 break 

658 elif lv <= dist < self.section_origins[i + 1]: 

659 break 

660 

661 # Get upper and lower perimeter 

662 lps = self.perimeters[i](s) 

663 ups = self.perimeters[i + 1](s) 

664 

665 # Get upper and lower bounding sections 

666 ls = self.section_origins[i] 

667 us = self.section_origins[i + 1] 

668 

669 # Linearly interpolate between cross-sectional perimeters 

670 args = { 

671 a: getattr(lps, a) 

672 + (dist - ls) * (getattr(ups, a) - getattr(lps, a)) / (us - ls) 

673 for a in self._other_axes 

674 } 

675 args[self._sweep_axis] = dist 

676 

677 return Vector3(**args) 

678 

679 

680class SweptPatchfromEdges(ParametricSurface): 

681 """Creates a swept patch from a series of cross sections. 

682 Cross sections do not need to be parallel or aligned.""" 

683 

684 __slots__ = ["edges"] 

685 

686 def __init__(self, edges: list) -> None: 

687 """Construct the SweptPatch object. 

688 

689 Parameters 

690 ----------- 

691 edges : list 

692 A list containing the edges through which to sweep. 

693 """ 

694 self.edges = edges 

695 self.n_edges = len(edges) 

696 

697 def __repr__(self): 

698 return "Swept Patch" 

699 

700 def __call__(self, r, s) -> Vector3: 

701 r_slices = np.linspace(0, 1, self.n_edges) 

702 if r in r_slices: 

703 # can evaluate exactly 

704 i = np.where(r_slices == r)[0][0] 

705 return self.edges[i](s) 

706 else: 

707 # need to interpolate 

708 index_above = np.where(r_slices > r)[0][0] 

709 index_below = index_above - 1 

710 above = self.edges[index_above](s) 

711 below = self.edges[index_below](s) 

712 alpha = (r - r_slices[index_below]) / ( 

713 r_slices[index_above] - r_slices[index_below] 

714 ) 

715 return (1 - alpha) * below + alpha * above 

716 

717 

718class RotatedPatch(ParametricSurface): 

719 """ 

720 Rotates a surface about a point in an axis-specified direction. 

721 """ 

722 

723 __slots__ = ["underlying_surf", "angle", "axis", "point"] 

724 

725 def __init__(self, underlying_surf, angle, axis="x", point=Vector3(x=0, y=0, z=0)): 

726 self.underlying_surf = underlying_surf 

727 self.angle = angle 

728 self.axis = axis.lower() 

729 self.point = point 

730 

731 def __repr__(self): 

732 str = f" (rotated by {np.rad2deg(self.angle)} degrees)" 

733 return self.underlying_surf.__repr__() + str 

734 

735 def __call__(self, r, s): 

736 pos = self.underlying_surf(r, s) - self.point 

737 

738 if self.axis == "x": 

739 x = pos.x 

740 y = pos.y * np.cos(self.angle) - pos.z * np.sin(self.angle) 

741 z = pos.y * np.sin(self.angle) + pos.z * np.cos(self.angle) 

742 elif self.axis == "y": 

743 x = pos.x * np.cos(self.angle) + pos.z * np.sin(self.angle) 

744 y = pos.y 

745 z = -pos.x * np.sin(self.angle) + pos.z * np.cos(self.angle) 

746 

747 elif self.axis == "z": 

748 x = pos.x * np.cos(self.angle) - pos.y * np.sin(self.angle) 

749 y = pos.x * np.sin(self.angle) + pos.y * np.cos(self.angle) 

750 z = pos.z 

751 

752 return Vector3(x=x, y=y, z=z) + self.point 

753 

754 

755class MirroredPatch(ParametricSurface): 

756 """Mirrors a surface in an axis-specified direction.""" 

757 

758 __slots__ = ["underlying_surf", "axis"] 

759 

760 def __init__(self, underlying_surf, axis="x"): 

761 self.underlying_surf = underlying_surf 

762 self.axis = axis.lower() 

763 

764 def __repr__(self): 

765 return self.underlying_surf.__repr__() + f" mirrored along {self.axis}-axis" 

766 

767 def __call__(self, r, s): 

768 pos = self.underlying_surf(r, s) 

769 x = pos.x 

770 y = pos.y 

771 z = pos.z 

772 if self.axis == "x": 

773 # Mirror about y-z plane 

774 x *= -1 

775 elif self.axis == "y": 

776 # Mirror about x-z plane 

777 y *= -1 

778 elif self.axis == "z": 

779 # Mirror about x-y plane 

780 z *= -1 

781 

782 return Vector3(x=x, y=y, z=z) 

783 

784 

785class CubePatch(ParametricSurface): 

786 """Creates a cube face patch for a cube of length 

787 2a about the Vector3 centre. 

788 """ 

789 

790 __slots__ = ["a", "centre", "face"] 

791 

792 def __init__(self, a, centre, face): 

793 self.a = a 

794 self.face = face 

795 self.centre = centre 

796 

797 def __repr__(self): 

798 return "Cube: a = {}, centre = {}".format( 

799 self.r, self.centre 

800 ) + "face = {}".format(self.face) 

801 

802 def __call__(self, r, s): 

803 if self.face == "east": 

804 x = 1.0 

805 y = -1.0 + 2.0 * r 

806 z = -1.0 + 2.0 * s 

807 elif self.face == "west": 

808 x = -1.0 

809 y = -1.0 + 2.0 * r 

810 z = -1.0 + 2.0 * s 

811 elif self.face == "south": 

812 x = -1.0 + 2.0 * r 

813 y = -1.0 

814 z = -1.0 + 2.0 * s 

815 elif self.face == "north": 

816 x = -1.0 + 2.0 * r 

817 y = 1.0 

818 z = -1.0 + 2.0 * s 

819 elif self.face == "bottom": 

820 x = -1.0 + 2.0 * r 

821 y = -1.0 + 2.0 * s 

822 z = -1.0 

823 elif self.face == "top": 

824 x = -1.0 + 2.0 * r 

825 y = -1.0 + 2.0 * s 

826 z = 1.0 

827 

828 else: 

829 raise ValueError( 

830 "Incorrect face name." 

831 + "Allowable faces are: east, west, south, north, bottom or top." 

832 ) 

833 

834 x_cube = x * self.a + self.centre.x 

835 y_cube = y * self.a + self.centre.y 

836 z_cube = z * self.a + self.centre.z 

837 

838 return Vector3(x=x_cube, y=y_cube, z=z_cube) 

839 

840 

841class SpherePatch(ParametricSurface): 

842 """Creates a sphere face patch for a cube of length 

843 2a about the Vector3 centre. 

844 """ 

845 

846 __slots__ = ["r", "centre", "face"] 

847 

848 def __init__(self, r, centre, face): 

849 self.r = r 

850 self.face = face 

851 self.centre = centre 

852 

853 def __repr__(self): 

854 return "Sphere: r = {}, centre = {}".format( 

855 self.r, self.centre 

856 ) + "face = {}".format(self.face) 

857 

858 def __call__(self, r, s): 

859 # First create a cube face, then map to a sphere. 

860 # Can change below to call a cube patch rather than 

861 # computing cube face here, it just adds a dependency. 

862 

863 if self.face == "east": 

864 x = 1.0 

865 y = -1.0 + 2.0 * r 

866 z = -1.0 + 2.0 * s 

867 elif self.face == "west": 

868 x = -1.0 

869 y = -1.0 + 2.0 * r 

870 z = -1.0 + 2.0 * s 

871 elif self.face == "south": 

872 x = -1.0 + 2.0 * r 

873 y = -1.0 

874 z = -1.0 + 2.0 * s 

875 elif self.face == "north": 

876 x = -1.0 + 2.0 * r 

877 y = 1.0 

878 z = -1.0 + 2.0 * s 

879 elif self.face == "bottom": 

880 x = -1.0 + 2.0 * r 

881 y = -1.0 + 2.0 * s 

882 z = -1.0 

883 elif self.face == "top": 

884 x = -1.0 + 2.0 * r 

885 y = -1.0 + 2.0 * s 

886 z = 1.0 

887 

888 else: 

889 raise ValueError( 

890 "Incorrect face name." 

891 + "Allowable faces are: east, west, south, north, bottom or top." 

892 ) 

893 

894 x_dash = x * np.sqrt(1.0 - 0.5 * z * z - 0.5 * y * y + y * y * z * z / 3.0) 

895 y_dash = y * np.sqrt(1.0 - 0.5 * z * z - 0.5 * x * x + x * x * z * z / 3.0) 

896 z_dash = z * np.sqrt(1.0 - 0.5 * y * y - 0.5 * x * x + x * x * y * y / 3.0) 

897 

898 x_sphere = x_dash * self.r + self.centre.x 

899 y_sphere = y_dash * self.r + self.centre.y 

900 z_sphere = z_dash * self.r + self.centre.z 

901 

902 return Vector3(x=x_sphere, y=y_sphere, z=z_sphere) 

903 

904 

905class SurfacePerimeter(Path): 

906 """Returns a path corresponding to the perimiter of an underlying 

907 surface.""" 

908 

909 __slots__ = ["underlying_surf"] 

910 

911 def __init__(self, underlying_surf: ParametricSurface) -> None: 

912 self.underlying_surf = underlying_surf 

913 

914 def __repr__(self): 

915 return "Surface Perimeter Path" 

916 

917 def __call__(self, t) -> Vector3: 

918 """Returns the point on the perimeter.""" 

919 f = 0.25 

920 rem = t % f 

921 

922 # Calculate r,s along perimeter 

923 if t < 0.25: 

924 r = round(rem / f, 8) 

925 s = 0.0 

926 elif t < 0.5: 

927 r = 1.0 

928 s = round(rem / f, 8) 

929 elif t < 0.75: 

930 r = round(1 - rem / f, 8) 

931 s = 1.0 

932 elif t < 1: 

933 r = 0.0 

934 s = round(1 - rem / f, 8) 

935 else: 

936 r = 0.0 

937 s = 0.0 

938 

939 return self.underlying_surf(r, s)