Coverage for src/beamme/core/mesh.py: 96%
463 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-28 15:20 +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, elements, sets,
23...) for a meshed geometry."""
25import copy as _copy
26import warnings as _warnings
27from pathlib import Path as _Path
28from typing import Any as _Any
29from typing import Self as _Self
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.geometric_search.find_close_points import (
65 find_close_points as _find_close_points,
66)
67from beamme.utils.environment import is_testing as _is_testing
68from beamme.utils.nodes import filter_nodes as _filter_nodes
69from beamme.utils.nodes import find_close_nodes as _find_close_nodes
70from beamme.utils.nodes import get_nodal_coordinates as _get_nodal_coordinates
71from beamme.utils.nodes import get_nodal_quaternions as _get_nodal_quaternions
72from beamme.utils.nodes import get_nodes_by_function as _get_nodes_by_function
73from beamme.utils.visualization import show_plotter as _show_plotter
76class Mesh:
77 """A class that contains a full mesh, i.e. Nodes, Elements, Boundary Conditions,
78 Sets, Couplings, Materials and Functions."""
80 def __init__(self):
81 """Initialize all empty containers."""
82 self.nodes = []
83 self.elements = []
84 self.materials = []
85 self.functions = []
86 self.geometry_sets = _GeometrySetContainer()
87 self.boundary_conditions = _BoundaryConditionContainer()
89 @staticmethod
90 def get_base_mesh_item_type(item):
91 """Return the base mesh type of the given item.
93 Amongst other things, we need this function so we can check if
94 all items of a list are of the same "base" type.
96 Args:
97 item: The object we want to get the base type from.
98 """
99 for cls in (
100 Mesh,
101 _Function,
102 _BoundaryConditionBase,
103 _Material,
104 _Node,
105 _Element,
106 _GeometrySetBase,
107 _GeometryName,
108 ):
109 if isinstance(item, cls):
110 return cls
111 return type(item)
113 def add(self, *args, **kwargs):
114 """Add an item to this mesh, depending on its type.
116 If an list is given each list element is added with this function. If multiple
117 arguments are given, each one is individually added with this function. Keyword
118 arguments are passed through to the adding function.
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)
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)
148 def add_mesh(self, mesh):
149 """Add the content of another mesh to this mesh."""
150 # Add each item from mesh to self.
151 self.add(mesh.nodes)
152 self.add(mesh.elements)
153 self.add(mesh.materials)
154 self.add(mesh.functions)
155 self.geometry_sets.extend(mesh.geometry_sets)
156 self.boundary_conditions.extend(mesh.boundary_conditions)
158 def add_bc(self, bc):
159 """Add a boundary condition to this mesh."""
160 bc_key = bc.bc_type
161 geom_key = bc.geometry_set.geometry_type
162 bc.geometry_set.check_replaced_nodes()
163 self.boundary_conditions.append((bc_key, geom_key), bc)
165 def add_function(self, function):
166 """Add a function to this mesh item.
168 Check that the function is only added once.
169 """
170 if function not in self.functions:
171 self.functions.append(function)
173 def add_material(self, material):
174 """Add a material to this mesh item.
176 Check that the material is only added once.
177 """
178 if material not in self.materials:
179 self.materials.append(material)
181 def add_node(self, node):
182 """Add a node to this mesh."""
183 if node in self.nodes:
184 raise ValueError("The node is already in this mesh!")
185 self.nodes.append(node)
187 def add_element(self, element):
188 """Add an element to this mesh."""
189 if element in self.elements:
190 raise ValueError("The element is already in this mesh!")
191 self.elements.append(element)
193 def add_geometry_set(self, geometry_set):
194 """Add a geometry set to this mesh."""
195 geometry_set.check_replaced_nodes()
196 self.geometry_sets.append(geometry_set.geometry_type, geometry_set)
198 def add_geometry_name(self, geometry_name):
199 """Add a set of geometry sets to this mesh.
201 Sort by the keys here to create a deterministic ordering, especially for testing
202 purposes
203 """
204 keys = list(geometry_name.keys())
205 keys.sort()
206 for key in keys:
207 self.add(geometry_name[key])
209 def add_list(self, add_list: list, **kwargs) -> None:
210 """Add a list of items to this mesh.
212 Args:
213 add_list:
214 List to be added to the mesh. This method checks that all
215 base types of the list items are the same.
217 In the special case of a node or element list, we add the whole list
218 at once. This avoids a check for duplicate entries for every
219 addition, which scales very badly. By doing it this way we only
220 check the final list for duplicate entries which is much more
221 performant.
223 For all other types of items, we add each element individually
224 via the Mesh.add method.
225 """
226 types = {self.get_base_mesh_item_type(item) for item in add_list}
227 if len(types) > 1:
228 raise TypeError(
229 f"You can only add lists with the same type of element. Got {types}"
230 )
231 elif len(types) == 1:
232 list_type = types.pop()
234 def extend_internal_list(self_list: list, new_list: list) -> None:
235 """Extend an internal list with the new list.
237 It is checked that the final list does not have duplicate entries.
238 """
239 self_list.extend(new_list)
240 if not len(set(self_list)) == len(self_list):
241 raise ValueError(
242 "The added list contains entries already existing in the Mesh"
243 )
245 if list_type == _Node:
246 extend_internal_list(self.nodes, add_list)
247 elif list_type == _Element:
248 extend_internal_list(self.elements, add_list)
249 else:
250 for item in add_list:
251 self.add(item, **kwargs)
253 def replace_nodes(self, replace_nodes: dict[_Node, _Node]) -> None:
254 """Replace all nodes in this mesh that are in the replacement map.
256 Args:
257 replace_nodes: A dictionary that maps source nodes to target nodes. The
258 source nodes will be replaced with the target nodes in the mesh.
259 """
260 # Nothing to do if the replacement map is empty.
261 if len(replace_nodes) == 0:
262 return
264 # Check that all source and target nodes are in the mesh.
265 for items, name in (
266 (replace_nodes.keys(), "source"),
267 (replace_nodes.values(), "target"),
268 ):
269 mesh_contains_all_nodes = set(items).issubset(self.nodes)
270 if not mesh_contains_all_nodes:
271 raise ValueError(f"Not all {name} nodes are in the mesh!")
273 # Remove source nodes from the mesh.
274 self.nodes = [node for node in self.nodes if node not in replace_nodes]
276 # Replace the nodes in the elements.
277 for element in self.elements:
278 for i, node in enumerate(element.nodes):
279 if node in replace_nodes:
280 element.nodes[i] = replace_nodes[node]
282 # Set the links to the target nodes in the source nodes, so they can be
283 # replaced in the geometry sets.
284 for source_node, target_node in replace_nodes.items():
285 source_node.target_node = target_node
287 # Replace the nodes in the geometry sets.
288 mesh_sets = self.get_unique_geometry_sets()
289 for geometry_set_list in mesh_sets.values():
290 for geometry_set in geometry_set_list:
291 geometry_set.check_replaced_nodes()
293 def get_unique_geometry_sets(
294 self, *, coupling_sets: bool = True
295 ) -> _GeometrySetContainer:
296 """Return a geometry set container that contains geometry sets explicitly added
297 to the mesh, as well as sets for boundary conditions.
299 The i_global values are set in the returned geometry sets.
301 Args:
302 coupling_sets:
303 If this is true, also sets for couplings will be added.
305 Returns:
306 A geometry set container that contains all geometry sets of this mesh.
307 """
308 # Make a copy of the sets in this mesh.
309 mesh_sets = self.geometry_sets.copy()
311 # Add sets from boundary conditions.
312 for (bc_key, geom_key), bc_list in self.boundary_conditions.items():
313 for bc in bc_list:
314 # Check if sets from couplings should be added.
315 is_coupling = bc_key in (
316 _bme.bc.point_coupling,
317 bc_key == _bme.bc.point_coupling_penalty,
318 )
319 if (is_coupling and coupling_sets) or (not is_coupling):
320 # Only add set if it is not already in the container.
321 # For example if multiple Neumann boundary conditions
322 # are applied on the same node set.
323 if bc.geometry_set not in mesh_sets[geom_key]:
324 mesh_sets[geom_key].append(bc.geometry_set)
326 return mesh_sets
328 def get_named_geometry_sets(self) -> dict[str, _GeometrySetBase]:
329 """Return a dictionary with the named geometry sets in this mesh.
331 Returns:
332 A dictionary that maps the name of a geometry set to the geometry set
333 object. Only named geometry sets are returned. This function throws an
334 error if there are multiple geometry sets with the same name.
335 """
336 named_geometry_sets = {}
337 for geometry_set_list in self.get_unique_geometry_sets().values():
338 for geometry_set in geometry_set_list:
339 if geometry_set.name is not None:
340 if geometry_set.name in named_geometry_sets:
341 raise ValueError(
342 f"Geometry set name {geometry_set.name} is not unique."
343 )
344 named_geometry_sets[geometry_set.name] = geometry_set
345 return named_geometry_sets
347 def set_node_links(self):
348 """Create a link of all elements to the nodes connected to them."""
349 for element in self.elements:
350 for node in element.nodes:
351 node.element_link.append(element)
353 def translate(self, vector: _NDArray | list[float]) -> None:
354 """Translate all beam nodes of this mesh.
356 Args:
357 vector: A 3D vector that will be added to all nodes.
358 """
359 for node in self.nodes:
360 node.coordinates += vector
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.
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 # Get array with all quaternions for the nodes.
379 rot1 = _get_nodal_quaternions(self.nodes)
381 # Apply the rotation to the rotation of all nodes.
382 rot_new = _add_rotations(rotation, rot1)
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)
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, :]
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.
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.
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.
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 """
420 # Normalize the normal vector.
421 normal_vector = _np.asarray(normal_vector) / _np.linalg.norm(normal_vector)
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)
427 # Check if origin has to be added.
428 if origin is not None:
429 pos -= origin
431 # Get the reflection matrix A.
432 A = _np.eye(3) - 2.0 * _np.outer(normal_vector, normal_vector)
434 # Calculate the new positions.
435 pos_new = _np.dot(pos, A)
437 # Move back from the origin.
438 if origin is not None:
439 pos_new += origin
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
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)
452 # Add to the existing rotations.
453 rot_new = _add_rotations(rot2, rot1)
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)
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()
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, :]
476 def wrap_around_cylinder(
477 self, radius: float | None = None, advanced_warning: bool = True
478 ) -> None:
479 """Wrap the geometry around a cylinder.
481 The y-z plane gets morphed into the z-axis of symmetry. If all nodes are
482 on the same y-z plane, the radius of the created cylinder is the x coordinate
483 of that plane. If the nodes are not on the same y-z plane, the radius has to
484 be given explicitly.
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 pos = _get_nodal_coordinates(self.nodes)
495 quaternions = _np.zeros([len(self.nodes), 4])
497 # The x coordinate is the radius, the y coordinate the arc length.
498 points_x = pos[:, 0].copy()
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 )
566 # Get the angle for all nodes.
567 phi = pos[:, 1] / radius_phi
569 # The rotation is about the z-axis.
570 quaternions[:, 0] = _np.cos(0.5 * phi)
571 quaternions[:, 3] = _np.sin(0.5 * phi)
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)
577 # Rotate the mesh
578 self.rotate(quaternions, only_rotate_triads=True)
580 # Set the new position for the nodes.
581 for i, node in enumerate(self.nodes):
582 node.coordinates = pos[i, :]
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 ) -> None:
592 """Search through nodes and connect all nodes with the same coordinates.
594 Args:
595 nodes:
596 List of nodes to couple. If None is given, all nodes of the mesh
597 are coupled (except middle nodes).
598 reuse_matching_nodes:
599 If two nodes have the same position and rotation, the nodes are
600 reduced to one node in the mesh. Be aware, that this might lead to
601 issues if not all DOFs of the nodes should be coupled.
602 coupling_type:
603 Type of point coupling.
604 coupling_dof_type:
605 `str`: The string that will be used in the input file.
606 `bme.coupling_dof.fix`: Fix all positional and rotational DOFs of the
607 nodes together.
608 `bme.coupling_dof.joint`: Fix all positional DOFs of the nodes
609 together.
610 """
611 # Check that a coupling BC is given.
612 if coupling_type not in (
613 _bme.bc.point_coupling,
614 _bme.bc.point_coupling_penalty,
615 ):
616 raise ValueError(
617 "Only coupling conditions can be applied in 'couple_nodes'!"
618 )
620 # Get the nodes that should be checked for coupling. Middle nodes are
621 # not checked, as coupling can only be applied to the boundary nodes.
622 if nodes is None:
623 node_list = self.nodes
624 else:
625 node_list = nodes
626 node_list = _filter_nodes(node_list, middle_nodes=False)
627 partner_nodes = _find_close_nodes(node_list)
628 if len(partner_nodes) == 0:
629 # If no partner nodes were found, end this function.
630 return
632 if reuse_matching_nodes:
633 # Check if there are nodes with the same rotation. If there are the
634 # nodes are reused, and no coupling is inserted.
636 # Go through partner nodes.
637 node_replacement_map: dict[_Node, _Node] = {}
638 for partner_node_list in partner_nodes:
639 # Get array with rotation vectors.
640 rotation_vectors = _np.zeros([len(partner_node_list), 3])
641 for i, node in enumerate(partner_node_list):
642 if isinstance(node, _NodeCosserat):
643 rotation_vectors[i, :] = node.rotation.get_rotation_vector()
644 else:
645 # For the case of nodes that belong to solid elements,
646 # we define the following default value:
647 rotation_vectors[i, :] = [4 * _np.pi, 0, 0]
649 # Use find close points function to find nodes with the
650 # same rotation.
651 partners, n_partners = _find_close_points(
652 rotation_vectors, tol=_bme.eps_quaternion
653 )
655 # Check if nodes with the same rotations were found.
656 if n_partners == 0:
657 self.add(
658 _coupling_factory(
659 partner_node_list, coupling_type, coupling_dof_type
660 )
661 )
662 else:
663 # There are nodes that need to be combined.
664 combining_nodes: list[list[_Node]] = []
665 coupling_nodes: list[_Node] = []
666 found_partner_id: list[int | None] = [
667 None for _i in range(n_partners)
668 ]
670 # Add the nodes that need to be combined and add the nodes
671 # that will be coupled.
672 for i, partner in enumerate(partners):
673 if partner == -1:
674 # This node does not have a partner with the same
675 # rotation.
676 coupling_nodes.append(partner_node_list[i])
678 elif found_partner_id[partner] is not None:
679 # This node has already a processed partner, add
680 # this one to the combining nodes.
681 combining_nodes[found_partner_id[partner]].append(
682 partner_node_list[i]
683 )
685 else:
686 # This is the first node of a partner set that was
687 # found. This one will remain, the other ones will
688 # be replaced with this one.
689 new_index = len(combining_nodes)
690 found_partner_id[partner] = new_index
691 combining_nodes.append([partner_node_list[i]])
692 coupling_nodes.append(partner_node_list[i])
694 # Add the coupling nodes.
695 if len(coupling_nodes) > 1:
696 self.add(
697 _coupling_factory(
698 coupling_nodes, coupling_type, coupling_dof_type
699 )
700 )
702 # Add to the replacement map.
703 for combine_list in combining_nodes:
704 target_node = combine_list[0]
705 for node in combine_list[1:]:
706 node_replacement_map[node] = target_node
708 # Replace the nodes in the elements and geometry sets.
709 self.replace_nodes(node_replacement_map)
711 else:
712 # Connect close nodes with a coupling.
713 for node_list in partner_nodes:
714 self.add(_coupling_factory(node_list, coupling_type, coupling_dof_type))
716 def unlink_nodes(self):
717 """Delete the linked arrays and global indices in all nodes."""
718 for node in self.nodes:
719 node.unlink()
721 def get_nodes_by_function(self, *args, **kwargs):
722 """Return all nodes for which the function evaluates to true."""
723 return _get_nodes_by_function(self.nodes, *args, **kwargs)
725 def get_mesh_representation(
726 self, material_to_i_global: dict[_Material, int] | None = None
727 ) -> tuple[
728 _MeshRepresentation,
729 dict[int, _Any],
730 dict[_GeometrySetBase, int],
731 dict[_NURBSPatch, int],
732 ]:
733 """Create a mesh representation for this mesh.
735 This function does not alter the mesh object. It assigns internal IDs to the
736 mesh object and returns mappings between the objects and the IDs.
738 Args:
739 material_to_i_global: A dictionary that maps materials to their global
740 index in the mesh representation.
742 Returns:
743 mesh_representation: `MeshRepresentation` object for this mesh.
744 element_type_id_to_data: A dictionary that maps the element type id to the
745 data of the element type.
746 geometry_sets_to_i_global: A dictionary that maps geometry sets to their
747 global index in the mesh representation.
748 nurbs_patch_to_i_global: A dictionary that maps each NURBS patch to the
749 global ID of that patch.
750 """
751 if material_to_i_global is None:
752 material_to_i_global = {}
754 # Get the global id mappings for geometry sets.
755 mesh_sets = self.get_unique_geometry_sets()
756 geometry_sets_to_i_global: dict[_GeometrySetBase, int] = {}
757 for geometry_type, geometry_list in mesh_sets.items():
758 for geometry_set in geometry_list:
759 geometry_sets_to_i_global[geometry_set] = len(geometry_sets_to_i_global)
761 # Extract nodes for the mesh representation and assign global IDs. We also extract
762 # other required information like the node type, the nodal rotation vector and
763 # control point weights information here. The weights are optional, so the array
764 # might be `None`.
765 n_nodes = len(self.nodes)
766 if n_nodes != len(set(self.nodes)):
767 raise ValueError("Nodes are not unique!")
769 points = _np.empty((n_nodes, 3), dtype=float)
770 point_types = _np.empty(n_nodes, dtype=int)
772 # Optional node information, which is only extracted if it is actually needed.
773 control_point_weights = None
774 point_arc_lengths = None
775 point_times = None
777 for i_node, node in enumerate(self.nodes):
778 node.i_global = i_node
779 node_type = type(node).node_type
780 point_types[i_node] = node_type.value
781 points[i_node] = node.coordinates
782 if node_type == _bme.node_type.control_point:
783 if control_point_weights is None:
784 control_point_weights = _np.full(n_nodes, -1.0)
785 control_point_weights[i_node] = node.weight
786 if node_type == _bme.node_type.space_time_cosserat:
787 if point_times is None:
788 point_times = _np.full(n_nodes, _np.nan)
789 point_times[i_node] = node.time
790 if node.arc_length is not None:
791 if point_arc_lengths is None:
792 point_arc_lengths = _np.full(n_nodes, _np.nan)
793 point_arc_lengths[i_node] = node.arc_length
795 # We don't get the rotation vectors in the loop, that would be very slow, instead
796 # we get the global quaternion array and convert that directly using the numpy
797 # quaternion library.
798 nodal_quaternions = _quaternion.from_float_array(
799 _get_nodal_quaternions(self.nodes)
800 )
801 nodal_rotation_vectors = _quaternion.as_rotation_vector(nodal_quaternions)
803 # Check that element are unique.
804 if len(self.elements) != len(set(self.elements)):
805 raise ValueError("Elements are not unique!")
807 # For the elements, we first have to loop over all elements, so we get the total
808 # number of elements, as NURBS patches can contain multiple elements.
809 # This is needed to initialize the numpy arrays with the correct size.
810 i_element = 0
811 nurbs_count = 0
812 nurbs_patch_to_i_global = {}
813 for element in self.elements:
814 # Perform consistency checks for the element.
815 element.check()
816 element.i_global = i_element
817 if isinstance(element, _NURBSPatch):
818 nurbs_patch_to_i_global[element] = nurbs_count
819 nurbs_count += 1
820 i_element += element.get_number_of_elements()
821 else:
822 i_element += 1
823 n_elements = i_element
825 # Now that we know the expected size, we can allocate the data arrays and
826 # actually gather the element data.
827 element_type_to_id: dict[type, int] = {}
828 element_type_id_to_data: dict[int, _Any] = {}
829 cell_connectivity = []
830 cell_types = _np.full(n_elements, -1)
831 cell_element_type_ids = _np.full(n_elements, -1)
832 cell_material_ids = _np.full(n_elements, -1)
833 cell_beamme_element_ids = _np.full(n_elements, -1)
834 for i_element_beamme, element in enumerate(self.elements):
835 # Get the element type id for this element.
836 if type(element) not in element_type_to_id:
837 element_type_id = len(element_type_to_id)
838 element_type_to_id[type(element)] = element_type_id
839 element_type_id_to_data[element_type_id] = _copy.deepcopy(
840 type(element).data
841 )
842 else:
843 element_type_id = element_type_to_id[type(element)]
845 # For elements which don't require a material we set the material id to -1.
846 # This is currently only the case for rigid sphere elements in 4C.
847 material_id = material_to_i_global.get(element.material, -1)
849 if isinstance(element, _NURBSPatch):
850 n_patch_elements = element.get_number_of_elements()
851 # To satisfy mypy, we do a cast here, since we know that we have set
852 # i_global previously for all elements.
853 element_i_global = _cast(int, element.i_global)
854 data_assignment_slice = slice(
855 element_i_global, element_i_global + n_patch_elements
856 )
858 for knot_span in element.get_knot_span_iterator():
859 element_cps_ids = element.get_ids_ctrlpts(*knot_span)
860 connectivity = [
861 element.nodes[index].i_global for index in element_cps_ids
862 ]
863 cell_connectivity.extend([len(connectivity), *connectivity])
865 else:
866 data_assignment_slice = element.i_global
868 reorder_indices = _MESH_REPRESENTATION_MAPPINGS[
869 "element_type_and_n_nodes_to_connectivity_mapping_beamme_to_vtk"
870 ].get((type(element).element_type, len(element.nodes)), None)
871 if reorder_indices is not None:
872 connectivity = [
873 element.nodes[index].i_global for index in reorder_indices
874 ]
875 else:
876 connectivity = [node.i_global for node in element.nodes]
878 cell_connectivity.extend([len(connectivity), *connectivity])
880 cell_material_ids[data_assignment_slice] = material_id
881 cell_element_type_ids[data_assignment_slice] = element_type_id
882 cell_types[data_assignment_slice] = type(element).vtk_cell_type
883 cell_beamme_element_ids[data_assignment_slice] = i_element_beamme
885 # Extract geometry sets.
886 geometry_sets = []
887 for geometry_type, geometry_list in mesh_sets.items():
888 for geometry_set in geometry_list:
889 node_set_flag = _np.zeros(n_nodes, dtype=int)
890 node_set_flag[
891 [node.i_global for node in geometry_set.get_all_nodes()]
892 ] = 1
893 if isinstance(geometry_set, _GeometrySet) and (
894 geometry_type == _bme.geo.line
895 or geometry_type == _bme.geo.surface
896 or geometry_type == _bme.geo.volume
897 ):
898 element_set_flag = _np.zeros(n_elements, dtype=int)
899 element_set_indices: list[int] = []
900 for element in geometry_set.get_geometry_objects():
901 if isinstance(element, _NURBSPatch):
902 # For NURBS, we have to set the flag for all elements that are part of the patch.
903 element_i_global = _cast(int, element.i_global)
904 element_set_indices.extend(
905 range(
906 element_i_global,
907 element_i_global + element.get_number_of_elements(),
908 )
909 )
910 else:
911 element_set_indices.append(_cast(int, element.i_global))
912 element_set_flag[element_set_indices] = 1
914 else:
915 element_set_flag = None
916 geometry_set_wrapper = _GeometrySetInfo(
917 geometry_type=geometry_type,
918 i_global=geometry_sets_to_i_global[geometry_set],
919 point_flag_vector=node_set_flag,
920 cell_flag_vector=element_set_flag,
921 name=geometry_set.name,
922 )
923 geometry_sets.append(geometry_set_wrapper)
925 # Reset the previously set indices.
926 for node in self.nodes:
927 node.i_global = None
928 for element in self.elements:
929 element.i_global = None
931 # Create the mesh representation.
932 mesh_representation = _MeshRepresentation(
933 cell_connectivity=cell_connectivity,
934 cell_types=cell_types,
935 points=points,
936 geometry_sets=geometry_sets,
937 cell_data={
938 "element_type_id": cell_element_type_ids,
939 "material_id": cell_material_ids,
940 "beamme_id": cell_beamme_element_ids,
941 },
942 point_data={
943 "point_type": point_types,
944 "arc_length": point_arc_lengths,
945 "control_point_weight": control_point_weights,
946 "rotation_vector": nodal_rotation_vectors,
947 "time": point_times,
948 },
949 )
951 return (
952 mesh_representation,
953 element_type_id_to_data,
954 geometry_sets_to_i_global,
955 nurbs_patch_to_i_global,
956 )
958 def get_vtu_representation(self) -> _pv.UnstructuredGrid:
959 """Return a vtu representation of this mesh.
961 Returns:
962 A pyvista UnstructuredGrid object that represents this mesh.
963 """
964 # Get mesh representation.
965 mesh_representation, _, _, _ = self.get_mesh_representation()
967 # Get the pyvista grid.
968 grid = mesh_representation.get_pyvista_grid(
969 cell_data_fields=["element_type_id", "beamme_id"],
970 point_data_fields=["point_type", "rotation_vector"],
971 add_geometry_sets=True,
972 )
974 # Get data arrays for visualization, i.e., element type and the cross-section radius for beams.
975 beamme_types = _np.empty(mesh_representation.n_cells, dtype=int)
976 cross_section_radii = _np.full(mesh_representation.n_cells, -1.0)
977 for i_cell, beamme_id in enumerate(
978 mesh_representation.data_iterator("cell_data", "beamme_id")
979 ):
980 element = self.elements[beamme_id]
981 element_type = type(element).element_type
982 beamme_types[i_cell] = element_type.value
983 if element_type == _bme.element_type.beam:
984 cross_section_radii[i_cell] = element.material.radius
985 node_value = _np.zeros(mesh_representation.n_points)
986 for i_node, node in enumerate(self.nodes):
987 if isinstance(node, _NodeCosserat):
988 if node.is_middle_node:
989 node_value[i_node] = 0.5
990 else:
991 node_value[i_node] = 1.0
993 # Add the data arrays created here.
994 grid.cell_data["beamme_type"] = beamme_types
995 grid.cell_data["cross_section_radius"] = cross_section_radii
996 grid.point_data["node_value"] = node_value
998 # Add the triads for each Cosserat node.
999 cosserat_mask = _np.isin(
1000 mesh_representation.point_data["point_type"],
1001 [_bme.node_type.cosserat.value, _bme.node_type.space_time_cosserat.value],
1002 )
1003 quaternions = _quaternion.from_rotation_vector(
1004 grid.point_data["rotation_vector"][cosserat_mask]
1005 )
1006 rotation_matrices = _quaternion.as_rotation_matrix(quaternions)
1007 for i in range(3):
1008 base_vector = _np.zeros((grid.n_points, 3))
1009 base_vector[cosserat_mask] = rotation_matrices[:, :, i]
1010 grid.point_data[f"base_vector_{i + 1}"] = base_vector
1012 return grid
1014 def write_vtu(self, file_name: _Path | str, binary=True):
1015 """Write the contents of this mesh to VTK files.
1017 Args
1018 ----
1019 file_name: The path or filename of the vtu file.
1020 binary: If the data should be written encoded in binary or in human readable text
1021 """
1022 path = _Path(file_name)
1023 if path.suffix == "":
1024 path = path.with_suffix(".vtu")
1025 elif path.suffix != ".vtu":
1026 raise ValueError(f"Expected file extension '.vtu', got '{path.suffix}'")
1027 grid = self.get_vtu_representation()
1028 grid.save(path, binary=binary)
1030 def display_pyvista(
1031 self,
1032 *,
1033 beam_nodes=True,
1034 beam_tube=True,
1035 beam_cross_section_directors=True,
1036 resolution=20,
1037 parallel_projection=False,
1038 ):
1039 """Display the mesh in pyvista.
1041 If this is called in a GitHub testing run, nothing will be shown, instead
1042 the _pv.plotter object will be returned.
1044 Args
1045 ----
1046 beam_nodes: bool
1047 If the beam nodes should be displayed. The start and end nodes of each
1048 beam will be shown in green, possible middle nodes inside the element
1049 are shown in cyan.
1050 beam_tube: bool
1051 If the beam should be rendered as a tube
1052 beam_cross_section_directors: bool
1053 If the cross section directors should be displayed (at each node)
1054 resolution: int
1055 Indicates how many triangulations will be performed to visualize arrows,
1056 tubes and spheres.
1057 parallel_projection: bool
1058 Flag to change camera view to parallel projection.
1059 """
1060 grid = self.get_vtu_representation()
1062 plotter = _pv.Plotter()
1063 plotter.renderer.add_axes()
1065 if parallel_projection:
1066 plotter.enable_parallel_projection()
1068 beam_mask = grid.cell_data["beamme_type"] == _bme.element_type.beam.value
1069 if _np.any(beam_mask):
1070 beam_grid = grid.extract_cells(beam_mask).cell_data_to_point_data()
1072 # Plot the nodes
1073 beam_finite_element_nodes = beam_grid.cast_to_poly_points()
1074 node_radius_scaling_factor = 1.5
1075 if beam_nodes:
1076 sphere = _pv.Sphere(
1077 radius=1.0,
1078 theta_resolution=resolution,
1079 phi_resolution=resolution,
1080 )
1081 nodes_glyph = beam_finite_element_nodes.glyph(
1082 geom=sphere,
1083 scale="cross_section_radius",
1084 factor=node_radius_scaling_factor,
1085 orient=False,
1086 )
1087 plotter.add_mesh(
1088 nodes_glyph.threshold(scalars="node_value", value=(0.9, 1.1)),
1089 color="green",
1090 )
1091 middle_nodes = nodes_glyph.threshold(
1092 scalars="node_value", value=(0.4, 0.6)
1093 )
1094 if len(middle_nodes.points) > 0:
1095 plotter.add_mesh(middle_nodes, color="cyan")
1097 # Plot the beams
1098 beam_color = [0.5, 0.5, 0.5]
1099 if beam_tube:
1100 # Check that all beams have a defined cross-section radius.
1101 if _np.any(beam_grid.point_data["cross_section_radius"] < 0):
1102 raise ValueError(
1103 "All beams must have a defined cross-section radius."
1104 )
1106 # TODO: The `algorithm=None` argument is used to avoid a warning from
1107 # pyvista regarding changed default arguments. For future pyvista
1108 # versions, this argument can be removed.
1109 beam_tube_grid = beam_grid.extract_surface(algorithm=None)
1110 beam_tube_grid = beam_tube_grid.tube(
1111 scalars="cross_section_radius",
1112 absolute=True,
1113 n_sides=resolution,
1114 )
1115 plotter.add_mesh(beam_tube_grid, color=beam_color)
1116 else:
1117 plotter.add_mesh(beam_grid, color=beam_color, line_width=4)
1119 # Plot the directors of the beam cross-section
1120 if beam_cross_section_directors:
1121 director_radius_scaling_factor = 3.5
1122 arrow = _pv.Arrow(
1123 tip_resolution=resolution, shaft_resolution=resolution
1124 )
1125 directors = [
1126 beam_finite_element_nodes.glyph(
1127 geom=arrow,
1128 orient=f"base_vector_{i + 1}",
1129 scale="cross_section_radius",
1130 factor=director_radius_scaling_factor,
1131 )
1132 for i in range(3)
1133 ]
1134 colors = ["white", "blue", "red"]
1135 for i, arrow in enumerate(directors):
1136 plotter.add_mesh(arrow, color=colors[i])
1138 solid_mask = grid.cell_data["beamme_type"] == _bme.element_type.solid.value
1139 if _np.any(solid_mask):
1140 solid_grid = grid.extract_cells(solid_mask)
1141 plotter.add_mesh(solid_grid, color="white", show_edges=True, opacity=0.5)
1143 if not _is_testing():
1144 _show_plotter(plotter)
1145 else:
1146 return plotter
1148 def copy(self) -> _Self:
1149 """Return a deep copy of this mesh.
1151 The internal mesh data (nodes, elements, boundary conditions, and
1152 geometry sets) are deep-copied. Materials and functions are not
1153 deep-copied.
1155 **Important:** Some mesh creation functions return geometry set
1156 containers (e.g., node or element sets) that hold a reference to the
1157 nodes or elements of the mesh they were created with. When using
1158 ``mesh.copy()``, these externally returned sets remain linked to the
1159 original mesh and are therefore not transferred to the copied mesh.
1161 To copy both the mesh and the corresponding geometry sets correctly,
1162 deep-copy them together.
1164 Returns:
1165 A deep copy of the mesh.
1166 """
1167 return _copy.deepcopy(self)