Coverage for src/meshpy/core/mesh.py: 95%

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 defines the Mesh class, which holds the content (nodes, 

23elements, sets, ...) for a meshed geometry.""" 

24 

25import copy as _copy 

26import os as _os 

27import warnings as _warnings 

28from typing import Dict as _Dict 

29from typing import List as _List 

30from typing import Optional as _Optional 

31 

32import numpy as _np 

33import pyvista as _pv 

34 

35from meshpy.core.boundary_condition import ( 

36 BoundaryConditionBase as _BoundaryConditionBase, 

37) 

38from meshpy.core.boundary_condition import ( 

39 BoundaryConditionContainer as _BoundaryConditionContainer, 

40) 

41from meshpy.core.conf import mpy as _mpy 

42from meshpy.core.coupling import coupling_factory as _coupling_factory 

43from meshpy.core.element import Element as _Element 

44from meshpy.core.element_beam import Beam as _Beam 

45from meshpy.core.geometry_set import GeometryName as _GeometryName 

46from meshpy.core.geometry_set import GeometrySetBase as _GeometrySetBase 

47from meshpy.core.geometry_set import GeometrySetContainer as _GeometrySetContainer 

48from meshpy.core.material import Material as _Material 

49from meshpy.core.node import Node as _Node 

50from meshpy.core.node import NodeCosserat as _NodeCosserat 

51from meshpy.core.rotation import Rotation as _Rotation 

52from meshpy.core.rotation import add_rotations as _add_rotations 

53from meshpy.core.rotation import rotate_coordinates as _rotate_coordinates 

54from meshpy.core.vtk_writer import VTKWriter as _VTKWriter 

55from meshpy.four_c.function import Function as _Function 

56from meshpy.geometric_search.find_close_points import ( 

57 find_close_points as _find_close_points, 

58) 

59from meshpy.geometric_search.find_close_points import ( 

60 point_partners_to_partner_indices as _point_partners_to_partner_indices, 

61) 

62from meshpy.utils.environment import is_testing as _is_testing 

63from meshpy.utils.nodes import filter_nodes as _filter_nodes 

64from meshpy.utils.nodes import find_close_nodes as _find_close_nodes 

65from meshpy.utils.nodes import get_min_max_nodes as _get_min_max_nodes 

66from meshpy.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates 

67from meshpy.utils.nodes import get_nodal_quaternions as _get_nodal_quaternions 

68from meshpy.utils.nodes import get_nodes_by_function as _get_nodes_by_function 

69 

70 

71class Mesh: 

72 """A class that contains a full mesh, i.e. Nodes, Elements, Boundary 

73 Conditions, Sets, Couplings, Materials and Functions.""" 

74 

75 def __init__(self): 

76 """Initialize all empty containers.""" 

77 

78 self.nodes = [] 

79 self.elements = [] 

80 self.materials = [] 

81 self.functions = [] 

82 self.geometry_sets = _GeometrySetContainer() 

83 self.boundary_conditions = _BoundaryConditionContainer() 

84 

85 @staticmethod 

86 def get_base_mesh_item_type(item): 

87 """Return the base mesh type of the given item. 

88 

89 Amongst other things, we need this function so we can check if 

90 all items of a list are of the same "base" type. 

91 

92 Args: 

93 item: The object we want to get the base type from. 

94 """ 

95 

96 for cls in ( 

97 Mesh, 

98 _Function, 

99 _BoundaryConditionBase, 

100 _Material, 

101 _Node, 

102 _Element, 

103 _GeometrySetBase, 

104 _GeometryName, 

105 ): 

106 if isinstance(item, cls): 

107 return cls 

108 return type(item) 

109 

110 def add(self, *args, **kwargs): 

111 """Add an item to this mesh, depending on its type. 

112 

113 If an list is given each list element is added with this 

114 function. If multiple arguments are given, each one is 

115 individually added with this function. Keyword arguments are 

116 passed through to the adding function. 

117 """ 

118 

119 match len(args): 

120 case 0: 

121 raise ValueError("At least one argument is required!") 

122 case 1: 

123 add_item = args[0] 

124 base_type = self.get_base_mesh_item_type(add_item) 

125 

126 base_type_to_method_map = { 

127 Mesh: self.add_mesh, 

128 _Function: self.add_function, 

129 _BoundaryConditionBase: self.add_bc, 

130 _Material: self.add_material, 

131 _Node: self.add_node, 

132 _Element: self.add_element, 

133 _GeometrySetBase: self.add_geometry_set, 

134 _GeometryName: self.add_geometry_name, 

135 list: self.add_list, 

136 } 

137 if base_type in base_type_to_method_map: 

138 base_type_to_method_map[base_type](add_item, **kwargs) 

139 else: 

140 raise TypeError( 

141 f'No Mesh.add case implemented for type: "{type(add_item)}" with base type "{base_type}"!' 

142 ) 

143 case _: 

144 for item in args: 

145 self.add(item, **kwargs) 

146 

147 def add_mesh(self, mesh): 

148 """Add the content of another mesh to this mesh.""" 

149 

150 # Add each item from mesh to self. 

151 self.add(mesh.nodes) 

152 self.add(mesh.elements) 

153 self.add(mesh.materials) 

154 self.add(mesh.functions) 

155 self.geometry_sets.extend(mesh.geometry_sets) 

156 self.boundary_conditions.extend(mesh.boundary_conditions) 

157 

158 def add_bc(self, bc): 

159 """Add a boundary condition to this mesh.""" 

160 bc_key = bc.bc_type 

161 geom_key = bc.geometry_set.geometry_type 

162 bc.geometry_set.check_replaced_nodes() 

163 self.boundary_conditions.append((bc_key, geom_key), bc) 

164 

165 def add_function(self, function): 

166 """Add a function to this mesh item. 

167 

168 Check that the function is only added once. 

169 """ 

170 if function not in self.functions: 

171 self.functions.append(function) 

172 

173 def add_material(self, material): 

174 """Add a material to this mesh item. 

175 

176 Check that the material is only added once. 

177 """ 

178 if material not in self.materials: 

179 self.materials.append(material) 

180 

181 def add_node(self, node): 

182 """Add a node to this mesh.""" 

183 if node in self.nodes: 

184 raise ValueError("The node is already in this mesh!") 

185 self.nodes.append(node) 

186 

187 def add_element(self, element): 

188 """Add an element to this mesh.""" 

189 if element in self.elements: 

190 raise ValueError("The element is already in this mesh!") 

191 self.elements.append(element) 

192 

193 def add_geometry_set(self, geometry_set): 

194 """Add a geometry set to this mesh.""" 

195 geometry_set.check_replaced_nodes() 

196 self.geometry_sets.append(geometry_set.geometry_type, geometry_set) 

197 

198 def add_geometry_name(self, geometry_name): 

199 """Add a set of geometry sets to this mesh. 

200 

201 Sort by the keys here to create a deterministic ordering, 

202 especially for testing purposes 

203 """ 

204 keys = list(geometry_name.keys()) 

205 keys.sort() 

206 for key in keys: 

207 self.add(geometry_name[key]) 

208 

209 def add_list(self, add_list: _List, **kwargs) -> None: 

210 """Add a list of items to this mesh. 

211 

212 Args: 

213 add_list: 

214 List to be added to the mesh. This method checks that all 

215 base types of the list items are the same. 

216 

217 In the special case of a node or element list, we add the whole list 

218 at once. This avoids a check for duplicate entries for every 

219 addition, which scales very badly. By doing it this way we only 

220 check the final list for duplicate entries which is much more 

221 performant. 

222 

223 For all other types of items, we add each element individually 

224 via the Mesh.add method. 

225 """ 

226 

227 types = {self.get_base_mesh_item_type(item) for item in add_list} 

228 if len(types) > 1: 

229 raise TypeError( 

230 f"You can only add lists with the same type of element. Got {types}" 

231 ) 

232 elif len(types) == 1: 

233 list_type = types.pop() 

234 

235 def extend_internal_list(self_list: _List, new_list: _List) -> None: 

236 """Extend an internal list with the new list. 

237 

238 It is checked that the final list does not have 

239 duplicate entries. 

240 """ 

241 self_list.extend(new_list) 

242 if not len(set(self_list)) == len(self_list): 

243 raise ValueError( 

244 "The added list contains entries already existing in the Mesh" 

245 ) 

246 

247 if list_type == _Node: 

248 extend_internal_list(self.nodes, add_list) 

249 elif list_type == _Element: 

250 extend_internal_list(self.elements, add_list) 

251 else: 

252 for item in add_list: 

253 self.add(item, **kwargs) 

254 

255 def replace_node(self, old_node, new_node): 

256 """Replace the first node with the second node.""" 

257 

258 # Check that the new node is in mesh. 

259 if new_node not in self.nodes: 

260 raise ValueError("The new node is not in the mesh!") 

261 

262 for i, node in enumerate(self.nodes): 

263 if node == old_node: 

264 del self.nodes[i] 

265 break 

266 else: 

267 raise ValueError("The node that should be replaced is not in the mesh") 

268 

269 def get_unique_geometry_sets( 

270 self, 

271 *, 

272 coupling_sets: bool = True, 

273 link_to_nodes: str = "no_link", 

274 geometry_set_start_indices: _Optional[_Dict] = None, 

275 ): 

276 """Return a geometry set container that contains geometry sets 

277 explicitly added to the mesh, as well as sets for boundary conditions. 

278 

279 The i_global values are set in the returned geometry sets. 

280 

281 Args: 

282 coupling_sets: 

283 If this is true, also sets for couplings will be added. 

284 link_to_nodes: 

285 "no_link": 

286 No link between the geometry set and the nodes is set 

287 "explicitly_contained_nodes": 

288 A link will be set for all nodes that are explicitly part of the geometry set 

289 "all_nodes": 

290 A link will be set for all nodes that are part of the geometry set, i.e., also 

291 nodes connected to elements of an element set. This is mainly used for vtk 

292 output so we can color the nodes which are part of element sets. 

293 geometry_set_start_indices: 

294 Dictionary where the keys are geometry type and the value is the starting 

295 index for those geometry sets. This can be used to define an offset in the 

296 i_global numbering of the geometry sets. The offsets default to 0. 

297 """ 

298 

299 is_link_nodes = not link_to_nodes == "no_link" 

300 if is_link_nodes: 

301 # First clear all links in existing nodes. 

302 for node in self.nodes: 

303 node.node_sets_link = [] 

304 

305 # Make a copy of the sets in this mesh. 

306 mesh_sets = self.geometry_sets.copy() 

307 

308 # Add sets from boundary conditions. 

309 for (bc_key, geom_key), bc_list in self.boundary_conditions.items(): 

310 for bc in bc_list: 

311 # Check if sets from couplings should be added. 

312 is_coupling = bc_key in ( 

313 _mpy.bc.point_coupling, 

314 bc_key == _mpy.bc.point_coupling_penalty, 

315 ) 

316 if (is_coupling and coupling_sets) or (not is_coupling): 

317 # Only add set if it is not already in the container. 

318 # For example if multiple Neumann boundary conditions 

319 # are applied on the same node set. 

320 if bc.geometry_set not in mesh_sets[geom_key]: 

321 mesh_sets[geom_key].append(bc.geometry_set) 

322 

323 for key in mesh_sets.keys(): 

324 i_global_offset = 0 

325 if geometry_set_start_indices is not None: 

326 if key in geometry_set_start_indices: 

327 i_global_offset = geometry_set_start_indices[key] 

328 else: 

329 raise KeyError("Could not find {key} in geometry_set_start_indices") 

330 for i, geometry_set in enumerate(mesh_sets[key]): 

331 # Add global indices to the geometry set. 

332 geometry_set.i_global = i + 1 + i_global_offset 

333 if is_link_nodes: 

334 geometry_set.link_to_nodes(link_to_nodes=link_to_nodes) 

335 

336 return mesh_sets 

337 

338 def set_node_links(self): 

339 """Create a link of all elements to the nodes connected to them. 

340 

341 Also add a link to this mesh. 

342 """ 

343 for element in self.elements: 

344 for node in element.nodes: 

345 node.element_link.append(element) 

346 for node in self.nodes: 

347 node.mesh = self 

348 

349 def translate(self, vector): 

350 """Translate all beam nodes of this mesh. 

351 

352 Args 

353 ---- 

354 vector: _np.array, list 

355 3D vector that will be added to all nodes. 

356 """ 

357 for node in self.nodes: 

358 node.coordinates += vector 

359 

360 def rotate(self, rotation, origin=None, only_rotate_triads=False): 

361 """Rotate all beam nodes of the mesh with rotation. 

362 

363 Args 

364 ---- 

365 rotation: Rotation, list(quaternions) (nx4) 

366 The rotation that will be applied to the nodes. Can also be an 

367 array with a quaternion for each node. 

368 origin: 3D vector 

369 If this is given, the mesh is rotated about this point. Default is 

370 (0,0,0) 

371 only_rotate_triads: bool 

372 If true the nodal positions are not changed. 

373 """ 

374 

375 # Get array with all quaternions for the nodes. 

376 rot1 = _get_nodal_quaternions(self.nodes) 

377 

378 # Apply the rotation to the rotation of all nodes. 

379 rot_new = _add_rotations(rotation, rot1) 

380 

381 if not only_rotate_triads: 

382 # Get array with all positions for the nodes. 

383 pos = _get_nodal_coordinates(self.nodes) 

384 pos_new = _rotate_coordinates(pos, rotation, origin=origin) 

385 

386 for i, node in enumerate(self.nodes): 

387 if isinstance(node, _NodeCosserat): 

388 node.rotation.q = rot_new[i, :] 

389 if not only_rotate_triads: 

390 node.coordinates = pos_new[i, :] 

391 

392 def reflect(self, normal_vector, origin=None, flip_beams=False): 

393 """Reflect all nodes of the mesh with respect to a plane defined by its 

394 normal_vector. Per default the plane goes through the origin, if not a 

395 point on the plane can be given with the parameter origin. 

396 

397 For the reflection we assume that e1' and e2' are mirrored with respect 

398 to the original frame and e3' is in the opposite direction than the 

399 mirrored e3. 

400 

401 With the defined mirroring strategy, the quaternion to be applied on 

402 the existing rotations can be calculated the following way: 

403 q[0] = e3 * n 

404 q[1,2,3] = e3 x n 

405 This constructs a rotation with the rotation axis on the plane, and 

406 normal to the vector e3. The rotation angle is twice the angle of e3 

407 to n. 

408 

409 Args 

410 ---- 

411 normal_3D vector 

412 The normal vector of the reflection plane. 

413 origin: 3D vector 

414 Per default the reflection plane goes through the origin. If this 

415 parameter is given, the point is on the plane. 

416 flip_beams: bool 

417 When True, the beams are flipped, so that the direction along the 

418 beam is reversed. 

419 """ 

420 

421 # Normalize the normal vector. 

422 normal_vector = _np.asarray(normal_vector) / _np.linalg.norm(normal_vector) 

423 

424 # Get array with all quaternions and positions for the nodes. 

425 pos = _get_nodal_coordinates(self.nodes) 

426 rot1 = _get_nodal_quaternions(self.nodes) 

427 

428 # Check if origin has to be added. 

429 if origin is not None: 

430 pos -= origin 

431 

432 # Get the reflection matrix A. 

433 A = _np.eye(3) - 2.0 * _np.outer(normal_vector, normal_vector) 

434 

435 # Calculate the new positions. 

436 pos_new = _np.dot(pos, A) 

437 

438 # Move back from the origin. 

439 if origin is not None: 

440 pos_new += origin 

441 

442 # First get all e3 vectors of the nodes. 

443 e3 = _np.zeros_like(pos) 

444 e3[:, 0] = 2 * (rot1[:, 0] * rot1[:, 2] + rot1[:, 1] * rot1[:, 3]) 

445 e3[:, 1] = 2 * (-1 * rot1[:, 0] * rot1[:, 1] + rot1[:, 2] * rot1[:, 3]) 

446 e3[:, 2] = rot1[:, 0] ** 2 - rot1[:, 1] ** 2 - rot1[:, 2] ** 2 + rot1[:, 3] ** 2 

447 

448 # Get the dot and cross product of e3 and the normal vector. 

449 rot2 = _np.zeros_like(rot1) 

450 rot2[:, 0] = _np.dot(e3, normal_vector) 

451 rot2[:, 1:] = _np.cross(e3, normal_vector) 

452 

453 # Add to the existing rotations. 

454 rot_new = _add_rotations(rot2, rot1) 

455 

456 if flip_beams: 

457 # To achieve the flip, the triads are rotated with the angle pi 

458 # around the e2 axis. 

459 rot_flip = _Rotation([0, 1, 0], _np.pi) 

460 rot_new = _add_rotations(rot_new, rot_flip) 

461 

462 # For solid elements we need to adapt the connectivity to avoid negative Jacobians. 

463 # For beam elements this is optional. 

464 for element in self.elements: 

465 if isinstance(element, _Beam): 

466 if flip_beams: 

467 element.flip() 

468 else: 

469 element.flip() 

470 

471 # Set the new positions and rotations. 

472 for i, node in enumerate(self.nodes): 

473 node.coordinates = pos_new[i, :] 

474 if isinstance(node, _NodeCosserat): 

475 node.rotation.q = rot_new[i, :] 

476 

477 def wrap_around_cylinder(self, radius=None, advanced_warning=True): 

478 """Wrap the geometry around a cylinder. The y-z plane gets morphed into 

479 the z-axis of symmetry. If all nodes are on the same y-z plane, the 

480 radius of the created cylinder is the x coordinate of that plane. If 

481 the nodes are not on the same y-z plane, the radius has to be given 

482 explicitly. 

483 

484 Args 

485 ---- 

486 radius: double 

487 If this value is given AND not all nodes are on the same y-z plane, 

488 then use this radius for the calculation of phi for all nodes. 

489 This might still lead to distorted elements!. 

490 advanced_warning: bool 

491 If each element should be checked if it is either parallel to the 

492 y-z or x-z plane. This is computationally expensive, but in most 

493 cases (up to 100,000 elements) this check can be left activated. 

494 """ 

495 

496 pos = _get_nodal_coordinates(self.nodes) 

497 quaternions = _np.zeros([len(self.nodes), 4]) 

498 

499 # The x coordinate is the radius, the y coordinate the arc length. 

500 points_x = pos[:, 0].copy() 

501 

502 # Check if all points are on the same y-z plane. 

503 if _np.abs(_np.min(points_x) - _np.max(points_x)) > _mpy.eps_pos: 

504 # The points are not all on the y-z plane, get the reference 

505 # radius. 

506 if radius is not None: 

507 if advanced_warning: 

508 # Here we check, if each element lays on a plane parallel 

509 # to the y-z plane, or parallel to the x-z plane. 

510 # 

511 # To be exactly sure, we could check the rotations here, 

512 # i.e. if they are also in plane. 

513 element_warning = [] 

514 for i_element, element in enumerate(self.elements): 

515 element_coordinates = _np.zeros([len(element.nodes), 3]) 

516 for i_node, node in enumerate(element.nodes): 

517 element_coordinates[i_node, :] = node.coordinates 

518 is_yz = ( 

519 _np.max( 

520 _np.abs( 

521 element_coordinates[:, 0] 

522 - element_coordinates[0, 0] 

523 ) 

524 ) 

525 < _mpy.eps_pos 

526 ) 

527 is_xz = ( 

528 _np.max( 

529 _np.abs( 

530 element_coordinates[:, 1] 

531 - element_coordinates[0, 1] 

532 ) 

533 ) 

534 < _mpy.eps_pos 

535 ) 

536 if not (is_yz or is_xz): 

537 element_warning.append(i_element) 

538 if len(element_warning) != 0: 

539 _warnings.warn( 

540 "There are elements which are not " 

541 "parallel to the y-z or x-y plane. This will lead " 

542 "to distorted elements!" 

543 ) 

544 else: 

545 _warnings.warn( 

546 "The nodes are not on the same y-z plane. " 

547 "This may lead to distorted elements!" 

548 ) 

549 else: 

550 raise ValueError( 

551 "The nodes that should be wrapped around a " 

552 "cylinder are not on the same y-z plane. This will give " 

553 "unexpected results. Give a reference radius!" 

554 ) 

555 radius_phi = radius 

556 radius_points = points_x 

557 elif radius is None or _np.abs(points_x[0] - radius) < _mpy.eps_pos: 

558 radius_points = radius_phi = points_x[0] 

559 else: 

560 raise ValueError( 

561 ( 

562 "The points are all on the same y-z plane with " 

563 "the x-coordinate {} but the given radius {} is different. " 

564 "This does not make sense." 

565 ).format(points_x[0], radius) 

566 ) 

567 

568 # Get the angle for all nodes. 

569 phi = pos[:, 1] / radius_phi 

570 

571 # The rotation is about the z-axis. 

572 quaternions[:, 0] = _np.cos(0.5 * phi) 

573 quaternions[:, 3] = _np.sin(0.5 * phi) 

574 

575 # Set the new positions in the global array. 

576 pos[:, 0] = radius_points * _np.cos(phi) 

577 pos[:, 1] = radius_points * _np.sin(phi) 

578 

579 # Rotate the mesh 

580 self.rotate(quaternions, only_rotate_triads=True) 

581 

582 # Set the new position for the nodes. 

583 for i, node in enumerate(self.nodes): 

584 node.coordinates = pos[i, :] 

585 

586 def couple_nodes( 

587 self, 

588 *, 

589 nodes=None, 

590 reuse_matching_nodes=False, 

591 coupling_type=_mpy.bc.point_coupling, 

592 coupling_dof_type=_mpy.coupling_dof.fix, 

593 ): 

594 """Search through nodes and connect all nodes with the same 

595 coordinates. 

596 

597 Args: 

598 ---- 

599 nodes: [Node] 

600 List of nodes to couple. If None is given, all nodes of the mesh 

601 are coupled (except middle nodes). 

602 reuse_matching_nodes: bool 

603 If two nodes have the same position and rotation, the nodes are 

604 reduced to one node in the mesh. Be aware, that this might lead to 

605 issues if not all DOFs of the nodes should be coupled. 

606 coupling_type: mpy.bc 

607 Type of point coupling. 

608 coupling_dof_type: str, mpy.coupling_dof 

609 str: The string that will be used in the input file. 

610 mpy.coupling_dof.fix: Fix all positional and rotational DOFs of the 

611 nodes together. 

612 mpy.coupling_dof.joint: Fix all positional DOFs of the nodes 

613 together. 

614 """ 

615 

616 # Check that a coupling BC is given. 

617 if coupling_type not in ( 

618 _mpy.bc.point_coupling, 

619 _mpy.bc.point_coupling_penalty, 

620 ): 

621 raise ValueError( 

622 "Only coupling conditions can be applied in 'couple_nodes'!" 

623 ) 

624 

625 # Get the nodes that should be checked for coupling. Middle nodes are 

626 # not checked, as coupling can only be applied to the boundary nodes. 

627 if nodes is None: 

628 node_list = self.nodes 

629 else: 

630 node_list = nodes 

631 node_list = _filter_nodes(node_list, middle_nodes=False) 

632 partner_nodes = _find_close_nodes(node_list) 

633 if len(partner_nodes) == 0: 

634 # If no partner nodes were found, end this function. 

635 return 

636 

637 if reuse_matching_nodes: 

638 # Check if there are nodes with the same rotation. If there are the 

639 # nodes are reused, and no coupling is inserted. 

640 

641 # Set the links to all nodes in the mesh. In this case we have to use 

642 # "all_nodes" since we also have to replace nodes that are in existing 

643 # GeometrySetNodes. 

644 self.unlink_nodes() 

645 self.get_unique_geometry_sets(link_to_nodes="explicitly_contained_nodes") 

646 self.set_node_links() 

647 

648 # Go through partner nodes. 

649 for node_list in partner_nodes: 

650 # Get array with rotation vectors. 

651 rotation_vectors = _np.zeros([len(node_list), 3]) 

652 for i, node in enumerate(node_list): 

653 if isinstance(node, _NodeCosserat): 

654 rotation_vectors[i, :] = node.rotation.get_rotation_vector() 

655 else: 

656 # For the case of nodes that belong to solid elements, 

657 # we define the following default value: 

658 rotation_vectors[i, :] = [4 * _np.pi, 0, 0] 

659 

660 # Use find close points function to find nodes with the 

661 # same rotation. 

662 partners, n_partners = _find_close_points( 

663 rotation_vectors, tol=_mpy.eps_quaternion 

664 ) 

665 

666 # Check if nodes with the same rotations were found. 

667 if n_partners == 0: 

668 self.add( 

669 _coupling_factory(node_list, coupling_type, coupling_dof_type) 

670 ) 

671 else: 

672 # There are nodes that need to be combined. 

673 combining_nodes = [] 

674 coupling_nodes = [] 

675 found_partner_id = [None for _i in range(n_partners)] 

676 

677 # Add the nodes that need to be combined and add the nodes 

678 # that will be coupled. 

679 for i, partner in enumerate(partners): 

680 if partner == -1: 

681 # This node does not have a partner with the same 

682 # rotation. 

683 coupling_nodes.append(node_list[i]) 

684 

685 elif found_partner_id[partner] is not None: 

686 # This node has already a processed partner, add 

687 # this one to the combining nodes. 

688 combining_nodes[found_partner_id[partner]].append( 

689 node_list[i] 

690 ) 

691 

692 else: 

693 # This is the first node of a partner set that was 

694 # found. This one will remain, the other ones will 

695 # be replaced with this one. 

696 new_index = len(combining_nodes) 

697 found_partner_id[partner] = new_index 

698 combining_nodes.append([node_list[i]]) 

699 coupling_nodes.append(node_list[i]) 

700 

701 # Add the coupling nodes. 

702 if len(coupling_nodes) > 1: 

703 self.add( 

704 _coupling_factory( 

705 coupling_nodes, coupling_type, coupling_dof_type 

706 ) 

707 ) 

708 

709 # Replace the identical nodes. 

710 for combine_list in combining_nodes: 

711 master_node = combine_list[0] 

712 for node in combine_list[1:]: 

713 node.replace_with(master_node) 

714 

715 else: 

716 # Connect close nodes with a coupling. 

717 for node_list in partner_nodes: 

718 self.add(_coupling_factory(node_list, coupling_type, coupling_dof_type)) 

719 

720 def unlink_nodes(self): 

721 """Delete the linked arrays and global indices in all nodes.""" 

722 for node in self.nodes: 

723 node.unlink() 

724 

725 def get_nodes_by_function(self, *args, **kwargs): 

726 """Return all nodes for which the function evaluates to true.""" 

727 return _get_nodes_by_function(self.nodes, *args, **kwargs) 

728 

729 def get_min_max_nodes(self, *args, **kwargs): 

730 """Return a geometry set with the max and min nodes in all 

731 directions.""" 

732 return _get_min_max_nodes(self.nodes, *args, **kwargs) 

733 

734 def check_overlapping_elements(self, raise_error=True): 

735 """Check if there are overlapping elements in the mesh. 

736 

737 This is done by checking if all middle nodes of beam elements 

738 have unique coordinates in the mesh. 

739 """ 

740 

741 # Number of middle nodes. 

742 middle_nodes = [node for node in self.nodes if node.is_middle_node] 

743 

744 # Only check if there are middle nodes. 

745 if len(middle_nodes) == 0: 

746 return 

747 

748 # Get array with middle nodes. 

749 coordinates = _np.zeros([len(middle_nodes), 3]) 

750 for i, node in enumerate(middle_nodes): 

751 coordinates[i, :] = node.coordinates 

752 

753 # Check if there are double entries in the coordinates. 

754 has_partner, partner = _find_close_points(coordinates) 

755 partner_indices = _point_partners_to_partner_indices(has_partner, partner) 

756 if partner > 0: 

757 if raise_error: 

758 raise ValueError( 

759 "There are multiple middle nodes with the " 

760 "same coordinates. Per default this raises an error! " 

761 "This check can be turned of with " 

762 "mpy.check_overlapping_elements=False" 

763 ) 

764 else: 

765 _warnings.warn( 

766 "There are multiple middle nodes with the same coordinates!" 

767 ) 

768 

769 # Add the partner index to the middle nodes. 

770 for i_partner, partners in enumerate(partner_indices): 

771 for i_node in partners: 

772 middle_nodes[i_node].element_partner_index = i_partner 

773 

774 def get_vtk_representation( 

775 self, *, overlapping_elements=True, coupling_sets=False, **kwargs 

776 ): 

777 """Return a vtk representation of the beams and solid in this mesh. 

778 

779 Args 

780 ---- 

781 overlapping_elements: bool 

782 I elements should be checked for overlapping. If they overlap, the 

783 output will mark them. 

784 coupling_sets: bool 

785 If coupling sets should also be displayed. 

786 """ 

787 

788 # Object to store VKT data (and write it to file) 

789 vtk_writer_beam = _VTKWriter() 

790 vtk_writer_solid = _VTKWriter() 

791 

792 # Get the set numbers of the mesh 

793 mesh_sets = self.get_unique_geometry_sets( 

794 coupling_sets=coupling_sets, link_to_nodes="all_nodes" 

795 ) 

796 

797 # Set the global value for digits in the VTK output. 

798 # Get highest number of node_sets. 

799 max_sets = max(len(geometry_list) for geometry_list in mesh_sets.values()) 

800 

801 # Set the mpy value. 

802 digits = len(str(max_sets)) 

803 _mpy.vtk_node_set_format = "{:0" + str(digits) + "}" 

804 

805 if overlapping_elements: 

806 # Check for overlapping elements. 

807 self.check_overlapping_elements(raise_error=False) 

808 

809 # Get representation of elements. 

810 for element in self.elements: 

811 element.get_vtk(vtk_writer_beam, vtk_writer_solid, **kwargs) 

812 

813 # Finish and return the writers 

814 vtk_writer_beam.complete_data() 

815 vtk_writer_solid.complete_data() 

816 return vtk_writer_beam, vtk_writer_solid 

817 

818 def write_vtk( 

819 self, output_name="meshpy", output_directory="", binary=True, **kwargs 

820 ): 

821 """Write the contents of this mesh to VTK files. 

822 

823 Args 

824 ---- 

825 output_name: str 

826 Base name of the output file. There will be a {}_beam.vtu and 

827 {}_solid.vtu file. 

828 output_directory: path 

829 Directory where the output files will be written. 

830 binary: bool 

831 If the data should be written encoded in binary or in human readable text 

832 

833 **kwargs 

834 For all of them look into: 

835 Mesh().get_vtk_representation 

836 Beam().get_vtk 

837 VolumeElement().get_vtk 

838 ---- 

839 beam_centerline_visualization_segments: int 

840 Number of segments to be used for visualization of beam centerline between successive 

841 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For 

842 Values greater than 1, a Hermite interpolation of the centerline is assumed for 

843 visualization purposes. 

844 """ 

845 

846 vtk_writer_beam, vtk_writer_solid = self.get_vtk_representation(**kwargs) 

847 

848 # Write to file, only if there is at least one point in the writer. 

849 if vtk_writer_beam.points.GetNumberOfPoints() > 0: 

850 filepath = _os.path.join(output_directory, output_name + "_beam.vtu") 

851 vtk_writer_beam.write_vtk(filepath, binary=binary) 

852 if vtk_writer_solid.points.GetNumberOfPoints() > 0: 

853 filepath = _os.path.join(output_directory, output_name + "_solid.vtu") 

854 vtk_writer_solid.write_vtk(filepath, binary=binary) 

855 

856 def display_pyvista( 

857 self, 

858 *, 

859 beam_nodes=True, 

860 beam_tube=True, 

861 beam_cross_section_directors=True, 

862 beam_radius_for_display=None, 

863 resolution=20, 

864 parallel_projection=False, 

865 **kwargs, 

866 ): 

867 """Display the mesh in pyvista. 

868 

869 If this is called in a GitHub testing run, nothing will be shown, instead 

870 the _pv.plotter object will be returned. 

871 

872 Args 

873 ---- 

874 beam_nodes: bool 

875 If the beam nodes should be displayed. The start and end nodes of each 

876 beam will be shown in green, possible middle nodes inside the element 

877 are shown in cyan. 

878 beam_tube: bool 

879 If the beam should be rendered as a tube 

880 beam_cross_section_directors: bool 

881 If the cross section directors should be displayed (at each node) 

882 beam_radius_for_display: float 

883 If not all beams have an explicitly given radius (in the material 

884 definition) this value will be used to approximate the beams radius 

885 for visualization 

886 resolution: int 

887 Indicates how many triangulations will be performed to visualize arrows, 

888 tubes and spheres. 

889 parallel_projection: bool 

890 Flag to change camera view to parallel projection. 

891 

892 **kwargs 

893 For all of them look into: 

894 Mesh().get_vtk_representation 

895 Beam().get_vtk 

896 VolumeElement().get_vtk 

897 ---- 

898 beam_centerline_visualization_segments: int 

899 Number of segments to be used for visualization of beam centerline between successive 

900 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For 

901 Values greater than 1, a Hermite interpolation of the centerline is assumed for 

902 visualization purposes. 

903 """ 

904 

905 plotter = _pv.Plotter() 

906 plotter.renderer.add_axes() 

907 

908 if parallel_projection: 

909 plotter.enable_parallel_projection() 

910 

911 vtk_writer_beam, vtk_writer_solid = self.get_vtk_representation(**kwargs) 

912 

913 if vtk_writer_beam.points.GetNumberOfPoints() > 0: 

914 beam_grid = _pv.UnstructuredGrid(vtk_writer_beam.grid) 

915 

916 # Check if all beams have a given cross-section radius, if not set the given input 

917 # value 

918 all_beams_have_cross_section_radius = (