Coverage for src/beamme/core/mesh.py: 95%
430 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-08 11:03 +0000
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."""
25import copy as _copy
26import os as _os
27import warnings as _warnings
28from typing import Any as _Any
29from typing import List as _List
30from typing import cast as _cast
32import numpy as _np
33import pyvista as _pv
34import quaternion as _quaternion
35from numpy.typing import NDArray as _NDArray
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 GeometrySet as _GeometrySet
50from beamme.core.geometry_set import GeometrySetBase as _GeometrySetBase
51from beamme.core.geometry_set import GeometrySetContainer as _GeometrySetContainer
52from beamme.core.material import Material as _Material
53from beamme.core.mesh_representation import (
54 MESH_REPRESENTATION_MAPPINGS as _MESH_REPRESENTATION_MAPPINGS,
55)
56from beamme.core.mesh_representation import GeometrySetInfo as _GeometrySetInfo
57from beamme.core.mesh_representation import MeshRepresentation as _MeshRepresentation
58from beamme.core.node import Node as _Node
59from beamme.core.node import NodeCosserat as _NodeCosserat
60from beamme.core.nurbs_patch import NURBSPatch as _NURBSPatch
61from beamme.core.rotation import Rotation as _Rotation
62from beamme.core.rotation import add_rotations as _add_rotations
63from beamme.core.rotation import rotate_coordinates as _rotate_coordinates
64from beamme.core.vtk_writer import VTKWriter as _VTKWriter
65from beamme.geometric_search.find_close_points import (
66 find_close_points as _find_close_points,
67)
68from beamme.utils.environment import is_testing as _is_testing
69from beamme.utils.nodes import filter_nodes as _filter_nodes
70from beamme.utils.nodes import find_close_nodes as _find_close_nodes
71from beamme.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates
72from beamme.utils.nodes import get_nodal_quaternions as _get_nodal_quaternions
73from beamme.utils.nodes import get_nodes_by_function as _get_nodes_by_function
76class Mesh:
77 """A class that contains a full mesh, i.e. Nodes, Elements, Boundary
78 Conditions, Sets, Couplings, Materials and Functions."""
80 def __init__(self):
81 """Initialize all empty containers."""
83 self.nodes = []
84 self.elements = []
85 self.materials = []
86 self.functions = []
87 self.geometry_sets = _GeometrySetContainer()
88 self.boundary_conditions = _BoundaryConditionContainer()
90 @staticmethod
91 def get_base_mesh_item_type(item):
92 """Return the base mesh type of the given item.
94 Amongst other things, we need this function so we can check if
95 all items of a list are of the same "base" type.
97 Args:
98 item: The object we want to get the base type from.
99 """
101 for cls in (
102 Mesh,
103 _Function,
104 _BoundaryConditionBase,
105 _Material,
106 _Node,
107 _Element,
108 _GeometrySetBase,
109 _GeometryName,
110 ):
111 if isinstance(item, cls):
112 return cls
113 return type(item)
115 def add(self, *args, **kwargs):
116 """Add an item to this mesh, depending on its type.
118 If an list is given each list element is added with this
119 function. If multiple arguments are given, each one is
120 individually added with this function. Keyword arguments are
121 passed through to the adding function.
122 """
124 match len(args):
125 case 0:
126 raise ValueError("At least one argument is required!")
127 case 1:
128 add_item = args[0]
129 base_type = self.get_base_mesh_item_type(add_item)
131 base_type_to_method_map = {
132 Mesh: self.add_mesh,
133 _Function: self.add_function,
134 _BoundaryConditionBase: self.add_bc,
135 _Material: self.add_material,
136 _Node: self.add_node,
137 _Element: self.add_element,
138 _GeometrySetBase: self.add_geometry_set,
139 _GeometryName: self.add_geometry_name,
140 list: self.add_list,
141 }
142 if base_type in base_type_to_method_map:
143 base_type_to_method_map[base_type](add_item, **kwargs)
144 else:
145 raise TypeError(
146 f'No Mesh.add case implemented for type: "{type(add_item)}" with base type "{base_type}"!'
147 )
148 case _:
149 for item in args:
150 self.add(item, **kwargs)
152 def add_mesh(self, mesh):
153 """Add the content of another mesh to this mesh."""
155 # Add each item from mesh to self.
156 self.add(mesh.nodes)
157 self.add(mesh.elements)
158 self.add(mesh.materials)
159 self.add(mesh.functions)
160 self.geometry_sets.extend(mesh.geometry_sets)
161 self.boundary_conditions.extend(mesh.boundary_conditions)
163 def add_bc(self, bc):
164 """Add a boundary condition to this mesh."""
165 bc_key = bc.bc_type
166 geom_key = bc.geometry_set.geometry_type
167 bc.geometry_set.check_replaced_nodes()
168 self.boundary_conditions.append((bc_key, geom_key), bc)
170 def add_function(self, function):
171 """Add a function to this mesh item.
173 Check that the function is only added once.
174 """
175 if function not in self.functions:
176 self.functions.append(function)
178 def add_material(self, material):
179 """Add a material to this mesh item.
181 Check that the material is only added once.
182 """
183 if material not in self.materials:
184 self.materials.append(material)
186 def add_node(self, node):
187 """Add a node to this mesh."""
188 if node in self.nodes:
189 raise ValueError("The node is already in this mesh!")
190 self.nodes.append(node)
192 def add_element(self, element):
193 """Add an element to this mesh."""
194 if element in self.elements:
195 raise ValueError("The element is already in this mesh!")
196 self.elements.append(element)
198 def add_geometry_set(self, geometry_set):
199 """Add a geometry set to this mesh."""
200 geometry_set.check_replaced_nodes()
201 self.geometry_sets.append(geometry_set.geometry_type, geometry_set)
203 def add_geometry_name(self, geometry_name):
204 """Add a set of geometry sets to this mesh.
206 Sort by the keys here to create a deterministic ordering,
207 especially for testing purposes
208 """
209 keys = list(geometry_name.keys())
210 keys.sort()
211 for key in keys:
212 self.add(geometry_name[key])
214 def add_list(self, add_list: _List, **kwargs) -> None:
215 """Add a list of items to this mesh.
217 Args:
218 add_list:
219 List to be added to the mesh. This method checks that all
220 base types of the list items are the same.
222 In the special case of a node or element list, we add the whole list
223 at once. This avoids a check for duplicate entries for every
224 addition, which scales very badly. By doing it this way we only
225 check the final list for duplicate entries which is much more
226 performant.
228 For all other types of items, we add each element individually
229 via the Mesh.add method.
230 """
232 types = {self.get_base_mesh_item_type(item) for item in add_list}
233 if len(types) > 1:
234 raise TypeError(
235 f"You can only add lists with the same type of element. Got {types}"
236 )
237 elif len(types) == 1:
238 list_type = types.pop()
240 def extend_internal_list(self_list: _List, new_list: _List) -> None:
241 """Extend an internal list with the new list.
243 It is checked that the final list does not have
244 duplicate entries.
245 """
246 self_list.extend(new_list)
247 if not len(set(self_list)) == len(self_list):
248 raise ValueError(
249 "The added list contains entries already existing in the Mesh"
250 )
252 if list_type == _Node:
253 extend_internal_list(self.nodes, add_list)
254 elif list_type == _Element:
255 extend_internal_list(self.elements, add_list)
256 else:
257 for item in add_list:
258 self.add(item, **kwargs)
260 def replace_nodes(self, replace_nodes: dict[_Node, _Node]) -> None:
261 """Replace all nodes in this mesh that are in the replacement map.
263 Args:
264 replace_nodes: A dictionary that maps source nodes to target nodes. The
265 source nodes will be replaced with the target nodes in the mesh.
266 """
268 # Nothing to do if the replacement map is empty.
269 if len(replace_nodes) == 0:
270 return
272 # Check that all source and target nodes are in the mesh.
273 for items, name in (
274 (replace_nodes.keys(), "source"),
275 (replace_nodes.values(), "target"),
276 ):
277 mesh_contains_all_nodes = set(items).issubset(self.nodes)
278 if not mesh_contains_all_nodes:
279 raise ValueError(f"Not all {name} nodes are in the mesh!")
281 # Remove source nodes from the mesh.
282 self.nodes = [node for node in self.nodes if node not in replace_nodes]
284 # Replace the nodes in the elements.
285 for element in self.elements:
286 for i, node in enumerate(element.nodes):
287 if node in replace_nodes:
288 element.nodes[i] = replace_nodes[node]
290 # Set the links to the target nodes in the source nodes, so they can be
291 # replaced in the geometry sets.
292 for source_node, target_node in replace_nodes.items():
293 source_node.target_node = target_node
295 # Replace the nodes in the geometry sets.
296 mesh_sets = self.get_unique_geometry_sets()
297 for geometry_set_list in mesh_sets.values():
298 for geometry_set in geometry_set_list:
299 geometry_set.check_replaced_nodes()
301 def get_unique_geometry_sets(
302 self, *, coupling_sets: bool = True
303 ) -> _GeometrySetContainer:
304 """Return a geometry set container that contains geometry sets
305 explicitly added to the mesh, as well as sets for boundary conditions.
307 The i_global values are set in the returned geometry sets.
309 Args:
310 coupling_sets:
311 If this is true, also sets for couplings will be added.
313 Returns:
314 A geometry set container that contains all geometry sets of this mesh.
315 """
317 # Make a copy of the sets in this mesh.
318 mesh_sets = self.geometry_sets.copy()
320 # Add sets from boundary conditions.
321 for (bc_key, geom_key), bc_list in self.boundary_conditions.items():
322 for bc in bc_list:
323 # Check if sets from couplings should be added.
324 is_coupling = bc_key in (
325 _bme.bc.point_coupling,
326 bc_key == _bme.bc.point_coupling_penalty,
327 )
328 if (is_coupling and coupling_sets) or (not is_coupling):
329 # Only add set if it is not already in the container.
330 # For example if multiple Neumann boundary conditions
331 # are applied on the same node set.
332 if bc.geometry_set not in mesh_sets[geom_key]:
333 mesh_sets[geom_key].append(bc.geometry_set)
335 return mesh_sets
337 def set_node_links(self):
338 """Create a link of all elements to the nodes connected to them."""
339 for element in self.elements:
340 for node in element.nodes:
341 node.element_link.append(element)
343 def translate(self, vector: _NDArray | list[float]) -> None:
344 """Translate all beam nodes of this mesh.
346 Args:
347 vector: A 3D vector that will be added to all nodes.
348 """
349 for node in self.nodes:
350 node.coordinates += vector
352 def rotate(
353 self,
354 rotation: _Rotation | _NDArray[_quaternion.quaternion],
355 origin=None,
356 only_rotate_triads: bool = False,
357 ) -> None:
358 """Rotate all beam nodes of the mesh with rotation.
360 Args:
361 rotation: The rotation(s) that will be applied to the nodes. If
362 this is an array, it has to hold a quaternion for each node.
363 origin (3D vector): If this is given, the mesh is rotated about
364 this point. Defaults to (0, 0, 0).
365 only_rotate_triads: If this is true, the nodal positions are not
366 changed.
367 """
369 # Get array with all quaternions for the nodes.
370 rot1 = _get_nodal_quaternions(self.nodes)
372 # Apply the rotation to the rotation of all nodes.
373 rot_new = _add_rotations(rotation, rot1)
375 if not only_rotate_triads:
376 # Get array with all positions for the nodes.
377 pos = _get_nodal_coordinates(self.nodes)
378 pos_new = _rotate_coordinates(pos, rotation, origin=origin)
380 for i, node in enumerate(self.nodes):
381 if isinstance(node, _NodeCosserat):
382 node.rotation.q = rot_new[i, :]
383 if not only_rotate_triads:
384 node.coordinates = pos_new[i, :]
386 def reflect(self, normal_vector, origin=None, flip_beams: bool = False) -> None:
387 """Reflect all nodes of the mesh with respect to a plane defined by its
388 normal_vector. Per default the plane goes through the origin, if not a
389 point on the plane can be given with the parameter origin.
391 For the reflection we assume that e1' and e2' are mirrored with respect
392 to the original frame and e3' is in the opposite direction than the
393 mirrored e3.
395 With the defined mirroring strategy, the quaternion to be applied on
396 the existing rotations can be calculated the following way:
397 q[0] = e3 * n
398 q[1,2,3] = e3 x n
399 This constructs a rotation with the rotation axis on the plane, and
400 normal to the vector e3. The rotation angle is twice the angle of e3
401 to n.
403 Args:
404 normal_vector (3D vector): The normal vector of the reflection plane.
405 origin (3D vector): The reflection plane goes through this point.
406 Defaults to (0, 0, 0).
407 flip_beams: When True, the beams are flipped, so that the direction
408 along the beam is reversed.
409 """
411 # Normalize the normal vector.
412 normal_vector = _np.asarray(normal_vector) / _np.linalg.norm(normal_vector)
414 # Get array with all quaternions and positions for the nodes.
415 pos = _get_nodal_coordinates(self.nodes)
416 rot1 = _get_nodal_quaternions(self.nodes)
418 # Check if origin has to be added.
419 if origin is not None:
420 pos -= origin
422 # Get the reflection matrix A.
423 A = _np.eye(3) - 2.0 * _np.outer(normal_vector, normal_vector)
425 # Calculate the new positions.
426 pos_new = _np.dot(pos, A)
428 # Move back from the origin.
429 if origin is not None:
430 pos_new += origin
432 # First get all e3 vectors of the nodes.
433 e3 = _np.zeros_like(pos)
434 e3[:, 0] = 2 * (rot1[:, 0] * rot1[:, 2] + rot1[:, 1] * rot1[:, 3])
435 e3[:, 1] = 2 * (-1 * rot1[:, 0] * rot1[:, 1] + rot1[:, 2] * rot1[:, 3])
436 e3[:, 2] = rot1[:, 0] ** 2 - rot1[:, 1] ** 2 - rot1[:, 2] ** 2 + rot1[:, 3] ** 2
438 # Get the dot and cross product of e3 and the normal vector.
439 rot2 = _np.zeros_like(rot1)
440 rot2[:, 0] = _np.dot(e3, normal_vector)
441 rot2[:, 1:] = _np.cross(e3, normal_vector)
443 # Add to the existing rotations.
444 rot_new = _add_rotations(rot2, rot1)
446 if flip_beams:
447 # To achieve the flip, the triads are rotated with the angle pi
448 # around the e2 axis.
449 rot_flip = _Rotation([0, 1, 0], _np.pi)
450 rot_new = _add_rotations(rot_new, rot_flip)
452 # For solid elements we need to adapt the connectivity to avoid negative Jacobians.
453 # For beam elements this is optional.
454 for element in self.elements:
455 if isinstance(element, _Beam):
456 if flip_beams:
457 element.flip()
458 else:
459 element.flip()
461 # Set the new positions and rotations.
462 for i, node in enumerate(self.nodes):
463 node.coordinates = pos_new[i, :]
464 if isinstance(node, _NodeCosserat):
465 node.rotation.q = rot_new[i, :]
467 def wrap_around_cylinder(
468 self, radius: float | None = None, advanced_warning: bool = True
469 ) -> None:
470 """Wrap the geometry around a cylinder. The y-z plane gets morphed into
471 the z-axis of symmetry. If all nodes are on the same y-z plane, the
472 radius of the created cylinder is the x coordinate of that plane. If
473 the nodes are not on the same y-z plane, the radius has to be given
474 explicitly.
476 Args:
477 radius: If this value is given AND not all nodes are on the same y-z
478 plane, then use this radius for the calculation of phi for all
479 nodes. This might still lead to distorted elements!
480 advanced_warning: If each element should be checked if it is either parallel
481 to the y-z or x-z plane. This is computationally expensive, but in most
482 cases (up to 100,000 elements) this check can be left activated.
483 """
485 pos = _get_nodal_coordinates(self.nodes)
486 quaternions = _np.zeros([len(self.nodes), 4])
488 # The x coordinate is the radius, the y coordinate the arc length.
489 points_x = pos[:, 0].copy()
491 # Check if all points are on the same y-z plane.
492 if _np.abs(_np.min(points_x) - _np.max(points_x)) > _bme.eps_pos:
493 # The points are not all on the y-z plane, get the reference
494 # radius.
495 if radius is not None:
496 if advanced_warning:
497 # Here we check, if each element lays on a plane parallel
498 # to the y-z plane, or parallel to the x-z plane.
499 #
500 # To be exactly sure, we could check the rotations here,
501 # i.e. if they are also in plane.
502 element_warning = []
503 for i_element, element in enumerate(self.elements):
504 element_coordinates = _np.zeros([len(element.nodes), 3])
505 for i_node, node in enumerate(element.nodes):
506 element_coordinates[i_node, :] = node.coordinates
507 is_yz = (
508 _np.max(
509 _np.abs(
510 element_coordinates[:, 0]
511 - element_coordinates[0, 0]
512 )
513 )
514 < _bme.eps_pos
515 )
516 is_xz = (
517 _np.max(
518 _np.abs(
519 element_coordinates[:, 1]
520 - element_coordinates[0, 1]
521 )
522 )
523 < _bme.eps_pos
524 )
525 if not (is_yz or is_xz):
526 element_warning.append(i_element)
527 if len(element_warning) != 0:
528 _warnings.warn(
529 "There are elements which are not "
530 "parallel to the y-z or x-y plane. This will lead "
531 "to distorted elements!"
532 )
533 else:
534 _warnings.warn(
535 "The nodes are not on the same y-z plane. "
536 "This may lead to distorted elements!"
537 )
538 else:
539 raise ValueError(
540 "The nodes that should be wrapped around a "
541 "cylinder are not on the same y-z plane. This will give "
542 "unexpected results. Give a reference radius!"
543 )
544 radius_phi = radius
545 radius_points = points_x
546 elif radius is None or _np.abs(points_x[0] - radius) < _bme.eps_pos:
547 radius_points = radius_phi = points_x[0]
548 else:
549 raise ValueError(
550 (
551 "The points are all on the same y-z plane with "
552 "the x-coordinate {} but the given radius {} is different. "
553 "This does not make sense."
554 ).format(points_x[0], radius)
555 )
557 # Get the angle for all nodes.
558 phi = pos[:, 1] / radius_phi
560 # The rotation is about the z-axis.
561 quaternions[:, 0] = _np.cos(0.5 * phi)
562 quaternions[:, 3] = _np.sin(0.5 * phi)
564 # Set the new positions in the global array.
565 pos[:, 0] = radius_points * _np.cos(phi)
566 pos[:, 1] = radius_points * _np.sin(phi)
568 # Rotate the mesh
569 self.rotate(quaternions, only_rotate_triads=True)
571 # Set the new position for the nodes.
572 for i, node in enumerate(self.nodes):
573 node.coordinates = pos[i, :]
575 def couple_nodes(
576 self,
577 *,
578 nodes=None,
579 reuse_matching_nodes=False,
580 coupling_type=_bme.bc.point_coupling,
581 coupling_dof_type=_bme.coupling_dof.fix,
582 ) -> None:
583 """Search through nodes and connect all nodes with the same
584 coordinates.
586 Args:
587 nodes:
588 List of nodes to couple. If None is given, all nodes of the mesh
589 are coupled (except middle nodes).
590 reuse_matching_nodes:
591 If two nodes have the same position and rotation, the nodes are
592 reduced to one node in the mesh. Be aware, that this might lead to
593 issues if not all DOFs of the nodes should be coupled.
594 coupling_type:
595 Type of point coupling.
596 coupling_dof_type:
597 `str`: The string that will be used in the input file.
598 `bme.coupling_dof.fix`: Fix all positional and rotational DOFs of the
599 nodes together.
600 `bme.coupling_dof.joint`: Fix all positional DOFs of the nodes
601 together.
602 """
604 # Check that a coupling BC is given.
605 if coupling_type not in (
606 _bme.bc.point_coupling,
607 _bme.bc.point_coupling_penalty,
608 ):
609 raise ValueError(
610 "Only coupling conditions can be applied in 'couple_nodes'!"
611 )
613 # Get the nodes that should be checked for coupling. Middle nodes are
614 # not checked, as coupling can only be applied to the boundary nodes.
615 if nodes is None:
616 node_list = self.nodes
617 else:
618 node_list = nodes
619 node_list = _filter_nodes(node_list, middle_nodes=False)
620 partner_nodes = _find_close_nodes(node_list)
621 if len(partner_nodes) == 0:
622 # If no partner nodes were found, end this function.
623 return
625 if reuse_matching_nodes:
626 # Check if there are nodes with the same rotation. If there are the
627 # nodes are reused, and no coupling is inserted.
629 # Go through partner nodes.
630 node_replacement_map: dict[_Node, _Node] = {}
631 for partner_node_list in partner_nodes:
632 # Get array with rotation vectors.
633 rotation_vectors = _np.zeros([len(partner_node_list), 3])
634 for i, node in enumerate(partner_node_list):
635 if isinstance(node, _NodeCosserat):
636 rotation_vectors[i, :] = node.rotation.get_rotation_vector()
637 else:
638 # For the case of nodes that belong to solid elements,
639 # we define the following default value:
640 rotation_vectors[i, :] = [4 * _np.pi, 0, 0]
642 # Use find close points function to find nodes with the
643 # same rotation.
644 partners, n_partners = _find_close_points(
645 rotation_vectors, tol=_bme.eps_quaternion
646 )
648 # Check if nodes with the same rotations were found.
649 if n_partners == 0:
650 self.add(
651 _coupling_factory(
652 partner_node_list, coupling_type, coupling_dof_type
653 )
654 )
655 else:
656 # There are nodes that need to be combined.
657 combining_nodes: list[list[_Node]] = []
658 coupling_nodes: list[_Node] = []
659 found_partner_id: list[int | None] = [
660 None for _i in range(n_partners)
661 ]
663 # Add the nodes that need to be combined and add the nodes
664 # that will be coupled.
665 for i, partner in enumerate(partners):
666 if partner == -1:
667 # This node does not have a partner with the same
668 # rotation.
669 coupling_nodes.append(partner_node_list[i])
671 elif found_partner_id[partner] is not None:
672 # This node has already a processed partner, add
673 # this one to the combining nodes.
674 combining_nodes[found_partner_id[partner]].append(
675 partner_node_list[i]
676 )
678 else:
679 # This is the first node of a partner set that was
680 # found. This one will remain, the other ones will
681 # be replaced with this one.
682 new_index = len(combining_nodes)
683 found_partner_id[partner] = new_index
684 combining_nodes.append([partner_node_list[i]])
685 coupling_nodes.append(partner_node_list[i])
687 # Add the coupling nodes.
688 if len(coupling_nodes) > 1:
689 self.add(
690 _coupling_factory(
691 coupling_nodes, coupling_type, coupling_dof_type
692 )
693 )
695 # Add to the replacement map.
696 for combine_list in combining_nodes:
697 target_node = combine_list[0]
698 for node in combine_list[1:]:
699 node_replacement_map[node] = target_node
701 # Replace the nodes in the elements and geometry sets.
702 self.replace_nodes(node_replacement_map)
704 else:
705 # Connect close nodes with a coupling.
706 for node_list in partner_nodes:
707 self.add(_coupling_factory(node_list, coupling_type, coupling_dof_type))
709 def unlink_nodes(self):
710 """Delete the linked arrays and global indices in all nodes."""
711 for node in self.nodes:
712 node.unlink()
714 def get_nodes_by_function(self, *args, **kwargs):
715 """Return all nodes for which the function evaluates to true."""
716 return _get_nodes_by_function(self.nodes, *args, **kwargs)
718 def get_mesh_representation(
719 self, material_to_i_global: dict[_Material, int]
720 ) -> tuple[
721 _MeshRepresentation,
722 dict[int, _Any],
723 dict[_GeometrySetBase, int],
724 dict[_NURBSPatch, int],
725 ]:
726 """Create a mesh representation for this mesh.
728 This function does not alter the mesh object. It assigns internal IDs to the
729 mesh object and returns mappings between the objects and the IDs.
731 Args:
732 material_to_i_global: A dictionary that maps materials to their global
733 index in the mesh representation.
735 Returns:
736 mesh_representation: `MeshRepresentation` object for this mesh.
737 element_type_id_to_data: A dictionary that maps the element type id to the
738 data of the element type.
739 geometry_sets_to_i_global: A dictionary that maps geometry sets to their
740 global index in the mesh representation.
741 nurbs_patch_to_i_global: A dictionary that maps each NURBS patch to the
742 global ID of that patch.
743 """
745 # Get the global id mappings for geometry sets.
746 mesh_sets = self.get_unique_geometry_sets()
747 geometry_sets_to_i_global: dict[_GeometrySetBase, int] = {}
748 for geometry_type, geometry_list in mesh_sets.items():
749 for geometry_set in geometry_list:
750 geometry_sets_to_i_global[geometry_set] = len(geometry_sets_to_i_global)
752 # Extract nodes for the mesh representation and assign global IDs. We also extract
753 # other required information like the node type, the nodal rotation vector and
754 # control point weights information here. The weights are optional, so the array
755 # might be `None`.
756 n_nodes = len(self.nodes)
757 if n_nodes != len(set(self.nodes)):
758 raise ValueError("Nodes are not unique!")
759 points = _np.zeros((n_nodes, 3))
760 point_types = _np.full(n_nodes, -1)
761 control_point_weights = None
762 for i_node, node in enumerate(self.nodes):
763 node.i_global = i_node
764 node_type = type(node).node_type
765 point_types[i_node] = node_type.value
766 points[i_node] = node.coordinates
767 if node_type == _bme.node_type.control_point:
768 if control_point_weights is None:
769 control_point_weights = _np.full(n_nodes, -1.0)
770 control_point_weights[i_node] = node.weight
771 # We don't get the rotation vectors in the loop, that would be very slow, instead
772 # we get the global quaternion array and convert that directly using the numpy
773 # quaternion library.
774 nodal_quaternions = _quaternion.from_float_array(
775 _get_nodal_quaternions(self.nodes)
776 )
777 nodal_rotation_vectors = _quaternion.as_rotation_vector(nodal_quaternions)
779 # Check that element are unique.
780 if len(self.elements) != len(set(self.elements)):
781 raise ValueError("Elements are not unique!")
783 # For the elements, we first have to loop over all elements, so we get the total
784 # number of elements, as NURBS patches can contain multiple elements.
785 # This is needed to initialize the numpy arrays with the correct size.
786 i_element = 0
787 nurbs_count = 0
788 nurbs_patch_to_i_global = {}
789 for element in self.elements:
790 # Perform consistency checks for the element.
791 element.check()
792 element.i_global = i_element
793 if isinstance(element, _NURBSPatch):
794 nurbs_patch_to_i_global[element] = nurbs_count
795 nurbs_count += 1
796 i_element += element.get_number_of_elements()
797 else:
798 i_element += 1
799 n_elements = i_element
801 # Now that we know the expected size, we can allocate the data arrays and
802 # actually gather the element data.
803 element_type_to_id: dict[type, int] = {}
804 element_type_id_to_data: dict[int, _Any] = {}
805 cell_connectivity = []
806 cell_types = _np.full(n_elements, -1)
807 cell_element_type_ids = _np.full(n_elements, -1)
808 cell_material_ids = _np.full(n_elements, -1)
809 cell_beamme_element_ids = _np.full(n_elements, -1)
810 for i_element_beamme, element in enumerate(self.elements):
811 # Get the element type id for this element.
812 if type(element) not in element_type_to_id:
813 element_type_id = len(element_type_to_id)
814 element_type_to_id[type(element)] = element_type_id
815 element_type_id_to_data[element_type_id] = _copy.deepcopy(
816 type(element).data
817 )
818 else:
819 element_type_id = element_type_to_id[type(element)]
821 # For elements which don't require a material we set the material id to -1.
822 # This is currently only the case for rigid sphere elements in 4C.
823 material_id = material_to_i_global.get(element.material, -1)
825 if isinstance(element, _NURBSPatch):
826 n_patch_elements = element.get_number_of_elements()
827 # To satisfy mypy, we do a cast here, since we know that we have set
828 # i_global previously for all elements.
829 element_i_global = _cast(int, element.i_global)
830 data_assignment_slice = slice(
831 element_i_global, element_i_global + n_patch_elements
832 )
834 for knot_span in element.get_knot_span_iterator():
835 element_cps_ids = element.get_ids_ctrlpts(*knot_span)
836 connectivity = [
837 element.nodes[index].i_global for index in element_cps_ids
838 ]
839 cell_connectivity.extend([len(connectivity), *connectivity])
841 else:
842 data_assignment_slice = element.i_global
844 reorder_indices = _MESH_REPRESENTATION_MAPPINGS[
845 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
846 ].get((type(element).element_type, len(element.nodes)), None)
847 if reorder_indices is not None:
848 connectivity = [
849 element.nodes[index].i_global for index in reorder_indices
850 ]
851 else:
852 connectivity = [node.i_global for node in element.nodes]
854 cell_connectivity.extend([len(connectivity), *connectivity])
856 cell_material_ids[data_assignment_slice] = material_id
857 cell_element_type_ids[data_assignment_slice] = element_type_id
858 cell_types[data_assignment_slice] = type(element).vtk_cell_type
859 cell_beamme_element_ids[data_assignment_slice] = i_element_beamme
861 # Extract geometry sets.
862 geometry_sets = []
863 for geometry_type, geometry_list in mesh_sets.items():
864 for geometry_set in geometry_list:
865 node_set_flag = _np.zeros(n_nodes, dtype=int)
866 node_set_flag[
867 [node.i_global for node in geometry_set.get_all_nodes()]
868 ] = 1
869 if isinstance(geometry_set, _GeometrySet) and (
870 geometry_type == _bme.geo.line
871 or geometry_type == _bme.geo.surface
872 or geometry_type == _bme.geo.volume
873 ):
874 element_set_flag = _np.zeros(n_elements, dtype=int)
875 element_set_indices: list[int] = []
876 for element in geometry_set.get_geometry_objects():
877 if isinstance(element, _NURBSPatch):
878 # For NURBS, we have to set the flag for all elements that are part of the patch.
879 element_i_global = _cast(int, element.i_global)
880 element_set_indices.extend(
881 range(
882 element_i_global,
883 element_i_global + element.get_number_of_elements(),
884 )
885 )
886 else:
887 element_set_indices.append(element.i_global)
888 element_set_flag[element_set_indices] = 1
890 else:
891 element_set_flag = None
892 geometry_set_wrapper = _GeometrySetInfo(
893 geometry_type=geometry_type,
894 i_global=geometry_sets_to_i_global[geometry_set],
895 point_flag_vector=node_set_flag,
896 cell_flag_vector=element_set_flag,
897 )
898 geometry_sets.append(geometry_set_wrapper)
900 # Reset the previously set indices.
901 for i, node in enumerate(self.nodes):
902 node.i_global = None
903 for element in self.elements:
904 element.i_global = None
906 # TODO: Let's sort the cells by their block ID here - in most solvers, the
907 # the elements are internally sorted by their block IDs anyway.
908 mesh_representation = _MeshRepresentation(
909 cell_connectivity=cell_connectivity,
910 cell_types=cell_types,
911 points=points,
912 geometry_sets=geometry_sets,
913 cell_data={
914 "element_type_id": cell_element_type_ids,
915 "material_id": cell_material_ids,
916 "beamme_id": cell_beamme_element_ids,
917 },
918 point_data={
919 "point_type": point_types,
920 "control_point_weight": control_point_weights,
921 "rotation_vector": nodal_rotation_vectors,
922 },
923 )
925 return (
926 mesh_representation,
927 element_type_id_to_data,
928 geometry_sets_to_i_global,
929 nurbs_patch_to_i_global,
930 )
932 def get_vtk_representation(self, *, coupling_sets=False, **kwargs):
933 """Return a vtk representation of the beams and solid in this mesh.
935 Args
936 ----
937 coupling_sets: bool
938 If coupling sets should also be displayed.
939 """
941 # Object to store VKT data (and write it to file)
942 vtk_writer_beam = _VTKWriter()
943 vtk_writer_solid = _VTKWriter()
945 # Get the set numbers of the mesh
946 mesh_sets = self.get_unique_geometry_sets(coupling_sets=coupling_sets)
948 # Set the global value for digits in the VTK output.
949 # Get highest number of node_sets.
950 max_sets = max(len(geometry_list) for geometry_list in mesh_sets.values())
952 # Set the bme value.
953 digits = len(str(max_sets))
954 _bme.vtk_node_set_format = "{:0" + str(digits) + "}"
956 # Get representation of elements.
957 for element in self.elements:
958 element.get_vtk(vtk_writer_beam, vtk_writer_solid, **kwargs)
960 # Finish and return the writers
961 vtk_writer_beam.complete_data()
962 vtk_writer_solid.complete_data()
963 return vtk_writer_beam, vtk_writer_solid
965 def write_vtk(
966 self, output_name="beamme", output_directory="", binary=True, **kwargs
967 ):
968 """Write the contents of this mesh to VTK files.
970 Args
971 ----
972 output_name: str
973 Base name of the output file. There will be a {}_beam.vtu and
974 {}_solid.vtu file.
975 output_directory: path
976 Directory where the output files will be written.
977 binary: bool
978 If the data should be written encoded in binary or in human readable text
980 **kwargs
981 For all of them look into:
982 Mesh().get_vtk_representation
983 Beam().get_vtk
984 VolumeElement().get_vtk
985 ----
986 beam_centerline_visualization_segments: int
987 Number of segments to be used for visualization of beam centerline between successive
988 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For
989 Values greater than 1, a Hermite interpolation of the centerline is assumed for
990 visualization purposes.
991 """
993 vtk_writer_beam, vtk_writer_solid = self.get_vtk_representation(**kwargs)
995 # Write to file, only if there is at least one point in the writer.
996 if vtk_writer_beam.points.GetNumberOfPoints() > 0:
997 filepath = _os.path.join(output_directory, output_name + "_beam.vtu")
998 vtk_writer_beam.write_vtk(filepath, binary=binary)
999 if vtk_writer_solid.points.GetNumberOfPoints() > 0:
1000 filepath = _os.path.join(output_directory, output_name + "_solid.vtu")
1001 vtk_writer_solid.write_vtk(filepath, binary=binary)
1003 def display_pyvista(
1004 self,
1005 *,
1006 beam_nodes=True,
1007 beam_tube=True,
1008 beam_cross_section_directors=True,
1009 beam_radius_for_display=None,
1010 resolution=20,
1011 parallel_projection=False,
1012 **kwargs,
1013 ):
1014 """Display the mesh in pyvista.
1016 If this is called in a GitHub testing run, nothing will be shown, instead
1017 the _pv.plotter object will be returned.
1019 Args
1020 ----
1021 beam_nodes: bool
1022 If the beam nodes should be displayed. The start and end nodes of each
1023 beam will be shown in green, possible middle nodes inside the element
1024 are shown in cyan.
1025 beam_tube: bool
1026 If the beam should be rendered as a tube
1027 beam_cross_section_directors: bool
1028 If the cross section directors should be displayed (at each node)
1029 beam_radius_for_display: float
1030 If not all beams have an explicitly given radius (in the material
1031 definition) this value will be used to approximate the beams radius
1032 for visualization
1033 resolution: int
1034 Indicates how many triangulations will be performed to visualize arrows,
1035 tubes and spheres.
1036 parallel_projection: bool
1037 Flag to change camera view to parallel projection.
1039 **kwargs
1040 For all of them look into:
1041 Mesh().get_vtk_representation
1042 Beam().get_vtk
1043 VolumeElement().get_vtk
1044 ----
1045 beam_centerline_visualization_segments: int
1046 Number of segments to be used for visualization of beam centerline between successive
1047 nodes. Default is 1, which means a straight line is drawn between the beam nodes. For
1048 Values greater than 1, a Hermite interpolation of the centerline is assumed for
1049 visualization purposes.
1050 """
1052 plotter = _pv.Plotter()
1053 plotter.renderer.add_axes()
1055 if parallel_projection:
1056 plotter.enable_parallel_projection()
1058 vtk_writer_beam, vtk_writer_solid = self.get_vtk_representation(**kwargs)
1060 if vtk_writer_beam.points.GetNumberOfPoints() > 0:
1061 beam_grid = _pv.UnstructuredGrid(vtk_writer_beam.grid)
1063 # Check if all beams have a given cross-section radius, if not set the given input
1064 # value
1065 all_beams_have_cross_section_radius = (
1066 min(beam_grid.cell_data["cross_section_radius"]) > 0
1067 )
1068 if not all_beams_have_cross_section_radius:
1069 if beam_radius_for_display is None:
1070 raise ValueError(
1071 "Not all beams have a radius, you need to set "
1072 "beam_radius_for_display to allow a display of the beams."
1073 )
1074 beam_grid.cell_data["cross_section_radius"] = beam_radius_for_display
1076 # Grid with beam polyline
1077 beam_grid = beam_grid.cell_data_to_point_data()
1079 # Poly data for nodes
1080 finite_element_nodes = beam_grid.cast_to_poly_points().threshold(
1081 scalars="node_value", value=(0.4, 1.1)
1082 )
1084 # Plot the nodes
1085 node_radius_scaling_factor = 1.5
1086 if beam_nodes:
1087 sphere = _pv.Sphere(
1088 radius=1.0,
1089 theta_resolution=resolution,
1090 phi_resolution=resolution,
1091 )
1092 nodes_glyph = finite_element_nodes.glyph(
1093 geom=sphere,
1094 scale="cross_section_radius",
1095 factor=node_radius_scaling_factor,
1096 orient=False,
1097 )
1098 plotter.add_mesh(
1099 nodes_glyph.threshold(scalars="node_value", value=(0.9, 1.1)),
1100 color="green",
1101 )
1102 middle_nodes = nodes_glyph.threshold(
1103 scalars="node_value", value=(0.4, 0.6)
1104 )
1105 if len(middle_nodes.points) > 0:
1106 plotter.add_mesh(middle_nodes, color="cyan")
1108 # Plot the beams
1109 beam_color = [0.5, 0.5, 0.5]
1110 if beam_tube:
1111 surface = beam_grid.extract_surface()
1112 if all_beams_have_cross_section_radius:
1113 tube = surface.tube(
1114 scalars="cross_section_radius",
1115 absolute=True,
1116 n_sides=resolution,
1117 )
1118 else:
1119 tube = surface.tube(
1120 radius=beam_radius_for_display, n_sides=resolution
1121 )
1122 plotter.add_mesh(tube, color=beam_color)
1123 else:
1124 plotter.add_mesh(beam_grid, color=beam_color, line_width=4)
1126 # Plot the directors of the beam cross-section
1127 director_radius_scaling_factor = 3.5
1128 if beam_cross_section_directors:
1129 arrow = _pv.Arrow(
1130 tip_resolution=resolution, shaft_resolution=resolution
1131 )
1132 directors = [
1133 finite_element_nodes.glyph(
1134 geom=arrow,
1135 orient=f"base_vector_{i + 1}",
1136 scale="cross_section_radius",
1137 factor=director_radius_scaling_factor,
1138 )
1139 for i in range(3)
1140 ]
1141 colors = ["white", "blue", "red"]
1142 for i, arrow in enumerate(directors):
1143 plotter.add_mesh(arrow, color=colors[i])
1145 if vtk_writer_solid.points.GetNumberOfPoints() > 0:
1146 solid_grid = _pv.UnstructuredGrid(vtk_writer_solid.grid).clean()
1147 plotter.add_mesh(solid_grid, color="white", show_edges=True, opacity=0.5)
1149 if not _is_testing():
1150 plotter.show()
1151 else:
1152 return plotter
1154 def copy(self) -> "Mesh":
1155 """Return a deep copy of this mesh.
1157 The internal mesh data (nodes, elements, boundary conditions, and
1158 geometry sets) are deep-copied. Materials and functions are not
1159 deep-copied.
1161 **Important:** Some mesh creation functions return geometry set
1162 containers (e.g., node or element sets) that hold a reference to the
1163 nodes or elements of the mesh they were created with. When using
1164 ``mesh.copy()``, these externally returned sets remain linked to the
1165 original mesh and are therefore not transferred to the copied mesh.
1167 To copy both the mesh and the corresponding geometry sets correctly,
1168 deep-copy them together.
1170 Returns:
1171 A deep copy of the mesh.
1172 """
1173 return _copy.deepcopy(self)