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

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 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 

34import quaternion as _quaternion 

35from numpy.typing import NDArray as _NDArray 

36 

37from beamme.core.boundary_condition import ( 

38 BoundaryConditionBase as _BoundaryConditionBase, 

39) 

40from beamme.core.boundary_condition import ( 

41 BoundaryConditionContainer as _BoundaryConditionContainer, 

42) 

43from beamme.core.conf import bme as _bme 

44from beamme.core.coupling import coupling_factory as _coupling_factory 

45from beamme.core.element import Element as _Element 

46from beamme.core.element_beam import Beam as _Beam 

47from beamme.core.function import Function as _Function 

48from beamme.core.geometry_set import GeometryName as _GeometryName 

49from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase 

50from beamme.core.geometry_set import GeometrySetContainer as _GeometrySetContainer 

51from beamme.core.material import Material as _Material 

52from beamme.core.node import Node as _Node 

53from beamme.core.node import NodeCosserat as _NodeCosserat 

54from beamme.core.rotation import Rotation as _Rotation 

55from beamme.core.rotation import add_rotations as _add_rotations 

56from beamme.core.rotation import rotate_coordinates as _rotate_coordinates 

57from beamme.core.vtk_writer import VTKWriter as _VTKWriter 

58from beamme.geometric_search.find_close_points import ( 

59 find_close_points as _find_close_points, 

60) 

61from beamme.geometric_search.find_close_points import ( 

62 point_partners_to_partner_indices as _point_partners_to_partner_indices, 

63) 

64from beamme.utils.environment import is_testing as _is_testing 

65from beamme.utils.nodes import filter_nodes as _filter_nodes 

66from beamme.utils.nodes import find_close_nodes as _find_close_nodes 

67from beamme.utils.nodes import get_min_max_nodes as _get_min_max_nodes 

68from beamme.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates 

69from beamme.utils.nodes import get_nodal_quaternions as _get_nodal_quaternions 

70from beamme.utils.nodes import get_nodes_by_function as _get_nodes_by_function 

71 

72 

73class Mesh: 

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

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

76 

77 def __init__(self): 

78 """Initialize all empty containers.""" 

79 

80 self.nodes = [] 

81 self.elements = [] 

82 self.materials = [] 

83 self.functions = [] 

84 self.geometry_sets = _GeometrySetContainer() 

85 self.boundary_conditions = _BoundaryConditionContainer() 

86 

87 @staticmethod 

88 def get_base_mesh_item_type(item): 

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

90 

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

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

93 

94 Args: 

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

96 """ 

97 

98 for cls in ( 

99 Mesh, 

100 _Function, 

101 _BoundaryConditionBase, 

102 _Material, 

103 _Node, 

104 _Element, 

105 _GeometrySetBase, 

106 _GeometryName, 

107 ): 

108 if isinstance(item, cls): 

109 return cls 

110 return type(item) 

111 

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

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

114 

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

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

117 individually added with this function. Keyword arguments are 

118 passed through to the adding function. 

119 """ 

120 

121 match len(args): 

122 case 0: 

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

124 case 1: 

125 add_item = args[0] 

126 base_type = self.get_base_mesh_item_type(add_item) 

127 

128 base_type_to_method_map = { 

129 Mesh: self.add_mesh, 

130 _Function: self.add_function, 

131 _BoundaryConditionBase: self.add_bc, 

132 _Material: self.add_material, 

133 _Node: self.add_node, 

134 _Element: self.add_element, 

135 _GeometrySetBase: self.add_geometry_set, 

136 _GeometryName: self.add_geometry_name, 

137 list: self.add_list, 

138 } 

139 if base_type in base_type_to_method_map: 

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

141 else: 

142 raise TypeError( 

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

144 ) 

145 case _: 

146 for item in args: 

147 self.add(item, **kwargs) 

148 

149 def add_mesh(self, mesh): 

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

151 

152 # Add each item from mesh to self. 

153 self.add(mesh.nodes) 

154 self.add(mesh.elements) 

155 self.add(mesh.materials) 

156 self.add(mesh.functions) 

157 self.geometry_sets.extend(mesh.geometry_sets) 

158 self.boundary_conditions.extend(mesh.boundary_conditions) 

159 

160 def add_bc(self, bc): 

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

162 bc_key = bc.bc_type 

163 geom_key = bc.geometry_set.geometry_type 

164 bc.geometry_set.check_replaced_nodes() 

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

166 

167 def add_function(self, function): 

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

169 

170 Check that the function is only added once. 

171 """ 

172 if function not in self.functions: 

173 self.functions.append(function) 

174 

175 def add_material(self, material): 

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

177 

178 Check that the material is only added once. 

179 """ 

180 if material not in self.materials: 

181 self.materials.append(material) 

182 

183 def add_node(self, node): 

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

185 if node in self.nodes: 

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

187 self.nodes.append(node) 

188 

189 def add_element(self, element): 

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

191 if element in self.elements: 

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

193 self.elements.append(element) 

194 

195 def add_geometry_set(self, geometry_set): 

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

197 geometry_set.check_replaced_nodes() 

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

199 

200 def add_geometry_name(self, geometry_name): 

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

202 

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

204 especially for testing purposes 

205 """ 

206 keys = list(geometry_name.keys()) 

207 keys.sort() 

208 for key in keys: 

209 self.add(geometry_name[key]) 

210 

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

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

213 

214 Args: 

215 add_list: 

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

217 base types of the list items are the same. 

218 

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

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

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

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

223 performant. 

224 

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

226 via the Mesh.add method. 

227 """ 

228 

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

230 if len(types) > 1: 

231 raise TypeError( 

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

233 ) 

234 elif len(types) == 1: 

235 list_type = types.pop() 

236 

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

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

239 

240 It is checked that the final list does not have 

241 duplicate entries. 

242 """ 

243 self_list.extend(new_list) 

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

245 raise ValueError( 

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

247 ) 

248 

249 if list_type == _Node: 

250 extend_internal_list(self.nodes, add_list) 

251 elif list_type == _Element: 

252 extend_internal_list(self.elements, add_list) 

253 else: 

254 for item in add_list: 

255 self.add(item, **kwargs) 

256 

257 def replace_node(self, old_node, new_node): 

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

259 

260 # Check that the new node is in mesh. 

261 if new_node not in self.nodes: 

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

263 

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

265 if node == old_node: 

266 del self.nodes[i] 

267 break 

268 else: 

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

270 

271 def get_unique_geometry_sets( 

272 self, 

273 *, 

274 coupling_sets: bool = True, 

275 link_to_nodes: str = "no_link", 

276 geometry_set_start_indices: _Optional[_Dict] = None, 

277 ): 

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

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

280 

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

282 

283 Args: 

284 coupling_sets: 

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

286 link_to_nodes: 

287 "no_link": 

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

289 "explicitly_contained_nodes": 

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

291 "all_nodes": 

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

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

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

295 geometry_set_start_indices: 

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

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

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

299 """ 

300 

301 is_link_nodes = not link_to_nodes == "no_link" 

302 if is_link_nodes: 

303 # First clear all links in existing nodes. 

304 for node in self.nodes: 

305 node.node_sets_link = [] 

306 

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

308 mesh_sets = self.geometry_sets.copy() 

309 

310 # Add sets from boundary conditions. 

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

312 for bc in bc_list: 

313 # Check if sets from couplings should be added. 

314 is_coupling = bc_key in ( 

315 _bme.bc.point_coupling, 

316 bc_key == _bme.bc.point_coupling_penalty, 

317 ) 

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

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

320 # For example if multiple Neumann boundary conditions 

321 # are applied on the same node set. 

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

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

324 

325 for key in mesh_sets.keys(): 

326 i_global_offset = 0 

327 if geometry_set_start_indices is not None: 

328 if key in geometry_set_start_indices: 

329 i_global_offset = geometry_set_start_indices[key] 

330 else: 

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

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

333 # Add global indices to the geometry set. 

334 geometry_set.i_global = i + i_global_offset 

335 if is_link_nodes: 

336 geometry_set.link_to_nodes(link_to_nodes=link_to_nodes) 

337 

338 return mesh_sets 

339 

340 def set_node_links(self): 

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

342 

343 Also add a link to this mesh. 

344 """ 

345 for element in self.elements: 

346 for node in element.nodes: 

347 node.element_link.append(element) 

348 for node in self.nodes: 

349 node.mesh = self 

350 

351 def translate(self, vector): 

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

353 

354 Args 

355 ---- 

356 vector: _np.array, list 

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

358 """ 

359 for node in self.nodes: 

360 node.coordinates += vector 

361 

362 def rotate( 

363 self, 

364 rotation: _Rotation | _NDArray[_quaternion.quaternion], 

365 origin=None, 

366 only_rotate_triads: bool = False, 

367 ) -> None: 

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

369 

370 Args: 

371 rotation: The rotation(s) that will be applied to the nodes. If 

372 this is an array, it has to hold a quaternion for each node. 

373 origin (3D vector): If this is given, the mesh is rotated about 

374 this point. Defaults to (0, 0, 0). 

375 only_rotate_triads: If this is true, the nodal positions are not 

376 changed. 

377 """ 

378 

379 # Get array with all quaternions for the nodes. 

380 rot1 = _get_nodal_quaternions(self.nodes) 

381 

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

383 rot_new = _add_rotations(rotation, rot1) 

384 

385 if not only_rotate_triads: 

386 # Get array with all positions for the nodes. 

387 pos = _get_nodal_coordinates(self.nodes) 

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

389 

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

391 if isinstance(node, _NodeCosserat): 

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

393 if not only_rotate_triads: 

394 node.coordinates = pos_new[i, :] 

395 

396 def reflect(self, normal_vector, origin=None, flip_beams: bool = False) -> None: 

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

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

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

400 

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

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

403 mirrored e3. 

404 

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

406 the existing rotations can be calculated the following way: 

407 q[0] = e3 * n 

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

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

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

411 to n. 

412 

413 Args: 

414 normal_vector (3D vector): The normal vector of the reflection plane. 

415 origin (3D vector): The reflection plane goes through this point. 

416 Defaults to (0, 0, 0). 

417 flip_beams: When True, the beams are flipped, so that the direction 

418 along the 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( 

478 self, radius: float | None = None, advanced_warning: bool = True 

479 ) -> None: 

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

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

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

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

484 explicitly. 

485 

486 Args: 

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

488 plane, then use this radius for the calculation of phi for all 

489 nodes. This might still lead to distorted elements! 

490 advanced_warning: If each element should be checked if it is either parallel 

491 to the y-z or x-z plane. This is computationally expensive, but in most 

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

493 """ 

494 

495 pos = _get_nodal_coordinates(self.nodes) 

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

497 

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

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

500 

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

502 if _np.abs(_np.min(points_x) - _np.max(points_x)) > _bme.eps_pos: 

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

504 # radius. 

505 if radius is not None: 

506 if advanced_warning: 

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

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

509 # 

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

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

512 element_warning = [] 

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

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

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

516 element_coordinates[i_node, :] = node.coordinates 

517 is_yz = ( 

518 _np.max( 

519 _np.abs( 

520 element_coordinates[:, 0] 

521 - element_coordinates[0, 0] 

522 ) 

523 ) 

524 < _bme.eps_pos 

525 ) 

526 is_xz = ( 

527 _np.max( 

528 _np.abs( 

529 element_coordinates[:, 1] 

530 - element_coordinates[0, 1] 

531 ) 

532 ) 

533 < _bme.eps_pos 

534 ) 

535 if not (is_yz or is_xz): 

536 element_warning.append(i_element) 

537 if len(element_warning) != 0: 

538 _warnings.warn( 

539 "There are elements which are not " 

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

541 "to distorted elements!" 

542 ) 

543 else: 

544 _warnings.warn( 

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

546 "This may lead to distorted elements!" 

547 ) 

548 else: 

549 raise ValueError( 

550 "The nodes that should be wrapped around a " 

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

552 "unexpected results. Give a reference radius!" 

553 ) 

554 radius_phi = radius 

555 radius_points = points_x 

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

557 radius_points = radius_phi = points_x[0] 

558 else: 

559 raise ValueError( 

560 ( 

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

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

563 "This does not make sense." 

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

565 ) 

566 

567 # Get the angle for all nodes. 

568 phi = pos[:, 1] / radius_phi 

569 

570 # The rotation is about the z-axis. 

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

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

573 

574 # Set the new positions in the global array. 

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

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

577 

578 # Rotate the mesh 

579 self.rotate(quaternions, only_rotate_triads=True) 

580 

581 # Set the new position for the nodes. 

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

583 node.coordinates = pos[i, :] 

584 

585 def couple_nodes( 

586 self, 

587 *, 

588 nodes=None, 

589 reuse_matching_nodes=False, 

590 coupling_type=_bme.bc.point_coupling, 

591 coupling_dof_type=_bme.coupling_dof.fix, 

592 ): 

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

594 coordinates. 

595 

596 Args: 

597 ---- 

598 nodes: [Node] 

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

600 are coupled (except middle nodes). 

601 reuse_matching_nodes: bool 

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

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

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

605 coupling_type: bme.bc 

606 Type of point coupling. 

607 coupling_dof_type: str, bme.coupling_dof 

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

609 bme.coupling_dof.fix: Fix all positional and rotational DOFs of the 

610 nodes together. 

611 bme.coupling_dof.joint: Fix all positional DOFs of the nodes 

612 together. 

613 """ 

614 

615 # Check that a coupling BC is given. 

616 if coupling_type not in ( 

617 _bme.bc.point_coupling, 

618 _bme.bc.point_coupling_penalty, 

619 ): 

620 raise ValueError( 

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

622 ) 

623 

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

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

626 if nodes is None: 

627 node_list = self.nodes 

628 else: 

629 node_list = nodes 

630 node_list = _filter_nodes(node_list, middle_nodes=False) 

631 partner_nodes = _find_close_nodes(node_list) 

632 if len(partner_nodes) == 0: 

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

634 return 

635 

636 if reuse_matching_nodes: 

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

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

639 

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

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

642 # GeometrySetNodes. 

643 self.unlink_nodes() 

644 self.get_unique_geometry_sets(link_to_nodes="explicitly_contained_nodes") 

645 self.set_node_links() 

646 

647 # Go through partner nodes. 

648 for node_list in partner_nodes: 

649 # Get array with rotation vectors. 

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

651 for i, node in enumerate(node_list): 

652 if isinstance(node, _NodeCosserat): 

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

654 else: 

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

656 # we define the following default value: 

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

658 

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

660 # same rotation. 

661 partners, n_partners = _find_close_points( 

662 rotation_vectors, tol=_bme.eps_quaternion 

663 ) 

664 

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

666 if n_partners == 0: 

667 self.add( 

668 _coupling_factory(node_list, coupling_type, coupling_dof_type) 

669 ) 

670 else: 

671 # There are nodes that need to be combined. 

672 combining_nodes = [] 

673 coupling_nodes = [] 

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

675 

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

677 # that will be coupled. 

678 for i, partner in enumerate(partners): 

679 if partner == -1: 

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

681 # rotation. 

682 coupling_nodes.append(node_list[i]) 

683 

684 elif found_partner_id[partner] is not None: 

685 # This node has already a processed partner, add 

686 # this one to the combining nodes. 

687 combining_nodes[found_partner_id[partner]].append( 

688 node_list[i] 

689 ) 

690 

691 else: 

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

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

694 # be replaced with this one. 

695 new_index = len(combining_nodes) 

696 found_partner_id[partner] = new_index 

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

698 coupling_nodes.append(node_list[i]) 

699 

700 # Add the coupling nodes. 

701 if len(coupling_nodes) > 1: 

702 self.add( 

703 _coupling_factory( 

704 coupling_nodes, coupling_type, coupling_dof_type 

705 ) 

706 ) 

707 

708 # Replace the identical nodes. 

709 for combine_list in combining_nodes: 

710 master_node = combine_list[0] 

711 for node in combine_list[1:]: 

712 node.replace_with(master_node) 

713 

714 else: 

715 # Connect close nodes with a coupling. 

716 for node_list in partner_nodes: 

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

718 

719 def unlink_nodes(self): 

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

721 for node in self.nodes: 

722 node.unlink() 

723 

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

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

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

727 

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

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

730 directions.""" 

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

732 

733 def check_overlapping_elements(self, raise_error=True): 

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

735 

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

737 have unique coordinates in the mesh. 

738 """ 

739 

740 # Number of middle nodes. 

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

742 

743 # Only check if there are middle nodes. 

744 if len(middle_nodes) == 0: 

745 return 

746 

747 # Get array with middle nodes. 

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

749 for i, node in enumerate(middle_nodes): 

750 coordinates[i, :] = node.coordinates 

751 

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

753 has_partner, partner = _find_close_points(coordinates) 

754 partner_indices = _point_partners_to_partner_indices(has_partner, partner) 

755 if partner > 0: 

756 if raise_error: 

757 raise ValueError( 

758 "There are multiple middle nodes with the " 

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

760 "This check can be turned of with " 

761 "bme.check_overlapping_elements=False" 

762 ) 

763 else: 

764 _warnings.warn( 

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

766 ) 

767 

768 # Add the partner index to the middle nodes. 

769 for i_partner, partners in enumerate(partner_indices): 

770 for i_node in partners: 

771 middle_nodes[i_node].element_partner_index = i_partner 

772 

773 def get_vtk_representation( 

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

775 ): 

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

777 

778 Args 

779 ---- 

780 overlapping_elements: bool 

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

782 output will mark them. 

783 coupling_sets: bool 

784 If coupling sets should also be displayed. 

785 """ 

786 

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

788 vtk_writer_beam = _VTKWriter() 

789 vtk_writer_solid = _VTKWriter() 

790 

791 # Get the set numbers of the mesh 

792 mesh_sets = self.get_unique_geometry_sets( 

793 coupling_sets=coupling_sets, link_to_nodes="all_nodes" 

794 ) 

795 

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

797 # Get highest number of node_sets. 

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

799 

800 # Set the bme value. 

801 digits = len(str(max_sets)) 

802 _bme.vtk_node_set_format = "{:0" + str(digits) + "}" 

803 

804 if overlapping_elements: 

805 # Check for overlapping elements. 

806 self.check_overlapping_elements(raise_error=False) 

807 

808 # Get representation of elements. 

809 for element in self.elements: 

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

811 

812 # Finish and return the writers 

813 vtk_writer_beam.complete_data() 

814 vtk_writer_solid.complete_data() 

815 return vtk_writer_beam, vtk_writer_solid 

816 

817 def write_vtk( 

818 self, output_name="beamme", output_directory="", binary=True, **kwargs 

819 ): 

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

821 

822 Args 

823 ---- 

824 output_name: str 

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

826 {}_solid.vtu file. 

827 output_directory: path 

828 Directory where the output files will be written. 

829 binary: bool 

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

831 

832 **kwargs 

833 For all of them look into: 

834 Mesh().get_vtk_representation 

835 Beam().get_vtk 

836 VolumeElement().get_vtk 

837 ---- 

838 beam_centerline_visualization_segments: int 

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

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

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

842 visualization purposes. 

843 """ 

844 

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

846 

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

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

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

850 vtk_writer_beam.write_vtk(filepath, binary=binary) 

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

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

853 vtk_writer_solid.write_vtk(filepath, binary=binary) 

854 

855 def display_pyvista( 

856 self, 

857 *, 

858 beam_nodes=True, 

859 beam_tube=True, 

860 beam_cross_section_directors=True, 

861 beam_radius_for_display=None, 

862 resolution=20, 

863 parallel_projection=False, 

864 **kwargs, 

865 ): 

866 """Display the mesh in pyvista. 

867 

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

869 the _pv.plotter object will be returned. 

870 

871 Args 

872 ---- 

873 beam_nodes: bool 

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

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

876 are shown in cyan. 

877 beam_tube: bool 

878 If the beam should be rendered as a tube 

879 beam_cross_section_directors: bool 

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

881 beam_radius_for_display: float 

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

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

884 for visualization 

885 resolution: int 

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

887 tubes and spheres. 

888 parallel_projection: bool 

889 Flag to change camera view to parallel projection. 

890 

891 **kwargs 

892 For all of them look into: 

893 Mesh().get_vtk_representation 

894 Beam().get_vtk 

895 VolumeElement().get_vtk 

896 ---- 

897 beam_centerline_visualization_segments: int 

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

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

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

901 visualization purposes. 

902 """ 

903 

904 plotter = _pv.Plotter() 

905 plotter.renderer.add_axes() 

906 

907 if parallel_projection: 

908 plotter.enable_parallel_projection() 

909 

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

911 

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

913 beam_grid = _pv.UnstructuredGrid(vtk_writer_beam.grid)