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

1# The MIT License (MIT) 

2# 

3# Copyright (c) 2018-2026 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_nodal_coordinates as _get_nodal_coordinates 

68from beamme.utils.nodes import get_nodal_quaternions as _get_nodal_quaternions 

69from beamme.utils.nodes import get_nodes_by_function as _get_nodes_by_function 

70 

71 

72class Mesh: 

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

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

75 

76 def __init__(self): 

77 """Initialize all empty containers.""" 

78 

79 self.nodes = [] 

80 self.elements = [] 

81 self.materials = [] 

82 self.functions = [] 

83 self.geometry_sets = _GeometrySetContainer() 

84 self.boundary_conditions = _BoundaryConditionContainer() 

85 

86 @staticmethod 

87 def get_base_mesh_item_type(item): 

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

89 

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

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

92 

93 Args: 

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

95 """ 

96 

97 for cls in ( 

98 Mesh, 

99 _Function, 

100 _BoundaryConditionBase, 

101 _Material, 

102 _Node, 

103 _Element, 

104 _GeometrySetBase, 

105 _GeometryName, 

106 ): 

107 if isinstance(item, cls): 

108 return cls 

109 return type(item) 

110 

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

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

113 

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

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

116 individually added with this function. Keyword arguments are 

117 passed through to the adding function. 

118 """ 

119 

120 match len(args): 

121 case 0: 

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

123 case 1: 

124 add_item = args[0] 

125 base_type = self.get_base_mesh_item_type(add_item) 

126 

127 base_type_to_method_map = { 

128 Mesh: self.add_mesh, 

129 _Function: self.add_function, 

130 _BoundaryConditionBase: self.add_bc, 

131 _Material: self.add_material, 

132 _Node: self.add_node, 

133 _Element: self.add_element, 

134 _GeometrySetBase: self.add_geometry_set, 

135 _GeometryName: self.add_geometry_name, 

136 list: self.add_list, 

137 } 

138 if base_type in base_type_to_method_map: 

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

140 else: 

141 raise TypeError( 

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

143 ) 

144 case _: 

145 for item in args: 

146 self.add(item, **kwargs) 

147 

148 def add_mesh(self, mesh): 

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

150 

151 # Add each item from mesh to self. 

152 self.add(mesh.nodes) 

153 self.add(mesh.elements) 

154 self.add(mesh.materials) 

155 self.add(mesh.functions) 

156 self.geometry_sets.extend(mesh.geometry_sets) 

157 self.boundary_conditions.extend(mesh.boundary_conditions) 

158 

159 def add_bc(self, bc): 

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

161 bc_key = bc.bc_type 

162 geom_key = bc.geometry_set.geometry_type 

163 bc.geometry_set.check_replaced_nodes() 

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

165 

166 def add_function(self, function): 

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

168 

169 Check that the function is only added once. 

170 """ 

171 if function not in self.functions: 

172 self.functions.append(function) 

173 

174 def add_material(self, material): 

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

176 

177 Check that the material is only added once. 

178 """ 

179 if material not in self.materials: 

180 self.materials.append(material) 

181 

182 def add_node(self, node): 

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

184 if node in self.nodes: 

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

186 self.nodes.append(node) 

187 

188 def add_element(self, element): 

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

190 if element in self.elements: 

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

192 self.elements.append(element) 

193 

194 def add_geometry_set(self, geometry_set): 

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

196 geometry_set.check_replaced_nodes() 

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

198 

199 def add_geometry_name(self, geometry_name): 

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

201 

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

203 especially for testing purposes 

204 """ 

205 keys = list(geometry_name.keys()) 

206 keys.sort() 

207 for key in keys: 

208 self.add(geometry_name[key]) 

209 

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

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

212 

213 Args: 

214 add_list: 

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

216 base types of the list items are the same. 

217 

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

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

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

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

222 performant. 

223 

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

225 via the Mesh.add method. 

226 """ 

227 

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

229 if len(types) > 1: 

230 raise TypeError( 

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

232 ) 

233 elif len(types) == 1: 

234 list_type = types.pop() 

235 

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

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

238 

239 It is checked that the final list does not have 

240 duplicate entries. 

241 """ 

242 self_list.extend(new_list) 

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

244 raise ValueError( 

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

246 ) 

247 

248 if list_type == _Node: 

249 extend_internal_list(self.nodes, add_list) 

250 elif list_type == _Element: 

251 extend_internal_list(self.elements, add_list) 

252 else: 

253 for item in add_list: 

254 self.add(item, **kwargs) 

255 

256 def replace_node(self, old_node, new_node): 

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

258 

259 # Check that the new node is in mesh. 

260 if new_node not in self.nodes: 

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

262 

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

264 if node == old_node: 

265 del self.nodes[i] 

266 break 

267 else: 

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

269 

270 def get_unique_geometry_sets( 

271 self, 

272 *, 

273 coupling_sets: bool = True, 

274 link_to_nodes: str = "no_link", 

275 geometry_set_start_indices: _Optional[_Dict] = None, 

276 ) -> _GeometrySetContainer: 

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

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

279 

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

281 

282 Args: 

283 coupling_sets: 

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

285 link_to_nodes: 

286 "no_link": 

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

288 "explicitly_contained_nodes": 

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

290 "all_nodes": 

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

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

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

294 geometry_set_start_indices: 

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

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

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

298 """ 

299 

300 is_link_nodes = not link_to_nodes == "no_link" 

301 if is_link_nodes: 

302 # First clear all links in existing nodes. 

303 for node in self.nodes: 

304 node.node_sets_link = [] 

305 

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

307 mesh_sets = self.geometry_sets.copy() 

308 

309 # Add sets from boundary conditions. 

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

311 for bc in bc_list: 

312 # Check if sets from couplings should be added. 

313 is_coupling = bc_key in ( 

314 _bme.bc.point_coupling, 

315 bc_key == _bme.bc.point_coupling_penalty, 

316 ) 

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

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

319 # For example if multiple Neumann boundary conditions 

320 # are applied on the same node set. 

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

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

323 

324 for key in mesh_sets.keys(): 

325 i_global_offset = 0 

326 if geometry_set_start_indices is not None: 

327 if key in geometry_set_start_indices: 

328 i_global_offset = geometry_set_start_indices[key] 

329 else: 

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

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

332 # Add global indices to the geometry set. 

333 geometry_set.i_global = i + i_global_offset 

334 if is_link_nodes: 

335 geometry_set.link_to_nodes(link_to_nodes=link_to_nodes) 

336 

337 return mesh_sets 

338 

339 def set_node_links(self): 

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

341 

342 Also add a link to this mesh. 

343 """ 

344 for element in self.elements: 

345 for node in element.nodes: 

346 node.element_link.append(element) 

347 for node in self.nodes: 

348 node.mesh = self 

349 

350 def translate(self, vector): 

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

352 

353 Args 

354 ---- 

355 vector: _np.array, list 

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

357 """ 

358 for node in self.nodes: 

359 node.coordinates += vector 

360 

361 def rotate( 

362 self, 

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

364 origin=None, 

365 only_rotate_triads: bool = False, 

366 ) -> None: 

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

368 

369 Args: 

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

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

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

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

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

375 changed. 

376 """ 

377 

378 # Get array with all quaternions for the nodes. 

379 rot1 = _get_nodal_quaternions(self.nodes) 

380 

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

382 rot_new = _add_rotations(rotation, rot1) 

383 

384 if not only_rotate_triads: 

385 # Get array with all positions for the nodes. 

386 pos = _get_nodal_coordinates(self.nodes) 

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

388 

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

390 if isinstance(node, _NodeCosserat): 

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

392 if not only_rotate_triads: 

393 node.coordinates = pos_new[i, :] 

394 

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

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

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

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

399 

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

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

402 mirrored e3. 

403 

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

405 the existing rotations can be calculated the following way: 

406 q[0] = e3 * n 

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

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

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

410 to n. 

411 

412 Args: 

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

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

415 Defaults to (0, 0, 0). 

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

417 along the beam is reversed. 

418 """ 

419 

420 # Normalize the normal vector. 

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

422 

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

424 pos = _get_nodal_coordinates(self.nodes) 

425 rot1 = _get_nodal_quaternions(self.nodes) 

426 

427 # Check if origin has to be added. 

428 if origin is not None: 

429 pos -= origin 

430 

431 # Get the reflection matrix A. 

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

433 

434 # Calculate the new positions. 

435 pos_new = _np.dot(pos, A) 

436 

437 # Move back from the origin. 

438 if origin is not None: 

439 pos_new += origin 

440 

441 # First get all e3 vectors of the nodes. 

442 e3 = _np.zeros_like(pos) 

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

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

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

446 

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

448 rot2 = _np.zeros_like(rot1) 

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

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

451 

452 # Add to the existing rotations. 

453 rot_new = _add_rotations(rot2, rot1) 

454 

455 if flip_beams: 

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

457 # around the e2 axis. 

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

459 rot_new = _add_rotations(rot_new, rot_flip) 

460 

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

462 # For beam elements this is optional. 

463 for element in self.elements: 

464 if isinstance(element, _Beam): 

465 if flip_beams: 

466 element.flip() 

467 else: 

468 element.flip() 

469 

470 # Set the new positions and rotations. 

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

472 node.coordinates = pos_new[i, :] 

473 if isinstance(node, _NodeCosserat): 

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

475 

476 def wrap_around_cylinder( 

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

478 ) -> None: 

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

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

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

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

483 explicitly. 

484 

485 Args: 

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

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

488 nodes. This might still lead to distorted elements! 

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

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

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

492 """ 

493 

494 pos = _get_nodal_coordinates(self.nodes) 

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

496 

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

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

499 

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

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

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

503 # radius. 

504 if radius is not None: 

505 if advanced_warning: 

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

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

508 # 

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

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

511 element_warning = [] 

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

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

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

515 element_coordinates[i_node, :] = node.coordinates 

516 is_yz = ( 

517 _np.max( 

518 _np.abs( 

519 element_coordinates[:, 0] 

520 - element_coordinates[0, 0] 

521 ) 

522 ) 

523 < _bme.eps_pos 

524 ) 

525 is_xz = ( 

526 _np.max( 

527 _np.abs( 

528 element_coordinates[:, 1] 

529 - element_coordinates[0, 1] 

530 ) 

531 ) 

532 < _bme.eps_pos 

533 ) 

534 if not (is_yz or is_xz): 

535 element_warning.append(i_element) 

536 if len(element_warning) != 0: 

537 _warnings.warn( 

538 "There are elements which are not " 

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

540 "to distorted elements!" 

541 ) 

542 else: 

543 _warnings.warn( 

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

545 "This may lead to distorted elements!" 

546 ) 

547 else: 

548 raise ValueError( 

549 "The nodes that should be wrapped around a " 

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

551 "unexpected results. Give a reference radius!" 

552 ) 

553 radius_phi = radius 

554 radius_points = points_x 

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

556 radius_points = radius_phi = points_x[0] 

557 else: 

558 raise ValueError( 

559 ( 

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

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

562 "This does not make sense." 

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

564 ) 

565 

566 # Get the angle for all nodes. 

567 phi = pos[:, 1] / radius_phi 

568 

569 # The rotation is about the z-axis. 

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

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

572 

573 # Set the new positions in the global array. 

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

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

576 

577 # Rotate the mesh 

578 self.rotate(quaternions, only_rotate_triads=True) 

579 

580 # Set the new position for the nodes. 

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

582 node.coordinates = pos[i, :] 

583 

584 def couple_nodes( 

585 self, 

586 *, 

587 nodes=None, 

588 reuse_matching_nodes=False, 

589 coupling_type=_bme.bc.point_coupling, 

590 coupling_dof_type=_bme.coupling_dof.fix, 

591 ): 

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

593 coordinates. 

594 

595 Args: 

596 ---- 

597 nodes: [Node] 

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

599 are coupled (except middle nodes). 

600 reuse_matching_nodes: bool 

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

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

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

604 coupling_type: bme.bc 

605 Type of point coupling. 

606 coupling_dof_type: str, bme.coupling_dof 

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

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

609 nodes together. 

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

611 together. 

612 """ 

613 

614 # Check that a coupling BC is given. 

615 if coupling_type not in ( 

616 _bme.bc.point_coupling, 

617 _bme.bc.point_coupling_penalty, 

618 ): 

619 raise ValueError( 

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

621 ) 

622 

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

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

625 if nodes is None: 

626 node_list = self.nodes 

627 else: 

628 node_list = nodes 

629 node_list = _filter_nodes(node_list, middle_nodes=False) 

630 partner_nodes = _find_close_nodes(node_list) 

631 if len(partner_nodes) == 0: 

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

633 return 

634 

635 if reuse_matching_nodes: 

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

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

638 

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

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

641 # GeometrySetNodes. 

642 self.unlink_nodes() 

643 self.get_unique_geometry_sets(link_to_nodes="explicitly_contained_nodes") 

644 self.set_node_links() 

645 

646 # Go through partner nodes. 

647 for node_list in partner_nodes: 

648 # Get array with rotation vectors. 

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

650 for i, node in enumerate(node_list): 

651 if isinstance(node, _NodeCosserat): 

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

653 else: 

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

655 # we define the following default value: 

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

657 

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

659 # same rotation. 

660 partners, n_partners = _find_close_points( 

661 rotation_vectors, tol=_bme.eps_quaternion 

662 ) 

663 

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

665 if n_partners == 0: 

666 self.add( 

667 _coupling_factory(node_list, coupling_type, coupling_dof_type) 

668 ) 

669 else: 

670 # There are nodes that need to be combined. 

671 combining_nodes = [] 

672 coupling_nodes = [] 

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

674 

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

676 # that will be coupled. 

677 for i, partner in enumerate(partners): 

678 if partner == -1: 

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

680 # rotation. 

681 coupling_nodes.append(node_list[i]) 

682 

683 elif found_partner_id[partner] is not None: 

684 # This node has already a processed partner, add 

685 # this one to the combining nodes. 

686 combining_nodes[found_partner_id[partner]].append( 

687 node_list[i] 

688 ) 

689 

690 else: 

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

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

693 # be replaced with this one. 

694 new_index = len(combining_nodes) 

695 found_partner_id[partner] = new_index 

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

697 coupling_nodes.append(node_list[i]) 

698 

699 # Add the coupling nodes. 

700 if len(coupling_nodes) > 1: 

701 self.add( 

702 _coupling_factory( 

703 coupling_nodes, coupling_type, coupling_dof_type 

704 ) 

705 ) 

706 

707 # Replace the identical nodes. 

708 for combine_list in combining_nodes: 

709 master_node = combine_list[0] 

710 for node in combine_list[1:]: 

711 node.replace_with(master_node) 

712 

713 else: 

714 # Connect close nodes with a coupling. 

715 for node_list in partner_nodes: 

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

717 

718 def unlink_nodes(self): 

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

720 for node in self.nodes: 

721 node.unlink() 

722 

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

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

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

726 

727 def check_overlapping_elements(self, raise_error=True): 

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

729 

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

731 have unique coordinates in the mesh. 

732 """ 

733 

734 # Number of middle nodes. 

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

736 

737 # Only check if there are middle nodes. 

738 if len(middle_nodes) == 0: 

739 return 

740 

741 # Get array with middle nodes. 

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

743 for i, node in enumerate(middle_nodes): 

744 coordinates[i, :] = node.coordinates 

745 

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

747 has_partner, partner = _find_close_points(coordinates) 

748 partner_indices = _point_partners_to_partner_indices(has_partner, partner) 

749 if partner > 0: 

750 if raise_error: 

751 raise ValueError( 

752 "There are multiple middle nodes with the " 

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

754 "This check can be turned of with " 

755 "bme.check_overlapping_elements=False" 

756 ) 

757 else: 

758 _warnings.warn( 

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

760 ) 

761 

762 # Add the partner index to the middle nodes. 

763 for i_partner, partners in enumerate(partner_indices): 

764 for i_node in partners: 

765 middle_nodes[i_node].element_partner_index = i_partner 

766 

767 def get_vtk_representation( 

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

769 ): 

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

771 

772 Args 

773 ---- 

774 overlapping_elements: bool 

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

776 output will mark them. 

777 coupling_sets: bool 

778 If coupling sets should also be displayed. 

779 """ 

780 

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

782 vtk_writer_beam = _VTKWriter() 

783 vtk_writer_solid = _VTKWriter() 

784 

785 # Get the set numbers of the mesh 

786 mesh_sets = self.get_unique_geometry_sets( 

787 coupling_sets=coupling_sets, link_to_nodes="all_nodes" 

788 ) 

789 

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

791 # Get highest number of node_sets. 

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

793 

794 # Set the bme value. 

795 digits = len(str(max_sets)) 

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

797 

798 if overlapping_elements: 

799 # Check for overlapping elements. 

800 self.check_overlapping_elements(raise_error=False) 

801 

802 # Get representation of elements. 

803 for element in self.elements: 

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

805 

806 # Finish and return the writers 

807 vtk_writer_beam.complete_data() 

808 vtk_writer_solid.complete_data() 

809 return vtk_writer_beam, vtk_writer_solid 

810 

811 def write_vtk( 

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

813 ): 

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

815 

816 Args 

817 ---- 

818 output_name: str 

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

820 {}_solid.vtu file. 

821 output_directory: path 

822 Directory where the output files will be written. 

823 binary: bool 

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

825 

826 **kwargs 

827 For all of them look into: 

828 Mesh().get_vtk_representation 

829 Beam().get_vtk 

830 VolumeElement().get_vtk 

831 ---- 

832 beam_centerline_visualization_segments: int 

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

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

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

836 visualization purposes. 

837 """ 

838 

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

840 

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

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

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

844 vtk_writer_beam.write_vtk(filepath, binary=binary) 

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

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

847 vtk_writer_solid.write_vtk(filepath, binary=binary) 

848 

849 def display_pyvista( 

850 self, 

851 *, 

852 beam_nodes=True, 

853 beam_tube=True, 

854 beam_cross_section_directors=True, 

855 beam_radius_for_display=None, 

856 resolution=20, 

857 parallel_projection=False, 

858 **kwargs, 

859 ): 

860 """Display the mesh in pyvista. 

861 

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

863 the _pv.plotter object will be returned. 

864 

865 Args 

866 ---- 

867 beam_nodes: bool 

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

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

870 are shown in cyan. 

871 beam_tube: bool 

872 If the beam should be rendered as a tube 

873 beam_cross_section_directors: bool 

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

875 beam_radius_for_display: float 

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

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

878 for visualization 

879 resolution: int 

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

881 tubes and spheres. 

882 parallel_projection: bool 

883 Flag to change camera view to parallel projection. 

884 

885 **kwargs 

886 For all of them look into: 

887 Mesh().get_vtk_representation 

888 Beam().get_vtk 

889 VolumeElement().get_vtk 

890 ---- 

891 beam_centerline_visualization_segments: int 

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

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

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

895 visualization purposes. 

896 """ 

897 

898 plotter = _pv.Plotter() 

899 plotter.renderer.add_axes() 

900 

901 if parallel_projection: 

902 plotter.enable_parallel_projection() 

903 

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

905 

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

907 beam_grid = _pv.UnstructuredGrid(vtk_writer_beam.grid) 

908 

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

910 # value 

911 all_beams_have_cross_section_radius = ( 

912 min(beam_grid.cell_data["cross_section_radius"]) > 0 

913 ) 

914 if not all_beams_have_cross_section_radius: 

915 if beam_radius_for_display is None: