Coverage for src / sdynpy / core / sdynpy_geometry.py: 19%
2100 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1"""
2Objects and procedures to handle operations on test or model geometries
4This module defines a Geometry object as well as all of the subcomponents of
5a geometry object: nodes, elements, tracelines and coordinate system. Geometry
6plotting is also handled in this module.
7"""
9"""
10Copyright 2022 National Technology & Engineering Solutions of Sandia,
11LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
12Government retains certain rights in this software.
14This program is free software: you can redistribute it and/or modify
15it under the terms of the GNU General Public License as published by
16the Free Software Foundation, either version 3 of the License, or
17(at your option) any later version.
19This program is distributed in the hope that it will be useful,
20but WITHOUT ANY WARRANTY; without even the implied warranty of
21MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22GNU General Public License for more details.
24You should have received a copy of the GNU General Public License
25along with this program. If not, see <https://www.gnu.org/licenses/>.
26"""
28import copy
29import numpy as np
30import pyvista as pv
31import pyvistaqt as pvqt
32import pyqtgraph as pg
33import vtk
34from qtpy import uic, QtCore
35from qtpy.QtGui import QKeySequence
36try:
37 from qtpy.QtGui import QAction
38except ImportError:
39 from qtpy.QtWidgets import QAction
40from qtpy.QtWidgets import QApplication, QMainWindow, QSizePolicy, QFileDialog, QMenu
41from .sdynpy_array import SdynpyArray
42from .sdynpy_colors import colormap, coord_colormap
43from .sdynpy_coordinate import CoordinateArray, coordinate_array, from_nodelist
44from ..signal_processing.sdynpy_rotation import R, lstsq_rigid_transform
45from ..signal_processing.sdynpy_camera import point_on_pixel
46from .sdynpy_matrix import matrix
47import time
48import os
49from PIL import Image
50from scipy.spatial import Delaunay
51from scipy.linalg import block_diag
52import pandas as pd
53import openpyxl
54import warnings
55from scipy.io import savemat as scipy_savemat
57try:
58 repr(pv.GPUInfo())
59 IGNORE_PLOTS = False
60 # Adding this to help out non-interactive consoles
61 if QApplication.instance() is None:
62 app = QApplication([''])
64except RuntimeError:
65 IGNORE_PLOTS = True
66 print('No GPU Found, Geometry plotting will not work!')
68# Enumerations
70_cs_type_strings = {0: 'Cartesian',
71 1: 'Polar',
72 2: 'Spherical'}
74_element_types = {
75 11: 'Rod',
76 21: 'Linear beam',
77 22: 'Tapered beam',
78 23: 'Curved beam',
79 24: 'Parabolic beam',
80 31: 'Straight pipe',
81 32: 'Curved pipe',
82 41: 'Plane Stress Linear Triangle',
83 42: 'Plane Stress Parabolic Triangle',
84 43: 'Plane Stress Cubic Triangle',
85 44: 'Plane Stress Linear Quadrilateral',
86 45: 'Plane Stress Parabolic Quadrilateral',
87 46: 'Plane Strain Cubic Quadrilateral',
88 51: 'Plane Strain Linear Triangle',
89 52: 'Plane Strain Parabolic Triangle',
90 53: 'Plane Strain Cubic Triangle',
91 54: 'Plane Strain Linear Quadrilateral',
92 55: 'Plane Strain Parabolic Quadrilateral',
93 56: 'Plane Strain Cubic Quadrilateral',
94 61: 'Plate Linear Triangle',
95 62: 'Plate Parabolic Triangle',
96 63: 'Plate Cubic Triangle',
97 64: 'Plate Linear Quadrilateral',
98 65: 'Plate Parabolic Quadrilateral',
99 66: 'Plate Cubic Quadrilateral',
100 71: 'Membrane Linear Quadrilateral',
101 72: 'Membrane Parabolic Triangle',
102 73: 'Membrane Cubic Triangle',
103 74: 'Membrane Linear Triangle',
104 75: 'Membrane Parabolic Quadrilateral',
105 76: 'Membrane Cubic Quadrilateral',
106 81: 'Axisymetric Solid Linear Triangle',
107 82: 'Axisymetric Solid Parabolic Triangle',
108 84: 'Axisymetric Solid Linear Quadrilateral',
109 85: 'Axisymetric Solid Parabolic Quadrilateral',
110 91: 'Thin Shell Linear Triangle',
111 92: 'Thin Shell Parabolic Triangle',
112 93: 'Thin Shell Cubic Triangle',
113 94: 'Thin Shell Linear Quadrilateral',
114 95: 'Thin Shell Parabolic Quadrilateral',
115 96: 'Thin Shell Cubic Quadrilateral',
116 101: 'Thick Shell Linear Wedge',
117 102: 'Thick Shell Parabolic Wedge',
118 103: 'Thick Shell Cubic Wedge',
119 104: 'Thick Shell Linear Brick',
120 105: 'Thick Shell Parabolic Brick',
121 106: 'Thick Shell Cubic Brick',
122 111: 'Solid Linear Tetrahedron',
123 112: 'Solid Linear Wedge',
124 113: 'Solid Parabolic Wedge',
125 114: 'Solid Cubic Wedge',
126 115: 'Solid Linear Brick',
127 116: 'Solid Parabolic Brick',
128 117: 'Solid Cubic Brick',
129 118: 'Solid Parabolic Tetrahedron',
130 121: 'Rigid Bar',
131 122: 'Rigid Element',
132 136: 'Node To Node Translational Spring',
133 137: 'Node To Node Rotational Spring',
134 138: 'Node To Ground Translational Spring',
135 139: 'Node To Ground Rotational Spring',
136 141: 'Node To Node Damper',
137 142: 'Node To Gound Damper',
138 151: 'Node To Node Gap',
139 152: 'Node To Ground Gap',
140 161: 'Lumped Mass',
141 171: 'Axisymetric Linear Shell',
142 172: 'Axisymetric Parabolic Shell',
143 181: 'Constraint',
144 191: 'Plastic Cold Runner',
145 192: 'Plastic Hot Runner',
146 193: 'Plastic Water Line',
147 194: 'Plastic Fountain',
148 195: 'Plastic Baffle',
149 196: 'Plastic Rod Heater',
150 201: 'Linear node-to-node interface',
151 202: 'Linear edge-to-edge interface',
152 203: 'Parabolic edge-to-edge interface',
153 204: 'Linear face-to-face interface',
154 208: 'Parabolic face-to-face interface',
155 212: 'Linear axisymmetric interface',
156 213: 'Parabolic axisymmetric interface',
157 221: 'Linear rigid surface',
158 222: 'Parabolic rigid surface',
159 231: 'Axisymetric linear rigid surface',
160 232: 'Axisymentric parabolic rigid surface'}
162_exodus_elem_type_map = {'triangle':41,
163 'quadrilateral':44,
164 'hexahedral':115,
165 'tet':111,
166 'tetrahedral':111,
167 'sphere':161,
168 'spring':136,
169 'bar':21,
170 'bar2':21,
171 'bar3':23,
172 'beam':21,
173 'beam2':21,
174 'beam3':23,
175 'truss':21,
176 'truss2':21,
177 'truss3':23,
178 'quad':94,
179 'quad4':94,
180 'quad5':94,
181 'quad8':95,
182 'quad9':95,
183 'shell':94,
184 'shell4':94,
185 'shell8':95,
186 'shell9':95,
187 'tri':91,
188 'tri3':91,
189 'tri6':92,
190 'tri7':92,
191 'trishell':91,
192 'trishell3':91,
193 'trishell6':92,
194 'trishell7':92,
195 'hex':115,
196 'hex8':115,
197 'hex9':116,
198 'hex20':116,
199 'hex27':117,
200 'tetra':111,
201 'tetra4':111,
202 'tetra8':118,
203 'tetra10':118,
204 'tetra14':118,
205 'wedge':112,
206 'hexshell':116,
207 # 'pyramid'
208 # 'pyramid5'
209 # 'pyramid8'
210 # 'pyramid13'
211 # 'pyramid18'
212 # Add new elements here
213 }
215_beam_elem_types = [11, 21, 22, 23, 24, 31, 32, 121]
216_face_element_types = [v for v in range(41, 107)]
217_solid_element_types = [v for v in range(111, 119)]
218_vtk_element_map = {
219 111: vtk.VTK_TETRA,
220 112: vtk.VTK_WEDGE,
221 113: vtk.VTK_QUADRATIC_WEDGE,
222 115: vtk.VTK_HEXAHEDRON,
223 116: vtk.VTK_QUADRATIC_HEXAHEDRON,
224 118: vtk.VTK_QUADRATIC_TETRA,
225}
227_vtk_connectivity_reorder = {
228 vtk.VTK_QUADRATIC_HEXAHEDRON: [0, 1, 2, 3, 4, 5, 6,
229 7, 8, 9, 10, 11, 16, 17, 18, 19, 12, 13, 14, 15]
230}
232MAX_NUMBER_REPR = 100
235class GeometryPlotter(pvqt.BackgroundPlotter):
236 """Class used to plot geometry
238 This class is essentially identical to PyVista's BackgroundPlotter;
239 however a small amount of additional functionality is added."""
241 def __init__(self, *args, **kwargs):
242 super().__init__(*args, **kwargs)
243 for action in self.main_menu.actions():
244 if action.menu() and action.text().lower() == 'view':
245 view_menu = action.menu()
246 break
247 for action in view_menu.actions():
248 if action.menu() and action.text().lower() == 'orientation marker':
249 orien_menu = action.menu()
250 break
251 orien_menu.addSeparator()
252 orien_menu.addAction('Show Origin', self.show_origin_marker)
253 orien_menu.addAction('Hide Origin', self.hide_origin_marker)
254 self._origin_marker = None
256 def show_origin_marker(self, arrow_scale=0.05):
257 if self._origin_marker is None:
258 bbox_diagonal = self.length
259 arrow_factor = bbox_diagonal * arrow_scale
260 self._origin_marker = pv.tools.create_axes_marker(label_size=np.array([0.25, 0.1])*arrow_scale)
261 self._origin_marker.SetTotalLength(arrow_factor, arrow_factor, arrow_factor)
262 self.renderer.AddActor(self._origin_marker)
263 self.renderer.Modified()
265 def hide_origin_marker(self):
266 if self._origin_marker is not None:
267 self.renderer.remove_actor(self._origin_marker)
268 self.renderer.Modified()
270 def save_rotation_animation(self, filename=None, frames=20, frame_rate=20,
271 individual_images=False):
272 """
273 Saves an animation of the rotating geometry to a file
275 Parameters
276 ----------
277 filename : str, optional
278 Path to the file being saved. The default is None, which returns
279 the images as an array rather than saves them to a file
280 frames : int, optional
281 Number of frames in the animation. The default is 20.
282 frame_rate : float, optional
283 Number of frames per second if the animation is saved to a GIF.
284 The default is 20.
285 individual_images : bool, optional
286 Boolean to specify whether the images are saved as individual PNG
287 files or a single GIF. The default is False.
289 Returns
290 -------
291 imgs : np.ndarray or None
292 Returns array of images if filename is None, otherwise returns None
294 """
295 imgs = []
296 phases = np.linspace(0, 2 * np.pi, frames, endpoint=False)
297 focus = np.array(self.camera.focal_point)
298 distance = np.array(self.camera.distance)
299 initial_view_from = np.array(self.camera.position) - focus
300 view_up = self.camera.up
301 for phase in phases:
302 # Compute a rotation matrix
303 rotmat = R(view_up, phase)
304 view_from = (rotmat @ initial_view_from[:, np.newaxis]).squeeze()
305 self.camera.position = np.array(focus) + distance * \
306 np.array(view_from) / np.linalg.norm(view_from)
307 self.camera.focal_point = focus
308 self.camera.up = view_up
309 self.render()
310 # Force GUI Update to ensure the OS doesn't delay events
311 QApplication.processEvents()
312 imgs.append(self.screenshot())
313 if filename is None:
314 return imgs
315 else:
316 if individual_images:
317 for i, img in enumerate(imgs):
318 fname = os.path.splitext(filename)
319 fname = fname[0] + '-{:d}'.format(i) + fname[1]
320 Image.fromarray(img).save(fname)
321 else:
322 imgs = [Image.fromarray(img).convert(
323 'P', palette=Image.ADAPTIVE, colors=256) for img in imgs]
324 imgs[0].save(fp=filename, format='GIF', append_images=imgs[1:],
325 save_all=True, duration=int(1 / frame_rate * 1000), loop=0)
328class TransientPlotter(GeometryPlotter):
329 """Class used to plot transient deformations"""
331 def __init__(self, geometry, displacement_data, displacement_scale=1.0,
332 frames_per_second=20,
333 undeformed_opacity=0.0, deformed_opacity=1.0, plot_kwargs={},
334 transformation_shapes=None, num_curves=50,
335 show: bool = True,
336 app=None,
337 window_size=None,
338 off_screen=None,
339 allow_quit_keypress=True,
340 toolbar=True,
341 menu_bar=True,
342 editor=False,
343 update_app_icon=None,
344 **kwargs):
345 """
346 Create a TransientPlotter object to plot displacements over time
348 Parameters
349 ----------
350 geometry : Geometry
351 Geometry on which the displacements will be plotted
352 displacement_data : TimeHistoryArray
353 Transient displacement data that will be applied
354 displacement_scale : float, optional
355 Scale factor applied to displacements. The default is 1.0.
356 frames_per_second : float, optional
357 Number of time steps to plot per second while the displacement is
358 animating. Default is 20.
359 undeformed_opacity : float, optional
360 Opacity of the undeformed geometry. The default is 0.0, or
361 completely transparent.
362 deformed_opacity : float, optional
363 Opacity of the deformed geometry. The default is 1.0, or completely
364 opaque.
365 plot_kwargs : dict, optional
366 Keyword arguments passed to the Geometry.plot function
367 transformation_shapes : ShapeArray
368 Shape matrix that will be used to expand the data. Must be the
369 same size as the `displacement_data`
370 num_curves : int, optional
371 Maximum number of curves to plot on the time selector. Default is
372 50
373 show : bool, optional
374 Show the plotting window. If ``False``, show this window by
375 running ``show()``. The default is True.
376 app : QApplication, optional
377 Creates a `QApplication` if left as `None`. The default is None.
378 window_size : list of int, optional
379 Window size in pixels. Defaults to ``[1024, 768]``
380 off_screen : TYPE, optional
381 Renders off screen when True. Useful for automated
382 screenshots or debug testing. The default is None.
383 allow_quit_keypress : bool, optional
384 Allow user to exit by pressing ``"q"``. The default is True.
385 toolbar : bool, optional
386 If True, display the default camera toolbar. Defaults to True.
387 menu_bar : bool, optional
388 If True, display the default main menu. Defaults to True.
389 editor : TYPE, optional
390 If True, display the VTK object editor. Defaults to False.
391 update_app_icon : bool, optional
392 If True, update_app_icon will be called automatically to update the
393 Qt app icon based on the current rendering output. If None, the
394 logo of PyVista will be used. If False, no icon will be set.
395 Defaults to None. The default is None.
396 title : str, optional
397 Title of plotting window.
398 multi_samples : int, optional
399 The number of multi-samples used to mitigate aliasing. 4 is a
400 good default but 8 will have better results with a potential
401 impact on performance.
402 line_smoothing : bool, optional
403 If True, enable line smothing
404 point_smoothing : bool, optional
405 If True, enable point smothing
406 polygon_smoothing : bool, optional
407 If True, enable polygon smothing
408 auto_update : float, bool, optional
409 Automatic update rate in seconds. Useful for automatically
410 updating the render window when actors are change without
411 being automatically ``Modified``. If set to ``True``, update
412 rate will be 1 second.
414 Returns
415 -------
416 None.
418 """
419 self.shape_select_toolbar: pvqt.plotting.QToolBar = None
420 self.animation_control_toolbar: pvqt.plotting.QToolBar = None
421 self.loop_action = None
422 self.plot_abscissa_action = None
423 self.timer: pvqt.plotting.QTimer = None
424 self.frames_per_second = frames_per_second
425 self.animation_speed = 1.0
426 self.last_time = time.time()
428 super(TransientPlotter, self).__init__(show, app, None if window_size is None else tuple(window_size),
429 off_screen,
430 allow_quit_keypress, toolbar,
431 menu_bar, editor, update_app_icon,
432 **kwargs)
434 self.text = self.add_text('', font_size=10)
435 self.loop_action.setChecked(False)
436 self.plot_abscissa_action.setChecked(True)
438 # Get the layout so we can add a new widget to it
439 layout = self.frame.layout()
440 self.time_selector_plot = pg.PlotWidget()
441 layout.addWidget(self.time_selector_plot)
442 sp = self.time_selector_plot.sizePolicy()
443 sp.setVerticalPolicy(QSizePolicy.Minimum)
444 self.time_selector_plot.setSizePolicy(sp)
445 self.time_selector_plot.sizeHint = lambda: QtCore.QSize(1006, 80)
446 # self.time_selector_plot.getPlotItem().hideAxis('left')
447 # self.time_selector_plot.getPlotItem().axes['left']['item'].setStyle(showValues=False)
449 # Now we need to plot the data in the plot
450 if num_curves < displacement_data.size:
451 indices = np.linspace(0, displacement_data.size - 1, num_curves).astype(int)
452 else:
453 indices = np.arange(displacement_data.size)
454 for data in displacement_data.flatten()[indices]:
455 self.time_selector_plot.plot(data.abscissa, data.ordinate)
456 # Plot the envelope
457 max_envelope = displacement_data.flatten().ordinate.max(axis=0)
458 min_envelope = displacement_data.flatten().ordinate.min(axis=0)
459 self.time_selector_plot.plot(displacement_data.flatten(
460 )[0].abscissa, max_envelope, pen=pg.mkPen(color='k', width=2))
461 self.time_selector_plot.plot(displacement_data.flatten(
462 )[0].abscissa, min_envelope, pen=pg.mkPen(color='k', width=2))
464 # Add an infinite line for a cursor
465 self.time_selector = pg.InfiniteLine(
466 0, movable=True, bounds=[displacement_data.abscissa.min(),
467 displacement_data.abscissa.max()])
468 self.time_selector_plot.addItem(self.time_selector)
469 self.time_selector.sigPositionChanged.connect(self.modify_abscissa)
471 # Set initial values
472 self.displacement_scale = displacement_scale
473 self.abscissa_index = 0
474 self.displacement_data = displacement_data.flatten()
475 self.transformation_shapes = None if transformation_shapes is None else transformation_shapes.flatten()
476 self.plot_geometry = geometry
477 self.num_abscissa = self.displacement_data.num_elements
478 self.abscissa = self.displacement_data[0].abscissa
479 self.abscissa_sig_figs = int(np.ceil(-np.log10(np.min(np.diff(self.abscissa))))) + 1
480 self.abscissa_format_spec = 'Time: {{:0.{:}f}}'.format(
481 self.abscissa_sig_figs if self.abscissa_sig_figs > 0 else 0)
483 # Plot the geometry
484 plotter, self.face_mesh_undeformed, self.point_mesh_undeformed, self.solid_mesh_undeformed = geometry.plot(plotter=self,
485 opacity=undeformed_opacity,
486 **plot_kwargs)
488 plotter, self.face_mesh_deformed, self.point_mesh_deformed, self.solid_mesh_deformed = geometry.plot(plotter=self,
489 opacity=deformed_opacity,
490 **plot_kwargs)
492 # Compute the node indices up front to make things faster when switching steps
493 if self.transformation_shapes is None:
494 coordinates = np.unique(self.displacement_data.coordinate)
495 else:
496 coordinates = np.unique(self.transformation_shapes.coordinate)
497 coordinates = coordinates[
498 np.isin(coordinates.node, self.plot_geometry.node.id)
499 & (np.abs(coordinates.direction) > 0)
500 & (np.abs(coordinates.direction) < 4)
501 ]
502 global_deflections = self.plot_geometry.global_deflection(coordinates).T
503 if self.transformation_shapes is None:
504 shape_array = self.displacement_data[coordinates[:, np.newaxis]].ordinate
505 else:
506 shape_array = self.transformation_shapes[coordinates].T
507 # Compute global displacements for each node
508 self.node_deflections = np.zeros((3 * self.plot_geometry.node.size, shape_array.shape[-1]))
509 for node_index, node in self.plot_geometry.node.ndenumerate():
510 coord_indices = coordinates.node == node.id
511 self.node_deflections[node_index[0] * 3:(node_index[0] + 1) *
512 3] = global_deflections[:, coord_indices] @ shape_array[coord_indices]
513 if self.transformation_shapes is None:
514 self.weighting = None
515 else:
516 self.weighting = self.displacement_data.ordinate
517 self.update_displacement()
519 def modify_abscissa(self):
520 abscissa = self.time_selector.value()
521 self.abscissa_index = np.argmin(np.abs(self.abscissa - abscissa))
522 self.update_displacement(update_selector=False)
524 def update_displacement(self, update_selector=True):
525 index = int(round(self.abscissa_index))
526 abscissa = self.abscissa[index]
527 if self.plot_abscissa_action.isChecked():
528 self.text.SetText(2, self.abscissa_format_spec.format(abscissa))
529 else:
530 self.text.SetText(2, '')
531 if self.weighting is None:
532 displacement = self.node_deflections[:, index] * self.displacement_scale
533 else:
534 displacement = self.node_deflections @ self.weighting[:,
535 np.newaxis, index] * self.displacement_scale
536 for mesh_deformed, mesh_undeformed in (
537 (self.face_mesh_deformed, self.face_mesh_undeformed),
538 (self.point_mesh_deformed, self.point_mesh_undeformed),
539 (self.solid_mesh_deformed, self.solid_mesh_undeformed)
540 ):
541 if mesh_deformed is not None:
542 mesh_deformed.points = (mesh_undeformed.points +
543 displacement.reshape(-1, 3) * self.displacement_scale
544 )
545 # Update the time selector
546 if update_selector:
547 self.time_selector.blockSignals(True)
548 self.time_selector.setValue(abscissa)
549 self.time_selector.blockSignals(False)
550 self.render()
552 def add_toolbars(self) -> None:
553 """
554 Adds toolbars to the BackgroundPlotter
556 """
557 super().add_toolbars()
558 self.shape_select_toolbar = self.app_window.addToolBar(
559 "Shape Selector"
560 )
562 self._add_action(
563 self.shape_select_toolbar, '|<', self.go_to_first_step
564 )
565 self._add_action(
566 self.shape_select_toolbar, '<', self.go_to_previous_step
567 )
568 self._add_action(
569 self.shape_select_toolbar, '>', self.go_to_next_step
570 )
571 self._add_action(
572 self.shape_select_toolbar, '>|', self.go_to_last_step
573 )
575 self.animation_control_toolbar = self.app_window.addToolBar(
576 "Animation Control"
577 )
579 self._add_action(
580 self.animation_control_toolbar, '< Play', self.play_animation_reverse
581 )
582 self._add_action(
583 self.animation_control_toolbar, 'Stop', self.stop_animation
584 )
585 self._add_action(
586 self.animation_control_toolbar, 'Play >', self.play_animation
587 )
589 def add_menu_bar(self) -> None:
590 """
591 Adds the menu bar to the BackgroundPlotter
593 """
594 super().add_menu_bar()
596 shape_menu = self.main_menu.addMenu("Shape")
597 scaling_menu = shape_menu.addMenu('Scaling')
598 scaling_menu.addAction('1/4x', self.select_scaling_0p25)
599 scaling_menu.addAction('1/2x', self.select_scaling_0p5)
600 self.scale_0p8_action = QAction('0.8x')
601 self.scale_0p8_action.triggered.connect(self.select_scaling_0p8)
602 self.scale_0p8_action.setShortcut(QKeySequence(','))
603 scaling_menu.addAction(self.scale_0p8_action)
604 self.scale_reset_action = QAction('Reset')
605 self.scale_reset_action.triggered.connect(self.select_scaling_1)
606 self.scale_reset_action.setShortcut(QKeySequence('/'))
607 scaling_menu.addAction(self.scale_reset_action)
608 self.scale_1p25_action = QAction('1.25x')
609 self.scale_1p25_action.triggered.connect(self.select_scaling_1p25)
610 self.scale_1p25_action.setShortcut(QKeySequence('.'))
611 scaling_menu.addAction(self.scale_1p25_action)
612 scaling_menu.addAction('2x', self.select_scaling_2p0)
613 scaling_menu.addAction('4x', self.select_scaling_4p0)
614 speed_menu = shape_menu.addMenu('Animation Speed')
615 self.speed_0p8_action = QAction('0.8x')
616 self.speed_0p8_action.triggered.connect(self.select_speed_0p8)
617 self.speed_0p8_action.setShortcut(QKeySequence(';'))
618 speed_menu.addAction(self.speed_0p8_action)
619 speed_menu.addAction('Reset', self.select_speed_1)
620 self.speed_1p25_action = QAction('1.25x')
621 self.speed_1p25_action.triggered.connect(self.select_speed_1p25)
622 self.speed_1p25_action.setShortcut(QKeySequence("'"))
623 speed_menu.addAction(self.speed_1p25_action)
624 shape_menu.addSeparator()
625 self.plot_abscissa_action = pvqt.plotting.QAction('Show Abscissa', checkable=True)
626 shape_menu.addAction(self.plot_abscissa_action)
627 self.loop_action = pvqt.plotting.QAction('Loop Animation', checkable=True)
628 shape_menu.addAction(self.loop_action)
630 def go_to_first_step(self):
631 self.stop_animation()
632 self.abscissa_index = 0
633 self.update_displacement()
635 def go_to_previous_step(self):
636 self.stop_animation()
637 self.abscissa_index -= 1
638 if self.abscissa_index < 0:
639 if self.loop_action.isChecked():
640 self.abscissa_index = self.num_abscissa - 1
641 else:
642 self.abscissa_index = 0
643 self.update_displacement()
645 def go_to_next_step(self):
646 self.stop_animation()
647 self.abscissa_index += 1
648 if self.abscissa_index > self.num_abscissa - 1:
649 if self.loop_action.isChecked():
650 self.abscissa_index = 0
651 else:
652 self.abscissa_index = self.num_abscissa - 1
653 self.update_displacement()
655 def go_to_last_step(self):
656 self.stop_animation()
657 self.abscissa_index = self.num_abscissa - 1
658 self.update_displacement()
660 def select_scaling_0p25(self):
661 """
662 Adjusts the current shape scaling by 0.25x
663 """
664 self.displacement_scale *= 0.25
665 self.last_time = time.time()
666 self.update_displacement()
668 def select_scaling_0p5(self):
669 """
670 Adjusts the current shape scaling by 0.5x
671 """
672 self.displacement_scale *= 0.5
673 self.last_time = time.time()
674 self.update_displacement()
676 def select_scaling_0p8(self):
677 """
678 Adjusts the current shape scaling by 0.8x
679 """
680 self.displacement_scale *= 0.8
681 self.last_time = time.time()
682 self.update_displacement()
684 def select_scaling_1(self):
685 """
686 Resets shape scaling to 1
687 """
688 self.displacement_scale = 1.0
689 self.last_time = time.time()
690 self.update_displacement()
692 def select_scaling_1p25(self):
693 """
694 Adjusts the current shape scaling by 1.25x
695 """
696 self.displacement_scale *= 1.25
697 self.last_time = time.time()
698 self.update_displacement()
700 def select_scaling_2p0(self):
701 """
702 Adjusts the current shape scaling by 2x
703 """
704 self.displacement_scale *= 2.0
705 self.last_time = time.time()
706 self.update_displacement()
708 def select_scaling_4p0(self):
709 """
710 Adjusts the current shape scaling by 4x
711 """
712 self.displacement_scale *= 4.0
713 self.last_time = time.time()
714 self.update_displacement()
716 def select_speed_1(self):
717 """
718 Resets animation speed to default quantities
719 """
720 self.animation_speed = 1.0
722 def select_speed_0p8(self):
723 """
724 Adjusts the current animation speed by 0.8
725 """
726 self.animation_speed *= 0.8
728 def select_speed_1p25(self):
729 """
730 Adjusts the current animation speed by 1.25
731 """
732 self.animation_speed *= 1.25
734 def play_animation(self):
735 self.timer = pvqt.plotting.QTimer()
736 self.timer.timeout.connect(self.update_abscissa)
737 self.timer.start(10)
738 self.last_time = time.time()
740 def play_animation_reverse(self):
741 self.timer = pvqt.plotting.QTimer()
742 self.timer.timeout.connect(self.update_abscissa_reverse)
743 self.timer.start(10)
744 self.last_time = time.time()
746 def stop_animation(self):
747 if self.timer is not None:
748 self.timer.stop()
749 self.timer = None
751 def update_abscissa(self):
752 now = time.time()
753 time_elapsed = now - self.last_time
754 index_change = self.frames_per_second * self.animation_speed * time_elapsed
755 self.last_time = now
756 self.abscissa_index += index_change
757 if self.abscissa_index > self.num_abscissa - 1:
758 if self.loop_action.isChecked():
759 self.abscissa_index = 0
760 else:
761 self.stop_animation()
762 self.abscissa_index = self.num_abscissa - 1
763 self.update_displacement()
765 def update_abscissa_reverse(self):
766 now = time.time()
767 time_elapsed = now - self.last_time
768 index_change = self.frames_per_second * self.animation_speed * time_elapsed
769 self.last_time = now
770 self.abscissa_index -= index_change
771 if self.abscissa_index < 0:
772 if self.loop_action.isChecked():
773 self.abscissa_index = self.num_abscissa - 1
774 else:
775 self.stop_animation()
776 self.abscissa_index = 0
777 self.update_displacement()
779 def save_animation(self, filename=None, frame_rate=20,
780 individual_images=False, start_time=None, stop_time=None,
781 step=None):
782 """
783 Saves the current shape animation to a file
785 Parameters
786 ----------
787 filename : str, optional
788 Path to the file being saved. The default is None, which returns
789 the images as an array rather than saves them to a file
790 frame_rate : float, optional
791 Number of frames per second if the animation is saved to a GIF.
792 The default is 20.
793 individual_images : bool, optional
794 Boolean to specify whether the images are saved as individual PNG
795 files or a single GIF. The default is False.
796 start_time : float, optional
797 Time value at which the animation will start. Default is first
798 time step
799 stop_time : float, optional
800 Time value at which the animation will end. Default is last time
801 step
802 step : int, optional
803 Only render every ``step`` time steps. Default is to render all
804 timesteps
805 Returns
806 -------
807 imgs : np.ndarray or None
808 Returns array of images if filename is None, otherwise returns None
810 """
811 self.stop_animation()
812 imgs = []
813 indices = np.arange(self.abscissa.size)
814 if start_time is not None:
815 indices = indices[self.abscissa[indices] >= start_time]
816 if stop_time is not None:
817 indices = indices[self.abscissa[indices] <= stop_time]
818 if step is not None:
819 indices = indices[::step]
820 for index in indices:
821 self.abscissa_index = index
822 self.update_displacement()
823 # Force GUI Update to ensure the OS doesn't delay events
824 QApplication.processEvents()
825 imgs.append(self.screenshot())
826 if filename is None:
827 return imgs
828 else:
829 if individual_images:
830 for i, img in enumerate(imgs):
831 fname = os.path.splitext(filename)
832 fname = fname[0] + '-{:d}'.format(i) + fname[1]
833 Image.fromarray(img).save(fname)
834 else:
835 imgs = [Image.fromarray(img).convert(
836 'P', palette=Image.ADAPTIVE, colors=256) for img in imgs]
837 imgs[0].save(fp=filename, format='GIF', append_images=imgs[1:],
838 save_all=True, duration=int(1 / frame_rate * 1000), loop=0)
840 def _close(self):
841 self.stop_animation()
842 super(TransientPlotter, self)._close()
845class ShapePlotter(GeometryPlotter):
846 """Class used to plot animated shapes"""
848 def __init__(self, geometry, shapes, plot_kwargs,
849 background_plotter_kwargs={'auto_update': 0.01},
850 undeformed_opacity=0.25, starting_scale=1.0,
851 deformed_opacity=1.0, shape_name='Mode',
852 show_damping=True):
853 """
854 Create a Shape Plotter object to plot shapes
856 Parameters
857 ----------
858 geometry : Geometry
859 Geometry on which the shapes will be plotted
860 shapes : ShapeArray
861 Shapes to plot on the geometry
862 plot_kwargs : dict
863 Keyword arguments passed to the Geometry.plot function
864 background_plotter_kwargs : dict, optional
865 Keyword arguments passed to the BackgroundPlotter constructor.
866 The default is {'auto_update':0.01}.
867 undeformed_opacity : float, optional
868 Opacity of the undeformed geometry. Set to zero for no undeformed
869 geometry. The default is 0.25.
870 starting_scale : float, optional
871 Starting scale of the shapes on the plot. The default is 1.0.
872 shape_name : str, optional
873 Name that will be used in the plotter to describe the shape. The
874 default is "Mode"
875 show_damping : bool, optional
876 Boolean that specifies whether the damping should be displayed in
877 the comment.
878 """
879 self.shape_select_toolbar: pvqt.plotting.QToolBar = None
880 self.animation_control_toolbar: pvqt.plotting.QToolBar = None
881 self.signed_amplitude_action = None
882 self.complex_action = None
883 self.real_action = None
884 self.imag_action = None
885 self.loop_action = None
886 self.plot_comments_action = None
887 self.save_animation_action = None
888 self.display_undeformed = None
889 self.geometry = geometry
890 self.shapes = shapes.ravel()
891 self.node_displacements = None
892 self.current_shape = 0
893 self.displacement_scale = 1.0
894 self.displacement_scale_modification = starting_scale
895 self.animation_speed = 1.0
896 self.last_time = time.time()
897 self.phase = 0.0
898 self.timer: pvqt.plotting.QTimer = None
899 self.shape_name = shape_name
900 self.show_damping = show_damping
901 self.shape_comment_table = None
903 super().__init__(**background_plotter_kwargs)
905 self.text = self.add_text('', font_size=10)
907 self.complex_action.setChecked(True)
908 self.loop_action.setChecked(True)
909 self.plot_comments_action.setChecked(True)
911 plotter, self.face_mesh_undeformed, self.point_mesh_undeformed, self.solid_mesh_undeformed = geometry.plot(plotter=self,
912 opacity=undeformed_opacity,
913 **plot_kwargs)
915 plotter, self.face_mesh_deformed, self.point_mesh_deformed, self.solid_mesh_deformed = geometry.plot(plotter=self,
916 opacity=deformed_opacity,
917 **plot_kwargs)
919 # Compute the node indices up front to make things faster when switching shapes
920 coordinates = self.shapes[self.current_shape].coordinate
921 nodes = coordinates.node
922 # Only do coordinates that are in the geometry node array
923 self.indices = np.isin(coordinates.node, self.geometry.node.id)
924 coordinates = coordinates[self.indices]
925 nodes = nodes[self.indices]
926 direction = coordinates.direction
927 self.coordinate_node_index_map = {}
928 for coordinate_index, node_id in np.ndenumerate(nodes):
929 if direction[coordinate_index] == 0 or abs(direction[coordinate_index]) > 3:
930 continue
931 if node_id not in self.coordinate_node_index_map:
932 self.coordinate_node_index_map[node_id] = []
933 self.coordinate_node_index_map[node_id].append(coordinate_index[0])
934 local_deformations = coordinates.local_direction()
935 ordered_nodes = self.geometry.node(nodes)
936 coordinate_systems = self.geometry.coordinate_system(ordered_nodes.disp_cs)
937 points = global_coord(self.geometry.coordinate_system(
938 ordered_nodes.def_cs), ordered_nodes.coordinate)
939 self.global_deflections = global_deflection(coordinate_systems, local_deformations, points)
941 self.compute_displacements()
942 self.show_comment()
943 self.render()
944 self.update_shape()
945 self.update_shape_mode(0)
946 QApplication.processEvents()
948 def compute_displacements(self, compute_scale=True) -> np.ndarray:
949 """
950 Computes displacements to apply to the geometry
952 Parameters
953 ----------
954 compute_scale : bool, optional
955 Renormalize the displacement scaling. The default is True.
957 """
958 shape_displacements = self.shapes[self.current_shape].shape_matrix[..., self.indices]
959 node_displacements = np.zeros(self.geometry.node.shape + (3,),
960 dtype=shape_displacements.dtype)
961 for node_index, node in self.geometry.node.ndenumerate():
962 try:
963 node_indices = self.coordinate_node_index_map[node.id]
964 except KeyError:
965 print('Node {:} not found in shape array'.format(node.id))
966 node_displacements[node_index] = 0
967 continue
968 node_deflection_scales = shape_displacements[node_indices]
969 node_deflection_directions = self.global_deflections[node_indices]
970 node_displacements[node_index] += np.sum(node_deflection_scales[:, np.newaxis] *
971 node_deflection_directions, axis=0)
972 self.node_displacements = node_displacements.reshape(-1, 3)
973 # Compute maximum displacement
974 if compute_scale:
975 max_displacement = np.max(np.linalg.norm(self.node_displacements, axis=-1))
976 global_coords = self.geometry.global_node_coordinate()
977 bbox_diagonal = np.linalg.norm(
978 np.max(global_coords, axis=0) - np.min(global_coords, axis=0))
979 self.displacement_scale = 0.05 * bbox_diagonal / \
980 max_displacement # Set to 5% of bounding box dimension
982 def update_shape_mode(self, phase=None):
983 """
984 Updates the mode that is being plotted
986 Parameters
987 ----------
988 phase : float, optional
989 Sets the current phase of the shape. The default is None, which
990 computes the new phase based on the time elapsed.
992 """
993 if phase is None:
994 now = time.time()
995 time_elapsed = now - self.last_time
996 phase_change = 2 * np.pi / self.animation_speed * time_elapsed
997 self.last_time = now
998 self.phase += phase_change
999 else:
1000 self.phase = phase
1001 if self.phase > 2 * np.pi:
1002 self.phase -= 2 * np.pi
1003 if not self.loop_action.isChecked():
1004 self.stop_animation()
1005 # print('Updating Phase: {:}'.format(self.phase))
1006 if self.complex_action.isChecked():
1007 # print('Complex')
1008 deformation = np.real(np.exp(self.phase * 1j) * self.node_displacements)
1009 elif self.real_action.isChecked():
1010 # print('Real')
1011 deformation = np.cos(self.phase) * np.real(self.node_displacements)
1012 elif self.imag_action.isChecked():
1013 # print('Imag')
1014 deformation = np.cos(self.phase) * np.imag(self.node_displacements)
1015 elif self.signed_amplitude_action.isChecked():
1016 # print('Signed Amplitude')
1017 imag = np.imag(self.node_displacements).flatten()[:, np.newaxis]
1018 real = np.real(self.node_displacements).flatten()[:, np.newaxis]
1019 best_fit = np.linalg.lstsq(real, imag)[0][0, 0]
1020 best_fit_phase = np.arctan(best_fit)
1021 rectified_displacements = self.node_displacements / np.exp(1j * best_fit_phase)
1022 deformation = np.cos(self.phase) * np.real(rectified_displacements)
1023 for mesh_deformed, mesh_undeformed in (
1024 (self.face_mesh_deformed, self.face_mesh_undeformed),
1025 (self.point_mesh_deformed, self.point_mesh_undeformed),
1026 (self.solid_mesh_deformed, self.solid_mesh_undeformed)
1027 ):
1028 if mesh_deformed is not None:
1029 mesh_deformed.points = (mesh_undeformed.points +
1030 deformation * self.displacement_scale * self.displacement_scale_modification
1031 )
1032 self.render()
1034 def update_shape(self):
1035 """
1036 Updates the shape that is being plotted
1037 """
1038 self.update_shape_mode()
1040 def save_animation_from_action(self):
1041 filename, file_filter = QFileDialog.getSaveFileName(
1042 self, 'Select File to Save Animation', filter='GIF (*.gif)')
1043 if filename == '':
1044 return
1045 self.save_animation(filename)
1047 def save_animation(self, filename=None, frames=20, frame_rate=20,
1048 individual_images=False):
1049 """
1050 Saves the current shape animation to a file
1052 Parameters
1053 ----------
1054 filename : str, optional
1055 Path to the file being saved. The default is None, which returns
1056 the images as an array rather than saves them to a file
1057 frames : int, optional
1058 Number of frames in the animation. The default is 20.
1059 frame_rate : float, optional
1060 Number of frames per second if the animation is saved to a GIF.
1061 The default is 20.
1062 individual_images : bool, optional
1063 Boolean to specify whether the images are saved as individual PNG
1064 files or a single GIF. The default is False.
1066 Returns
1067 -------
1068 imgs : np.ndarray or None
1069 Returns array of images if filename is None, otherwise returns None
1071 """
1072 self.stop_animation()
1073 imgs = []
1074 phases = np.linspace(0, 2 * np.pi, frames, endpoint=False)
1075 for phase in phases:
1076 self.update_shape_mode(phase)
1077 # Force GUI Update to ensure the OS doesn't delay events
1078 QApplication.processEvents()
1079 imgs.append(self.screenshot())
1080 if filename is None:
1081 return imgs
1082 else:
1083 if individual_images:
1084 for i, img in enumerate(imgs):
1085 fname = os.path.splitext(filename)
1086 fname = fname[0] + '-{:d}'.format(i) + fname[1]
1087 Image.fromarray(img).save(fname)
1088 else:
1089 imgs = [Image.fromarray(img).convert(
1090 'P', palette=Image.ADAPTIVE, colors=256) for img in imgs]
1091 imgs[0].save(fp=filename, format='GIF', append_images=imgs[1:],
1092 save_all=True, duration=int(1 / frame_rate * 1000), loop=0)
1094 def save_animation_all_shapes(self, filename_base='Shape_{:}.gif',
1095 frames=20, frame_rate=20,
1096 individual_images=False):
1097 """
1098 Helper function to easily save animations for all shapes
1100 Parameters
1101 ----------
1102 filename_base : str, optional
1103 Filename that will be passed to the format function to produce the
1104 actual file name for each shape. The default is 'Shape_{:}.gif'.
1105 frames : int, optional
1106 Number of frames in the animation. The default is 20.
1107 frame_rate : float, optional
1108 Number of frames per second if the animation is saved to a GIF.
1109 The default is 20.
1110 individual_images : bool, optional
1111 Boolean to specify whether the images are saved as individual PNG
1112 files or a single GIF. The default is False.
1114 """
1115 for shape in range(self.shapes.size):
1116 self.current_shape = shape
1117 self.compute_displacements()
1118 self.show_comment()
1119 self.save_animation(filename_base.format(shape + 1),
1120 frames, frame_rate, individual_images)
1122 def add_menu_bar(self) -> None:
1123 """
1124 Adds the menu bar to the BackgroundPlotter
1126 """
1127 super().add_menu_bar()
1129 # Get the file menu
1130 file_menu = [child for child in self.main_menu.children()
1131 if child.__class__.__name__ == 'QMenu'
1132 and child.title() == 'File'][0]
1133 self.save_animation_action = pvqt.plotting.QAction('Save Animation')
1134 self.save_animation_action.triggered.connect(self.save_animation_from_action)
1135 file_menu.insertAction(file_menu.actions()[1], self.save_animation_action)
1137 shape_menu = self.main_menu.addMenu("Shape")
1138 complex_display_menu = shape_menu.addMenu("Complex Display")
1139 self.signed_amplitude_action = pvqt.plotting.QAction('Signed Amplitude', checkable=True)
1140 self.signed_amplitude_action.triggered.connect(self.select_complex)
1141 complex_display_menu.addAction(self.signed_amplitude_action)
1142 self.complex_action = pvqt.plotting.QAction('Complex', checkable=True)
1143 self.complex_action.triggered.connect(self.select_complex)
1144 complex_display_menu.addAction(self.complex_action)
1145 self.real_action = pvqt.plotting.QAction('Real', checkable=True)
1146 self.real_action.triggered.connect(self.select_complex)
1147 complex_display_menu.addAction(self.real_action)
1148 self.imag_action = pvqt.plotting.QAction('Imaginary', checkable=True)
1149 self.imag_action.triggered.connect(self.select_complex)
1150 complex_display_menu.addAction(self.imag_action)
1151 scaling_menu = shape_menu.addMenu('Scaling')
1152 scaling_menu.addAction('1/4x', self.select_scaling_0p25)
1153 scaling_menu.addAction('1/2x', self.select_scaling_0p5)
1154 self.scale_0p8_action = QAction('0.8x')
1155 self.scale_0p8_action.triggered.connect(self.select_scaling_0p8)
1156 self.scale_0p8_action.setShortcut(QKeySequence(','))
1157 scaling_menu.addAction(self.scale_0p8_action)
1158 self.scale_reset_action = QAction('Reset')
1159 self.scale_reset_action.triggered.connect(self.select_scaling_1)
1160 self.scale_reset_action.setShortcut(QKeySequence('/'))
1161 scaling_menu.addAction(self.scale_reset_action)
1162 self.scale_1p25_action = QAction('1.25x')
1163 self.scale_1p25_action.triggered.connect(self.select_scaling_1p25)
1164 self.scale_1p25_action.setShortcut(QKeySequence('.'))
1165 scaling_menu.addAction(self.scale_1p25_action)
1166 scaling_menu.addAction('2x', self.select_scaling_2p0)
1167 scaling_menu.addAction('4x', self.select_scaling_4p0)
1168 speed_menu = shape_menu.addMenu('Animation Speed')
1169 self.speed_0p8_action = QAction('0.8x')
1170 self.speed_0p8_action.triggered.connect(self.select_speed_0p8)
1171 self.speed_0p8_action.setShortcut(QKeySequence(';'))
1172 speed_menu.addAction(self.speed_0p8_action)
1173 speed_menu.addAction('Reset', self.select_speed_1)
1174 self.speed_1p25_action = QAction('1.25x')
1175 self.speed_1p25_action.triggered.connect(self.select_speed_1p25)
1176 self.speed_1p25_action.setShortcut(QKeySequence("'"))
1177 speed_menu.addAction(self.speed_1p25_action)
1178 shape_menu.addSeparator()
1179 self.plot_comments_action = pvqt.plotting.QAction('Show Comments', checkable=True)
1180 self.plot_comments_action.triggered.connect(self.show_comment)
1181 shape_menu.addAction(self.plot_comments_action)
1182 self.loop_action = pvqt.plotting.QAction('Loop Animation', checkable=True)
1183 self.loop_action.triggered.connect(self.select_loop)
1184 shape_menu.addAction(self.loop_action)
1186 display_menu = self.main_menu.addMenu('Display')
1187 self.display_undeformed = pvqt.plotting.QAction('Show Undeformed', checkable=True)
1188 self.display_undeformed.triggered.connect(self.toggle_undeformed)
1189 display_menu.addAction(self.display_undeformed)
1191 def add_toolbars(self) -> None:
1192 """
1193 Adds toolbars to the BackgroundPlotter
1195 """
1196 super().add_toolbars()
1197 self.shape_select_toolbar = self.app_window.addToolBar(
1198 "Shape Selector"
1199 )
1201 self._add_action(
1202 self.shape_select_toolbar, '<<', self.prev_shape
1203 )
1204 self._add_action(
1205 self.shape_select_toolbar, '?', self.select_shape
1206 )
1207 self._add_action(
1208 self.shape_select_toolbar, '>>', self.next_shape
1209 )
1211 self.animation_control_toolbar = self.app_window.addToolBar(
1212 "Animation Control"
1213 )
1215 self._add_action(
1216 self.animation_control_toolbar, 'Play', self.play_animation
1217 )
1218 self._add_action(
1219 self.animation_control_toolbar, 'Stop', self.stop_animation
1220 )
1222 def play_animation(self):
1223 """Starts the animation playing"""
1224 self.timer = pvqt.plotting.QTimer()
1225 self.timer.timeout.connect(self.update_shape)
1226 self.timer.start(10)
1227 self.phase = 0.0
1228 self.last_time = time.time()
1230 def stop_animation(self):
1231 """Stops the animation from playing"""
1232 if self.timer is not None:
1233 self.timer.stop()
1234 self.timer = None
1236 def toggle_undeformed(self) -> None:
1237 pass
1239 def select_scaling_0p25(self):
1240 """
1241 Adjusts the current shape scaling by 0.25x
1242 """
1243 self.displacement_scale_modification *= 0.25
1244 self.phase = 0.0
1245 self.last_time = time.time()
1246 self.update_shape()
1248 def select_scaling_0p5(self):
1249 """
1250 Adjusts the current shape scaling by 0.5x
1251 """
1252 self.displacement_scale_modification *= 0.5
1253 self.phase = 0.0
1254 self.last_time = time.time()
1255 self.update_shape()
1257 def select_scaling_0p8(self):
1258 """
1259 Adjusts the current shape scaling by 0.8x
1260 """
1261 self.displacement_scale_modification *= 0.8
1262 self.phase = 0.0
1263 self.last_time = time.time()
1264 self.update_shape()
1266 def select_scaling_1(self):
1267 """
1268 Resets shape scaling to 1
1269 """
1270 self.displacement_scale_modification = 1.0
1271 self.phase = 0.0
1272 self.last_time = time.time()
1273 self.update_shape()
1275 def select_scaling_1p25(self):
1276 """
1277 Adjusts the current shape scaling by 1.25x
1278 """
1279 self.displacement_scale_modification *= 1.25
1280 self.phase = 0.0
1281 self.last_time = time.time()
1282 self.update_shape()
1284 def select_scaling_2p0(self):
1285 """
1286 Adjusts the current shape scaling by 2x
1287 """
1288 self.displacement_scale_modification *= 2.0
1289 self.phase = 0.0
1290 self.last_time = time.time()
1291 self.update_shape()
1293 def select_scaling_4p0(self):
1294 """
1295 Adjusts the current shape scaling by 4x
1296 """
1297 self.displacement_scale_modification *= 4.0
1298 self.phase = 0.0
1299 self.last_time = time.time()
1300 self.update_shape()
1302 def select_speed_1(self):
1303 """
1304 Resets animation speed to default quantities
1305 """
1306 self.animation_speed = 1.0
1308 def select_speed_0p8(self):
1309 """
1310 Adjusts the current animation speed by 0.8
1311 """
1312 self.animation_speed /= 0.8
1314 def select_speed_1p25(self):
1315 """
1316 Adjusts the current animation speed by 1.25
1317 """
1318 self.animation_speed /= 1.25
1320 def select_complex(self) -> None:
1321 """
1322 Adjusts how complex shapes are plotted
1323 """
1324 complex_actions = [
1325 self.signed_amplitude_action,
1326 self.complex_action,
1327 self.real_action,
1328 self.imag_action
1329 ]
1330 sender = self.sender()
1331 sender_index = complex_actions.index(sender)
1332 complex_actions.pop(sender_index)
1333 if not sender.isChecked():
1334 sender.setChecked(True)
1335 else:
1336 for action in complex_actions:
1337 action.setChecked(False)
1338 self.phase = 0.0
1339 self.last_time = time.time()
1340 self.update_shape()
1342 def select_loop(self) -> None:
1343 pass
1345 def next_shape(self) -> None:
1346 """
1347 Increases the index of the shape being plotted
1348 """
1349 self.current_shape += 1
1350 if self.current_shape >= self.shapes.size:
1351 self.current_shape = 0
1352 self.compute_displacements()
1353 self.phase = 0.0
1354 self.last_time = time.time()
1355 self.update_shape()
1356 self.show_comment()
1358 def prev_shape(self) -> None:
1359 """
1360 Decreases the index of the shape being plotted
1361 """
1362 self.current_shape -= 1
1363 if self.current_shape < 0:
1364 self.current_shape = self.shapes.size - 1
1365 self.compute_displacements()
1366 self.phase = 0.0
1367 self.last_time = time.time()
1368 self.update_shape()
1369 self.show_comment()
1371 def reset_shape(self) -> None:
1372 """
1373 Resets a shape if the value of the current shape is changed
1374 """
1375 self.compute_displacements()
1376 self.phase = 0.0
1377 self.last_time = time.time()
1378 self.update_shape()
1379 self.show_comment()
1381 def select_shape(self) -> None:
1382 from .sdynpy_shape import ShapeCommentTable
1383 self.shape_comment_table = ShapeCommentTable(self.shapes, self)
1385 def show_comment(self):
1386 """
1387 Shows or hides the shape comment
1388 """
1389 if self.plot_comments_action.isChecked():
1390 self.text.SetText(2, '\n'.join([
1391 '{:} {:}'.format(self.shape_name, self.current_shape + 1),
1392 ' Frequency: {:0.2f}'.format(self.shapes[self.current_shape].frequency)] +
1393 ([' Damping: {:.2f}%'.format(self.shapes[self.current_shape].damping * 100)] if self.show_damping else []) +
1394 [' ' + self.shapes[self.current_shape].comment1,
1395 ' ' + self.shapes[self.current_shape].comment2,
1396 ' ' + self.shapes[self.current_shape].comment3,
1397 ' ' + self.shapes[self.current_shape].comment4,
1398 ' ' + self.shapes[self.current_shape].comment5]))
1399 else:
1400 self.text.SetText(2, '')
1402 def _close(self):
1403 self.stop_animation()
1404 super(ShapePlotter, self)._close()
1407class DeflectionShapePlotter(ShapePlotter):
1408 """Class used to plot animated deflection shapes from spectra"""
1410 def __init__(self, geometry, deflection_shape_data, plot_kwargs,
1411 background_plotter_kwargs={'auto_update': 0.01},
1412 undeformed_opacity=0.25, starting_scale=1.0,
1413 deformed_opacity=1.0, num_curves=50):
1414 """
1415 Create a DeflectionShapePlotter object to plot shapes
1417 Parameters
1418 ----------
1419 geometry : Geometry
1420 Geometry on which the shapes will be plotted
1421 deflection_shape_data : NDDataArray
1422 Data array containing the deflection shapes to plot
1423 plot_kwargs : dict
1424 Keyword arguments passed to the Geometry.plot function
1425 background_plotter_kwargs : dict, optional
1426 Keyword arguments passed to the BackgroundPlotter constructor.
1427 The default is {'auto_update':0.01}.
1428 undeformed_opacity : float, optional
1429 Opacity of the undeformed geometry. Set to zero for no undeformed
1430 geometry. The default is 0.25.
1431 starting_scale : float, optional
1432 Starting scale of the shapes on the plot. The default is 1.0.
1433 num_curves : int, optional
1434 Maximum number of curves to plot on the frequency selector.
1435 Default is 50
1436 """
1437 # Convert function data into shapes
1438 from .sdynpy_shape import shape_array
1439 flattened_data = deflection_shape_data.flatten()
1440 response_coordinates = flattened_data.response_coordinate
1441 if not np.unique(response_coordinates).size == response_coordinates.size:
1442 raise ValueError('all response coordinates for operating data must be unique')
1443 shapes = shape_array(response_coordinates, flattened_data.ordinate.T,
1444 flattened_data[0].abscissa)
1445 super().__init__(geometry, shapes, plot_kwargs, background_plotter_kwargs,
1446 undeformed_opacity, starting_scale, deformed_opacity, shape_name='Shape',
1447 show_damping=False)
1449 # Get the layout so we can add a new widget to it
1450 layout = self.frame.layout()
1451 self.frequency_selector_plot = pg.PlotWidget()
1452 layout.addWidget(self.frequency_selector_plot)
1453 sp = self.frequency_selector_plot.sizePolicy()
1454 sp.setVerticalPolicy(QSizePolicy.Minimum)
1455 self.frequency_selector_plot.setSizePolicy(sp)
1456 self.frequency_selector_plot.sizeHint = lambda: QtCore.QSize(1006, 120)
1458 # Now we need to plot the data in the plot
1459 if num_curves < flattened_data.size:
1460 indices = np.linspace(0, flattened_data.size - 1, num_curves).astype(int)
1461 else:
1462 indices = np.arange(flattened_data.size)
1463 for data in flattened_data[indices]:
1464 self.frequency_selector_plot.plot(data.abscissa, np.abs(data.ordinate))
1465 # Plot the envelope
1466 max_envelope = np.linalg.svd(
1467 flattened_data.ordinate.T[..., np.newaxis], False, False, False).squeeze()
1468 self.frequency_selector_plot.plot(
1469 flattened_data[0].abscissa, max_envelope, pen=pg.mkPen(color='k', width=2))
1471 # Add an infinite line for a cursor
1472 self.frequency_selector = pg.InfiniteLine(
1473 0, movable=True, bounds=[flattened_data.abscissa.min(),
1474 flattened_data.abscissa.max()])
1475 self.frequency_selector_plot.addItem(self.frequency_selector)
1476 self.frequency_selector.sigPositionChanged.connect(self.modify_abscissa)
1478 self.frequency_selector_plot.setLogMode(False, True)
1479 self.frequency_selector_plot.setRange(
1480 yRange=(np.log10(max_envelope.min()), np.log10(max_envelope.max())))
1482 def modify_abscissa(self):
1483 frequency = self.frequency_selector.value()
1484 self.current_shape = np.argmin(np.abs(self.shapes.frequency - frequency))
1485 self.compute_displacements()
1486 self.update_shape()
1487 self.show_comment()
1489 def save_multiline_animation(self, filename=None, frame_rate=20,
1490 phase_change_per_frame=0.087,
1491 frequency_change_per_frame=0.5,
1492 start_phase=0.0,
1493 start_frequency=None, end_frequency=None,
1494 individual_images=False):
1495 self.stop_animation()
1496 imgs = []
1497 if start_frequency is None:
1498 start_frequency = self.shapes.frequency[0]
1499 if end_frequency is None:
1500 end_frequency = self.shapes.frequency[-1]
1501 frequency_frames = np.arange(start_frequency, end_frequency, frequency_change_per_frame)
1502 phase_frames = start_phase + (np.arange(frequency_frames.size) * phase_change_per_frame)
1503 for frequency, phase in zip(frequency_frames, phase_frames):
1504 self.current_shape = np.argmin(np.abs(self.shapes.frequency - frequency))
1505 self.compute_displacements()
1506 self.update_shape_mode(phase)
1507 self.show_comment()
1508 # Force GUI Update to ensure the OS doesn't delay events
1509 QApplication.processEvents()
1510 imgs.append(self.screenshot())
1511 if filename is None:
1512 return imgs
1513 else:
1514 if individual_images:
1515 for i, img in enumerate(imgs):
1516 fname = os.path.splitext(filename)
1517 fname = fname[0] + '-{:d}'.format(i) + fname[1]
1518 Image.fromarray(img).save(fname)
1519 else:
1520 imgs = [Image.fromarray(img).convert(
1521 'P', palette=Image.ADAPTIVE, colors=256) for img in imgs]
1522 imgs[0].save(fp=filename, format='GIF', append_images=imgs[1:],
1523 save_all=True, duration=int(1 / frame_rate * 1000), loop=0)
1526def split_list(seq, value):
1527 """
1528 Splits a list by a value that exists in the list
1530 Parameters
1531 ----------
1532 seq : iterable
1533 The sequence to split up
1534 value :
1535 Value within the list to use to split up the list
1537 Yields
1538 ------
1539 group : list
1540 Sublist of seq separated into groups.
1542 """
1543 group = []
1544 for val in seq:
1545 if val != value:
1546 group.append(val)
1547 elif group:
1548 yield group
1549 group = []
1550 yield group
1553def global_coord(coordinate_system, local_coord):
1554 '''Compute global coordinates from local coordinates
1556 Compute the global coordinates of a node from its local coordinates in the
1557 corresponding coordinate system.
1559 Parameters
1560 ----------
1562 coordinate_system : CoordinateSystemArray
1563 An array of coordinate systems of shape [...]
1564 local_coord : ndarray
1565 An array of coordinates of shape [...,3]
1567 Returns
1568 ------
1569 global_coord : ndarray
1570 An array of coordinates of shape [...,3]
1571 '''
1572 # Make sure we don't change local_coord
1573 local_coord = np.array(copy.deepcopy(local_coord))
1575 # First convert to cartesian if necessary
1576 cylindrical_css = coordinate_system.cs_type == 1
1577 cyl_coord = local_coord[cylindrical_css]
1578 cyl_coord[..., 0], cyl_coord[..., 1] = (cyl_coord[..., 0] * np.cos(cyl_coord[..., 1] * np.pi / 180),
1579 cyl_coord[..., 0] * np.sin(cyl_coord[..., 1] * np.pi / 180))
1580 local_coord[cylindrical_css] = cyl_coord
1582 spherical_css = coordinate_system.cs_type == 2
1583 sph_coord = local_coord[spherical_css]
1584 sph_coord[..., 0], sph_coord[..., 1], sph_coord[..., 2] = (
1585 sph_coord[..., 0] * np.sin(sph_coord[..., 1] * np.pi / 180) *
1586 np.cos(sph_coord[..., 2] * np.pi / 180),
1587 sph_coord[..., 0] * np.sin(sph_coord[..., 1] * np.pi / 180) *
1588 np.sin(sph_coord[..., 2] * np.pi / 180),
1589 sph_coord[..., 0] * np.cos(sph_coord[..., 1] * np.pi / 180))
1590 local_coord[spherical_css] = sph_coord
1591 # Convert to global cartesian from local cartesian
1592 return np.einsum('...ij,...i->...j', coordinate_system.matrix[..., :3, :], local_coord) + coordinate_system.matrix[..., 3, :]
1595def local_coord(coordinate_system, global_coord):
1596 """Compute local coordinates from global coordinates
1598 Compute the local coordinates of a node in the
1599 corresponding coordinate system from its global coordinates.
1601 Parameters
1602 ----------
1604 coordinate_system : CoordinateSystemArray
1605 An array of coordinate systems of shape [...]
1606 global_coord : ndarray
1607 An array of coordinates of shape [...,3]
1609 Returns
1610 ------
1611 local_coord : ndarray
1612 An array of coordinates of shape [...,3]
1613 """
1614 # Make sure we don't change global_coord
1615 global_coord = np.array(copy.deepcopy(global_coord))
1616 # Make sure we're working with a numpy array
1617 # Convert to the local cartesian
1618 local_coordinate = np.einsum(
1619 '...ji,...i->...j', coordinate_system.matrix[..., :3, :], global_coord - coordinate_system.matrix[..., 3, :])
1621 # Convert to cylindrical or spherical if necessary.
1622 cylindrical_css = coordinate_system.cs_type == 1
1623 cyl_coord = local_coordinate[cylindrical_css]
1624 cyl_coord[..., 0], cyl_coord[..., 1] = (np.sqrt(cyl_coord[..., 0]**2 + cyl_coord[..., 1]**2),
1625 np.arctan2(cyl_coord[..., 1], cyl_coord[..., 0]) * 180 / np.pi)
1626 local_coordinate[cylindrical_css] = cyl_coord
1628 spherical_css = coordinate_system.cs_type == 2
1629 sph_coord = local_coordinate[spherical_css]
1630 sph_coord[..., 0], sph_coord[..., 1] = (
1631 np.sqrt(sph_coord[..., 0]**2 + sph_coord[..., 1]**2 + sph_coord[..., 2]**2),
1632 np.arctan2(sph_coord[..., 1], sph_coord[..., 0]) * 180 / np.pi)
1633 # At this point local_coordinate[2,:] is still z and local_coordinate[0,:] is
1634 # now r instead of x
1635 sph_coord[..., 2] = np.arccos(sph_coord[..., 2] / sph_coord[..., 0]) * 180 / np.pi
1636 # Flip the 1 and 2 rows to be consistent with IMAT
1637 sph_coord = sph_coord[..., [0, 2, 1]]
1638 local_coordinate[spherical_css] = sph_coord
1639 return local_coordinate
1642def global_deflection(coordinate_system, local_deflection, global_point=None):
1643 """
1644 Compute deflections in the global coordinate system
1646 Parameters
1647 ----------
1648 coordinate_system : CoordinateSystemArray
1649 An array of coordinate systems of shape [...]
1650 local_deflection : ndarray
1651 An array of deflections of shape [...,3]
1652 global_point : np.ndarray, optional
1653 An array of coordinates in the global coordinate system of shape
1654 [...,3]. Must be specified if there are spherical or cylindrical
1655 coordinate systems. The default is None.
1657 Raises
1658 ------
1659 ValueError
1660 If cylindrical or spherical coordinate systems exist and global_point
1661 is not specified
1663 Returns
1664 -------
1665 global_deflection : ndarray
1666 An array of coordinates of shape [...,3] specifying the deflection
1667 directions
1669 """
1670 # Make sure we don't change local_deflection or point
1671 local_deflection = np.array(copy.deepcopy(local_deflection))
1672 if global_point is not None:
1673 global_point = np.array(copy.deepcopy(global_point))
1674 # Convert cylindrical and spherical to a local cartesian AT THE POINT
1675 # OF INTEREST (i.e. rotated in the coordinate system)
1676 local_rotation = np.zeros(coordinate_system.shape + (3, 3))
1677 for key, cs in coordinate_system.ndenumerate():
1678 if cs.cs_type == 0:
1679 this_local_rotation = np.eye(3)
1680 elif cs.cs_type == 1:
1681 if global_point is None:
1682 raise ValueError('point must be specified for Cylindrical coordinate systems')
1683 local_point = local_coord(cs, global_point[key])
1684 this_local_rotation = R(2, local_point[1], degrees=True)
1685 elif cs.cs_type == 2:
1686 if global_point is None:
1687 raise ValueError('point must be specified for Spherical coordinate systems')
1688 local_point = local_coord(cs, global_point[key])
1689 this_local_rotation = R(
1690 2, local_point[2], degrees=True) @ R(1, local_point[1], degrees=True)
1691 else:
1692 raise ValueError('invalid coordinate system type {:}'.format(cs.cs_type))
1693 local_rotation[key] = this_local_rotation
1694 return np.einsum('...ij,...j->...i', np.einsum('...ij,...ik->...jk', coordinate_system.matrix[..., :3, :], local_rotation), local_deflection)
1697class CoordinateSystemArray(SdynpyArray):
1698 """Coordinate System Information specifying position and directions.
1700 Use the coordinate_system_array helper function to create the array.
1701 """
1702 data_dtype = [('id', 'uint64'), ('name', '<U40'), ('color', 'uint16'),
1703 ('cs_type', 'uint16'), ('matrix', 'float64', (4, 3))]
1705 def __new__(subtype, shape, buffer=None, offset=0,
1706 strides=None, order=None):
1707 # Create the ndarray instance of our type, given the usual
1708 # ndarray input arguments. This will call the standard
1709 # ndarray constructor, but return an object of our type.
1710 # It also triggers a call to InfoArray.__array_finalize__
1711 obj = super(CoordinateSystemArray, subtype).__new__(subtype, shape,
1712 CoordinateSystemArray.data_dtype, buffer, offset, strides, order)
1713 # Finally, we must return the newly created object:
1714 return obj
1716 def __repr__(self):
1717 string_out = '{:>8s}, {:>6s}, {:>20s}, {:>5s}, {:>10s}\n'.format(
1718 'Index', 'ID', 'Name', 'Color', 'Type')
1719 if self.size == 0:
1720 string_out += '----------- Empty -------------\n'
1721 for i, (key, val) in enumerate(self.ndenumerate()):
1722 if i >= MAX_NUMBER_REPR:
1723 string_out += ' .\n .\n .\n'
1724 break
1725 string_out += '{:>8}, {:>6}, {:>20}, {:>5}, {:>10}\n'.format(
1726 str(key), val.id, val.name, val.color, _cs_type_strings[val.cs_type[()]])
1727 return string_out
1729 def __call__(self, index_by_id):
1730 """
1731 Select coordinate system by id rather than by index
1733 Parameters
1734 ----------
1735 index_by_id : int or np.ndarray
1736 ID number(s) to use to select coordinate systems.
1738 Raises
1739 ------
1740 ValueError
1741 If specified ID(s) not found in array.
1743 Returns
1744 -------
1745 output : CoordinateSystemArray
1746 Subset of CoordinateSystemArray with the specified IDs.
1748 """
1749 index_dict = {node.id[()]: index for index, node in self.ndenumerate()}
1750 ids = np.array(index_by_id)
1751 output = np.empty(ids.shape, dtype=self.dtype).view(self.__class__)
1752 for key, val in np.ndenumerate(ids):
1753 try:
1754 index = index_dict[val]
1755 except KeyError:
1756 raise ValueError('ID {:} not found in array'.format(val))
1757 output[key] = self[index]
1758 return output
1760 @staticmethod
1761 def from_unv(unv_data_dict, combine=True):
1762 """
1763 Load CoordinateSystemArrays from universal file data from read_unv
1765 Parameters
1766 ----------
1767 unv_data_dict : dict
1768 Dictionary containing data from read_unv
1769 combine : bool, optional
1770 If True, return as a single CoordinateSytemArray
1772 Returns
1773 -------
1774 output_arrays : CoordinateSystemArray
1775 Coordinate Systems read from unv
1777 """
1778 try:
1779 datasets = unv_data_dict[2420]
1780 except KeyError:
1781 if combine:
1782 return CoordinateSystemArray(0)
1783 else:
1784 return []
1785 output_arrays = []
1786 for dataset in datasets:
1787 output_arrays.append(coordinate_system_array(dataset.cs_labels,
1788 dataset.cs_names,
1789 dataset.cs_colors,
1790 dataset.cs_types,
1791 dataset.cs_matrices))
1792 if combine:
1793 output_arrays = np.concatenate(output_arrays)
1794 return output_arrays
1797def coordinate_system_array(id=1, name='', color=1, cs_type=0, matrix=np.concatenate((np.eye(3), np.zeros((1, 3))), axis=0),
1798 structured_array=None):
1799 """
1800 Creates an array that specifies coordinate systems in the geometry
1802 Creates an array of coordinate systems that specify directions of sensors
1803 in a test or analysis. Coordinate system arrays can be created using a
1804 numpy structured array or individual arrays for each attribute.
1805 Multidimensional arrays can be used.
1807 Parameters
1808 ----------
1809 id : ndarray
1810 Integer array corresponding to the id of the coordinate systems. Input
1811 will be cast to an integer (i.e. 2.0 -> 2, 1.9 -> 1)
1812 name : ndarray
1813 name of the coordinate systems
1814 color : ndarray
1815 color of the coordinate systems as integers
1816 cs_type : ndarray
1817 Coordinate system types (0 = Cartesian, 1 = Polar, 2 = Spherical)
1818 matrix : ndarray
1819 Coordinate system transformation matrix with shape [... x 4 x 3]
1820 structured_array : ndarray (structured)
1821 Alternatively to the individual attributes, a single numpy structured
1822 array can be passed, which should have the same name as the inputs to
1823 the function listed above.
1825 Returns
1826 -------
1827 coordinate_system_array : CoordinateSystemArray
1828 """
1829 if structured_array is not None:
1830 try:
1831 id = structured_array['id']
1832 name = structured_array['name']
1833 color = structured_array['color']
1834 cs_type = structured_array['cs_type']
1835 matrix = structured_array['matrix']
1836 except (ValueError, TypeError):
1837 raise ValueError(
1838 'structured_array must be numpy.ndarray with dtype names "id", "name", "color", "cs_type" and "matrix"')
1839 else:
1840 id = np.array(id)
1841 name = np.array(name)
1842 color = np.array(color)
1843 cs_type = np.array(cs_type)
1844 matrix = np.array(matrix)
1846 # Don't check shapes because we want them to be broadcastable, but id have to
1847 # be unique so we will use that for the shape
1848 coord_sys_array = CoordinateSystemArray(id.shape)
1849 coord_sys_array.id = id
1850 coord_sys_array.name = name
1851 coord_sys_array.color = color
1852 coord_sys_array.cs_type = cs_type
1853 coord_sys_array.matrix = matrix
1855 return coord_sys_array
1858class ElementArray(SdynpyArray):
1859 """Element information array
1861 Use the element_array helper function to create the array.
1862 """
1864 data_dtype = [('id', 'uint64'), ('type', 'uint8'), ('color', 'uint16'),
1865 ('connectivity', 'object')]
1867 def __new__(subtype, shape, buffer=None, offset=0,
1868 strides=None, order=None):
1869 # Create the ndarray instance of our type, given the usual
1870 # ndarray input arguments. This will call the standard
1871 # ndarray constructor, but return an object of our type.
1872 # It also triggers a call to InfoArray.__array_finalize__
1874 obj = super(ElementArray, subtype).__new__(subtype, shape,
1875 ElementArray.data_dtype, buffer, offset, strides, order)
1876 # Finally, we must return the newly created object:
1877 return obj
1879 def __repr__(self):
1880 string_out = '{:>8s}, {:>6s}, {:>4s}, {:>5s}, {:>7s}\n'.format(
1881 'Index', 'ID', 'Type', 'Color', '# Nodes')
1882 if self.size == 0:
1883 string_out += '----------- Empty -------------\n'
1884 for i, (key, val) in enumerate(self.ndenumerate()):
1885 if i >= MAX_NUMBER_REPR:
1886 string_out += ' .\n .\n .\n'
1887 break
1888 string_out += '{:>8}, {:>6}, {:>4}, {:>5}, {:>7}\n'.format(
1889 str(key), val.id, val.type, val.color, len(val.connectivity))
1890 return string_out
1892 def __call__(self, index_by_id):
1893 """
1894 Select elements by id rather than by index
1896 Parameters
1897 ----------
1898 index_by_id : int or np.ndarray
1899 ID number(s) to use to select elements.
1901 Raises
1902 ------
1903 ValueError
1904 If specified ID(s) not found in array.
1906 Returns
1907 -------
1908 output : ElementArray
1909 Subset of ElementArray with the specified IDs.
1911 """
1912 index_dict = {node.id[()]: index for index, node in self.ndenumerate()}
1913 ids = np.array(index_by_id)
1914 output = np.empty(ids.shape, dtype=self.dtype).view(self.__class__)
1915 for key, val in np.ndenumerate(ids):
1916 try:
1917 index = index_dict[val]
1918 except KeyError:
1919 raise ValueError('ID {:} not found in array'.format(val))
1920 output[key] = self[index]
1921 return output
1923 def reduce(self, node_list):
1924 """
1925 Keep only elements that have all nodes contained in node_list
1927 Parameters
1928 ----------
1929 node_list : iterable
1930 Iterable containing nodes to keep.
1932 Returns
1933 -------
1934 ElementArray
1935 ElementArray containing only elements with all nodes in node_list.
1937 """
1938 output_list = []
1939 for key, element in self.ndenumerate():
1940 if np.all(np.isin(element.connectivity, node_list)):
1941 output_list.append(element)
1942 return np.array(output_list, self.dtype).view(ElementArray)
1944 def remove(self, element_list):
1945 """
1946 Removes elements with id numbers in traceline_list
1948 Parameters
1949 ----------
1950 traceline_list : interable
1951 Iterable containing tracelines to discard.
1953 Returns
1954 -------
1955 TracelineArray
1956 TracelineArray containing tracelines with id numbers not it
1957 traceline_list.
1959 """
1960 output_list = []
1961 for key,element in self.ndenumerate():
1962 if element.id not in element_list:
1963 output_list.append(element)
1964 return np.array(output_list, self.dtype).view(ElementArray)
1966 @staticmethod
1967 def from_unv(unv_data_dict, combine=True):
1968 """
1969 Load ElementArrays from universal file data from read_unv
1971 Parameters
1972 ----------
1973 unv_data_dict : dict
1974 Dictionary containing data from read_unv
1975 combine : bool, optional
1976 If True, return as a single ElementArray
1978 Returns
1979 -------
1980 output_arrays : ElementArray
1981 Elements read from unv
1983 """
1984 try:
1985 datasets = unv_data_dict[2412]
1986 except KeyError:
1987 if combine:
1988 return ElementArray(0)
1989 else:
1990 return []
1991 output_arrays = []
1992 for dataset in datasets:
1993 output_arrays.append(element_array(dataset.element_labels,
1994 dataset.fe_descriptor_ids,
1995 dataset.colors,
1996 dataset.connectivities))
1997 if combine:
1998 output_arrays = np.concatenate(output_arrays)
1999 return output_arrays
2002def element_array(id=1, type=None, color=1,
2003 connectivity=None, structured_array=None):
2004 """
2005 Creates an array that specifies elements in the geometry
2007 Creates an array of elements that specify connectivity of sensors
2008 in a test or analysis. Element arrays can be created using a
2009 numpy structured array or individual arrays for each attribute.
2010 Multidimensional arrays can be used.
2012 Parameters
2013 ----------
2014 id : ndarray
2015 Integer array corresponding to the id of the node. Input
2016 will be cast to an integer (i.e. 2.0 -> 2, 1.9 -> 1)
2017 type : ndarray
2018 Element types. See notes for type mapping
2019 color : ndarray
2020 color of the elements as integers
2021 connectivity : ndarray
2022 An object array of iterables defining the connectivity for each element.
2023 structured_array : ndarray (structured)
2024 Alternatively to the individual attributes, a single numpy structured
2025 array can be passed, which should have the same name as the inputs to
2026 the function listed above.
2029 Returns
2030 -------
2031 element_array : ElementArray
2034 Notes
2035 -----
2036 Here is a list of element types
2038 11 : 'Rod',
2039 21 : 'Linear beam',
2040 22 : 'Tapered beam',
2041 23 : 'Curved beam',
2042 24 : 'Parabolic beam',
2043 31 : 'Straight pipe',
2044 32 : 'Curved pipe',
2045 41 : 'Plane Stress Linear Triangle',
2046 42 : 'Plane Stress Parabolic Triangle',
2047 43 : 'Plane Stress Cubic Triangle',
2048 44 : 'Plane Stress Linear Quadrilateral',
2049 45 : 'Plane Stress Parabolic Quadrilateral',
2050 46 : 'Plane Strain Cubic Quadrilateral',
2051 51 : 'Plane Strain Linear Triangle',
2052 52 : 'Plane Strain Parabolic Triangle',
2053 53 : 'Plane Strain Cubic Triangle',
2054 54 : 'Plane Strain Linear Quadrilateral',
2055 55 : 'Plane Strain Parabolic Quadrilateral',
2056 56 : 'Plane Strain Cubic Quadrilateral',
2057 61 : 'Plate Linear Triangle',
2058 62 : 'Plate Parabolic Triangle',
2059 63 : 'Plate Cubic Triangle',
2060 64 : 'Plate Linear Quadrilateral',
2061 65 : 'Plate Parabolic Quadrilateral',
2062 66 : 'Plate Cubic Quadrilateral',
2063 71 : 'Membrane Linear Quadrilateral',
2064 72 : 'Membrane Parabolic Triangle',
2065 73 : 'Membrane Cubic Triangle',
2066 74 : 'Membrane Linear Triangle',
2067 75 : 'Membrane Parabolic Quadrilateral',
2068 76 : 'Membrane Cubic Quadrilateral',
2069 81 : 'Axisymetric Solid Linear Triangle',
2070 82 : 'Axisymetric Solid Parabolic Triangle',
2071 84 : 'Axisymetric Solid Linear Quadrilateral',
2072 85 : 'Axisymetric Solid Parabolic Quadrilateral',
2073 91 : 'Thin Shell Linear Triangle',
2074 92 : 'Thin Shell Parabolic Triangle',
2075 93 : 'Thin Shell Cubic Triangle',
2076 94 : 'Thin Shell Linear Quadrilateral',
2077 95 : 'Thin Shell Parabolic Quadrilateral',
2078 96 : 'Thin Shell Cubic Quadrilateral',
2079 101: 'Thick Shell Linear Wedge',
2080 102: 'Thick Shell Parabolic Wedge',
2081 103: 'Thick Shell Cubic Wedge',
2082 104: 'Thick Shell Linear Brick',
2083 105: 'Thick Shell Parabolic Brick',
2084 106: 'Thick Shell Cubic Brick',
2085 111: 'Solid Linear Tetrahedron',
2086 112: 'Solid Linear Wedge',
2087 113: 'Solid Parabolic Wedge',
2088 114: 'Solid Cubic Wedge',
2089 115: 'Solid Linear Brick',
2090 116: 'Solid Parabolic Brick',
2091 117: 'Solid Cubic Brick',
2092 118: 'Solid Parabolic Tetrahedron',
2093 121: 'Rigid Bar',
2094 122: 'Rigid Element',
2095 136: 'Node To Node Translational Spring',
2096 137: 'Node To Node Rotational Spring',
2097 138: 'Node To Ground Translational Spring',
2098 139: 'Node To Ground Rotational Spring',
2099 141: 'Node To Node Damper',
2100 142: 'Node To Gound Damper',
2101 151: 'Node To Node Gap',
2102 152: 'Node To Ground Gap',
2103 161: 'Lumped Mass',
2104 171: 'Axisymetric Linear Shell',
2105 172: 'Axisymetric Parabolic Shell',
2106 181: 'Constraint',
2107 191: 'Plastic Cold Runner',
2108 192: 'Plastic Hot Runner',
2109 193: 'Plastic Water Line',
2110 194: 'Plastic Fountain',
2111 195: 'Plastic Baffle',
2112 196: 'Plastic Rod Heater',
2113 201: 'Linear node-to-node interface',
2114 202: 'Linear edge-to-edge interface',
2115 203: 'Parabolic edge-to-edge interface',
2116 204: 'Linear face-to-face interface',
2117 208: 'Parabolic face-to-face interface',
2118 212: 'Linear axisymmetric interface',
2119 213: 'Parabolic axisymmetric interface',
2120 221: 'Linear rigid surface',
2121 222: 'Parabolic rigid surface',
2122 231: 'Axisymetric linear rigid surface',
2123 232: 'Axisymentric parabolic rigid surface'
2124 """
2125 if structured_array is not None:
2126 try:
2127 id = structured_array['id']
2128 type = structured_array['type']
2129 color = structured_array['color']
2130 connectivity = structured_array['connectivity']
2131 except (ValueError, TypeError):
2132 raise ValueError(
2133 'structured_array must be numpy.ndarray with dtype names "id", "type", "color", and "connectivity"')
2134 else:
2135 id = np.array(id)
2136 type = np.array(type)
2137 color = np.array(color)
2138 # Need to make connectivity an object array so don't do this!
2139 # connectivity = np.array(connectivity,dtype=object)
2141 # Don't check shapes because we want them to be broadcastable, but id have to
2142 # be unique so we will use that for the shape
2143 earray = ElementArray(id.shape)
2144 earray.id = id
2145 earray.type = type
2146 earray.color = color
2147 if earray.ndim == 0:
2148 # This is gross, but it works!
2149 np.ndarray.__getitem__(earray, 'connectivity')[()] = np.array(connectivity, dtype='uint32')
2150 else:
2151 connectivity = np.array(connectivity, dtype=object)
2152 for key, val in np.ndenumerate(id):
2153 earray.connectivity[key] = np.array(connectivity[key], dtype='uint32')
2154 return earray
2157class NodeArray(SdynpyArray):
2158 """Node information array
2160 Use the node_array helper function to create the array.
2161 """
2163 data_dtype = [('id', 'uint64'), ('coordinate', 'float64', (3,)), ('color', 'uint16'),
2164 ('def_cs', 'uint64'), ('disp_cs', 'uint64')]
2166 def __new__(subtype, shape, buffer=None, offset=0,
2167 strides=None, order=None):
2168 # Create the ndarray instance of our type, given the usual
2169 # ndarray input arguments. This will call the standard
2170 # ndarray constructor, but return an object of our type.
2171 # It also triggers a call to InfoArray.__array_finalize__
2173 obj = super(NodeArray, subtype).__new__(subtype, shape,
2174 NodeArray.data_dtype, buffer, offset, strides, order)
2175 # Finally, we must return the newly created object:
2176 return obj
2178 def __repr__(self):
2179 string_out = '{:>8s}, {:>6s}, {:>8s}, {:>8s}, {:>8s}, {:>5s}, {:>5s}\n'.format(
2180 'Index', 'ID', 'X', 'Y', 'Z', 'DefCS', 'DisCS')
2181 if self.size == 0:
2182 string_out += '----------- Empty -------------\n'
2183 for i, (key, val) in enumerate(self.ndenumerate()):
2184 if i >= MAX_NUMBER_REPR:
2185 string_out += ' .\n .\n .\n'
2186 break
2187 string_out += '{:>8}, {:>6}, {:>8.3f}, {:>8.3f}, {:>8.3f}, {:>5}, {:>5}\n'.format(
2188 str(key), val.id, *val.coordinate, val.def_cs, val.disp_cs)
2189 return string_out
2191 def __call__(self, index_by_id):
2192 """
2193 Select nodes by id rather than by index
2195 Parameters
2196 ----------
2197 index_by_id : int or np.ndarray
2198 ID number(s) to use to select nodes.
2200 Raises
2201 ------
2202 ValueError
2203 If specified ID(s) not found in array.
2205 Returns
2206 -------
2207 output : NodeArray
2208 Subset of NodeArray with the specified IDs.
2210 """
2211 index_dict = {node.id[()]: index for index, node in self.ndenumerate()}
2212 ids = np.array(index_by_id)
2213 output = np.empty(ids.shape, dtype=self.dtype).view(self.__class__)
2214 for key, val in np.ndenumerate(ids):
2215 try:
2216 index = index_dict[val]
2217 except KeyError:
2218 raise ValueError('ID {:} not found in array'.format(val))
2219 output[key] = self[index]
2220 return output
2222 def reduce(self, node_list):
2223 """
2224 Keep only nodes that are contained in node_list
2226 Parameters
2227 ----------
2228 node_list : iterable
2229 Iterable containing nodes to keep.
2231 Returns
2232 -------
2233 NodeArray
2234 NodeArray containing only nodes in node_list.
2236 """
2237 return self(node_list)
2239 @staticmethod
2240 def from_unv(unv_data_dict, combine=True):
2241 """
2242 Load NodeArrays from universal file data from read_unv
2244 Parameters
2245 ----------
2246 unv_data_dict : dict
2247 Dictionary containing data from read_unv
2248 combine : bool, optional
2249 If True, return as a single NodeArray
2251 Returns
2252 -------
2253 output_arrays : NodeArray
2254 Nodes read from unv
2256 """
2257 try:
2258 datasets = unv_data_dict[2411]
2259 except KeyError:
2260 if combine:
2261 return NodeArray(0)
2262 else:
2263 return []
2264 output_arrays = []
2265 for dataset in datasets:
2266 output_arrays.append(node_array(dataset.node_labels,
2267 dataset.coordinates,
2268 dataset.colors,
2269 dataset.export_coordinate_systems,
2270 dataset.displacement_coordinate_systems))
2271 if combine:
2272 output_arrays = np.concatenate(output_arrays)
2273 return output_arrays
2275 def by_position(self, position_array):
2276 """
2277 Select node by closest position
2279 Parameters
2280 ----------
2281 position_array : np.ndarray
2282 A (...,3) shape array containing positions of nodes to keep
2284 Returns
2285 -------
2286 NodeArray
2287 NodArray containing nodes that were closest to the positions in
2288 position_array.
2290 """
2291 node_differences = np.linalg.norm(
2292 self.coordinate - position_array[(Ellipsis,) + (np.newaxis,) * self.ndim + (slice(None),)], axis=-1)
2293 node_differences = node_differences.reshape(*position_array.shape[:-1], -1)
2294 node_indices = np.argmin(
2295 node_differences,
2296 axis=-1)
2297 return self.flatten()[node_indices]
2299 @staticmethod
2300 def project_to_minimum_plane(coordinates, return_3D=True):
2301 """
2302 Projects coordinates to a single best-fit plane
2304 Parameters
2305 ----------
2306 coordinates : np.ndarray
2307 A (...,3) coordinate array.
2308 return_3D : bool, optional
2309 If True, return the 3D coordinates of the projected points.
2310 Otherwise return the projected 2D coordinate. The default is True.
2312 Returns
2313 -------
2314 np.ndarray
2315 Points projected to a best-fit plane.
2317 """
2318 coords = coordinates.T.copy()
2319 mean = np.mean(coords, axis=-1, keepdims=True)
2320 coords -= mean
2321 U, S, V = np.linalg.svd(coords)
2322 coords = U.T @ coords
2323 coords[2, :] = 0
2324 if return_3D:
2325 return (U @ coords + mean).T
2326 else:
2327 return coords[:2].T
2329 def global_coordinate(self, cs_array):
2330 """
2331 Get the global coordinate of each node
2333 Parameters
2334 ----------
2335 cs_array : CoordinateSystemArray
2336 CoordinateSystemArray consisting of the local coordinate systems for
2337 each node
2339 Returns
2340 -------
2341 points : np.ndarray
2342 Coordinates of the nodes in the global coordinate system.
2344 """
2345 points = global_coord(cs_array(self.def_cs), self.coordinate)
2346 return points
2348 def triangulate(self, cs_array, projection_function=None,
2349 return_element_array=True, element_color=1,
2350 condition_threshold=None):
2351 """
2352 Creates elements for a node set using Delaunay Triangulation
2354 Parameters
2355 ----------
2356 cs_array : CoordinateSystemArray
2357 CoordinateSystemArray containing coordinate systems for each node.
2358 projection_function : function, optional
2359 Function to use to project 3D coordinates to 2D coordinates for
2360 triangulation. The default is None.
2361 return_element_array : bool, optional
2362 Returns an ElementArray if True, otherwise it simply returns the
2363 triangle simplices. The default is True.
2364 element_color : np.ndarray or int, optional
2365 Integers representing colors applied to the elements.
2366 The default is 1.
2367 condition_threshold : float, optional
2368 Condition number threshold used to remove triangles that are poorly
2369 shaped. The default is None.
2371 Returns
2372 -------
2373 ElementArray or np.ndarray
2374 ElementArray containing elements or np.ndarray containing triangle
2375 simplices.
2377 """
2378 if projection_function is None:
2379 def projection_function(coords): return self.project_to_minimum_plane(coords, False)
2380 # Compute global positions of the nodes
2381 global_positions = self.global_coordinate(cs_array)
2382 # Now compute the projection function
2383 projected_positions = projection_function(global_positions)
2384 # Now triangulate
2385 tri = Delaunay(projected_positions)
2386 simplices = tri.simplices
2387 if condition_threshold is not None:
2388 # Now remove small triangles via condition number
2389 tri_points = projected_positions[tri.simplices, :]
2390 tri_points -= np.mean(tri_points, axis=1, keepdims=True)
2391 simplices = simplices[np.linalg.cond(tri_points) < condition_threshold]
2392 if not return_element_array:
2393 return simplices
2394 else:
2395 # Create an element array
2396 connectivity = self.id[simplices]
2397 return element_array(id=np.arange(connectivity.shape[0]) + 1,
2398 type=61, color=element_color,
2399 connectivity=connectivity)
2401 def by_grid(self, grid_spacing, x_min=None, y_min=None, z_min=None,
2402 x_max=None, y_max=None, z_max=None):
2403 """
2404 Selects nodes in a grid
2406 Parameters
2407 ----------
2408 grid_spacing : float
2409 Approximate grid spacing between selected nodes
2410 x_min : float, optional
2411 Minimum X-value of the grid.
2412 The default is the minimum node x-coordinate.
2413 y_min : float, optional
2414 Minimun Y-value of the grid.
2415 The default is the minimum node y-coordinate.
2416 z_min : float, optional
2417 Minimum Z-value of the grid.
2418 The default is the minimum node z-coordinate.
2419 x_max : float, optional
2420 Maximum X-value of the grid.
2421 The default is the maximum node x-coordinate.
2422 y_max : float, optional
2423 Maximum Y-value of the grid.
2424 The default is the maximum node y-coordinate.
2425 z_max : float, optional
2426 Maximum Z-value of the grid.
2427 The default is the maximum node z-coordinate.
2429 Returns
2430 -------
2431 NodeArray
2432 NodeArray containing the closest nodes to the grid points
2433 """
2434 if x_min is None:
2435 x_min = self.coordinate[:, 0].min()
2436 if y_min is None:
2437 y_min = self.coordinate[:, 1].min()
2438 if z_min is None:
2439 z_min = self.coordinate[:, 2].min()
2440 if x_max is None:
2441 x_max = self.coordinate[:, 0].max()
2442 if y_max is None:
2443 y_max = self.coordinate[:, 1].max()
2444 if z_max is None:
2445 z_max = self.coordinate[:, 2].max()
2446 min_coords = np.array((x_min, y_min, z_min))
2447 max_coords = np.array((x_max, y_max, z_max))
2448 range_coords = max_coords - min_coords
2449 num_grids = np.ceil(range_coords / grid_spacing).astype(int)
2450 # Now we will create the grid in each dimension using linspace
2451 grid_arrays = [np.linspace(min_coord, max_coord, num_coord)
2452 for min_coord, max_coord, num_coord
2453 in zip(min_coords, max_coords, num_grids)]
2454 # We will then use meshgrid to assemble the array grid. We will flatten the
2455 # point dimension and transpose so the array has shape (n_points x 3)
2456 grid_points = np.array(np.meshgrid(*grid_arrays, indexing='ij')).reshape(3, -1).T
2457 # Now we will select nodes that are closest to the points in the grid
2458 candidate_nodes = self.by_position(grid_points)
2459 # Reduce to only nodes that are within one grid spacing of their requested point
2460 # this reduces the nodes that are selected by points far from the model
2461 candidate_nodes = candidate_nodes[
2462 np.linalg.norm(candidate_nodes.coordinate - grid_points, axis=-1) < grid_spacing]
2463 # Remove any duplicates
2464 candidate_node_ids = np.unique(candidate_nodes.id)
2465 return candidate_nodes(candidate_node_ids)
2468def node_array(id=1, coordinate=np.array((0, 0, 0)), color=1, def_cs=1, disp_cs=1,
2469 structured_array=None):
2470 """
2471 Creates an array that specifies nodes in the geometry
2473 Creates an array of nodes that specify positions of sensors
2474 in a test or analysis. Node arrays can be created using a
2475 numpy structured array or individual arrays for each attribute.
2476 Multidimensional arrays can be used.
2478 Parameters
2479 ----------
2480 id : ndarray
2481 Integer array corresponding to the id of the node. Input
2482 will be cast to an integer (i.e. 2.0 -> 2, 1.9 -> 1)
2483 coordinate : ndarray
2484 Positions of the nodes in space
2485 color : ndarray
2486 color of the nodes as integers
2487 def_cs : ndarray
2488 Coordinate system where the node's position is defined
2489 disp_cs : ndarray
2490 Coordinate system where the node's displacements are defined
2491 structured_array : ndarray (structured)
2492 Alternatively to the individual attributes, a single numpy structured
2493 array can be passed, which should have the same name as the inputs to
2494 the function listed above.
2497 Returns
2498 -------
2499 node_array : NodeArray
2500 """
2501 if structured_array is not None:
2502 try:
2503 id = structured_array['id']
2504 coordinate = structured_array['coordinate']
2505 color = structured_array['color']
2506 def_cs = structured_array['def_cs']
2507 disp_cs = structured_array['disp_cs']
2508 except (ValueError, TypeError):
2509 raise ValueError(
2510 'structured_array must be numpy.ndarray with dtype names "id", "coordinate", "color", "def_cs" and "disp_cs"')
2511 else:
2512 id = np.array(id)
2513 coordinate = np.array(coordinate)
2514 color = np.array(color)
2515 def_cs = np.array(def_cs)
2516 disp_cs = np.array(disp_cs)
2518 # Don't check shapes because we want them to be broadcastable, but id have to
2519 # be unique so we will use that for the shape
2520 narray = NodeArray(id.shape)
2521 narray.id = id
2522 narray.coordinate = coordinate
2523 narray.color = color
2524 narray.def_cs = def_cs
2525 narray.disp_cs = disp_cs
2527 return narray
2530class TracelineArray(SdynpyArray):
2531 """Traceline information array
2533 Use the traceline_array helper function to create the array.
2534 """
2535 data_dtype = [('id', 'uint64'), ('color', 'uint16'), ('description', '<U40'),
2536 ('connectivity', 'object')]
2538 def __new__(subtype, shape, buffer=None, offset=0,
2539 strides=None, order=None):
2540 # Create the ndarray instance of our type, given the usual
2541 # ndarray input arguments. This will call the standard
2542 # ndarray constructor, but return an object of our type.
2543 # It also triggers a call to InfoArray.__array_finalize__
2545 obj = super(TracelineArray, subtype).__new__(subtype, shape,
2546 TracelineArray.data_dtype, buffer, offset, strides, order)
2547 # Finally, we must return the newly created object:
2548 return obj
2550 def __repr__(self):
2551 string_out = '{:>8s}, {:>6s}, {:>20s}, {:>5s}, {:>7s}\n'.format(
2552 'Index', 'ID', 'Description', 'Color', '# Nodes')
2553 if self.size == 0:
2554 string_out += '----------- Empty -------------\n'
2555 for i, (key, val) in enumerate(self.ndenumerate()):
2556 if i >= MAX_NUMBER_REPR:
2557 string_out += ' .\n .\n .\n'
2558 break
2559 string_out += '{:>8s}, {:>6d}, {:>20s}, {:>5d}, {:>7d}\n'.format(
2560 str(key), val.id, val.description, val.color, len(val.connectivity))
2561 return string_out
2563 def __call__(self, index_by_id):
2564 """
2565 Select nodes by id rather than by index
2567 Parameters
2568 ----------
2569 index_by_id : int or np.ndarray
2570 ID number(s) to use to select nodes.
2572 Raises
2573 ------
2574 ValueError
2575 If specified ID(s) not found in array.
2577 Returns
2578 -------
2579 output : TracelineArray
2580 Subset of TracelineArray with the specified IDs.
2582 """
2583 index_dict = {node.id[()]: index for index, node in self.ndenumerate()}
2584 ids = np.array(index_by_id)
2585 output = np.empty(ids.shape, dtype=self.dtype).view(self.__class__)
2586 for key, val in np.ndenumerate(ids):
2587 try:
2588 index = index_dict[val]
2589 except KeyError:
2590 raise ValueError('ID {:} not found in array'.format(val))
2591 output[key] = self[index]
2592 return output
2594 def reduce(self, node_list):
2595 """
2596 Keep only tracelines fully contain nodes in node_list
2598 Parameters
2599 ----------
2600 node_list : iterable
2601 Iterable containing nodes to keep.
2603 Returns
2604 -------
2605 TracelineArray
2606 TracelineArray containing lines that contain only nodes in node_list.
2608 """
2609 output_list = []
2610 for key, traceline in self.ndenumerate():
2611 if np.all(np.isin(traceline.connectivity, node_list) | (traceline.connectivity == 0)):
2612 output_list.append(traceline)
2613 return np.array(output_list, self.dtype).view(TracelineArray)
2615 def remove(self, traceline_list):
2616 """
2617 Removes tracelines with id numbers in traceline_list
2619 Parameters
2620 ----------
2621 traceline_list : interable
2622 Iterable containing tracelines to discard.
2624 Returns
2625 -------
2626 TracelineArray
2627 TracelineArray containing tracelines with id numbers not it
2628 traceline_list.
2630 """
2631 output_list = []
2632 for key,traceline in self.ndenumerate():
2633 if traceline.id not in traceline_list:
2634 output_list.append(traceline)
2635 return np.array(output_list, self.dtype).view(TracelineArray)
2637 @staticmethod
2638 def from_unv(unv_data_dict, combine=True):
2639 """
2640 Load TracelineArrays from universal file data from read_unv
2642 Parameters
2643 ----------
2644 unv_data_dict : dict
2645 Dictionary containing data from read_unv
2646 combine : bool, optional
2647 If True, return as a single TracelineArray object
2649 Returns
2650 -------
2651 output_arrays : TracelineArray
2652 Tracelines read from unv
2654 """
2655 try:
2656 datasets = unv_data_dict[82]
2657 except KeyError:
2658 if combine:
2659 return TracelineArray(0)
2660 else:
2661 return []
2662 output_arrays = []
2663 for dataset in datasets:
2664 output_arrays.append(traceline_array(dataset.traceline_number,
2665 dataset.identification,
2666 dataset.color,
2667 dataset.nodes)[np.newaxis])
2668 if combine:
2669 output_arrays = np.concatenate(output_arrays)
2670 return output_arrays
2673def traceline_array(id=1, description='', color=1,
2674 connectivity=None, structured_array=None):
2675 """
2676 Creates an array that specifies tracelines in the geometry
2678 Creates an array of tracelines that specify connectivity of sensors
2679 in a test or analysis. Traceline arrays can be created using a
2680 numpy structured array or individual arrays for each attribute.
2681 Multidimensional arrays can be used.
2683 Parameters
2684 ----------
2685 id : ndarray
2686 Integer array corresponding to the id of the node. Input
2687 will be cast to an integer (i.e. 2.0 -> 2, 1.9 -> 1)
2688 Description : ndarray
2689 Traceline name
2690 color : ndarray
2691 color of the elements as integers
2692 connectivity : ndarray
2693 An object array of iterables defining the connectivity for each traceline.
2694 structured_array : ndarray (structured)
2695 Alternatively to the individual attributes, a single numpy structured
2696 array can be passed, which should have the same name as the inputs to
2697 the function listed above.
2700 Returns
2701 -------
2702 traceline_array : TracelineArray
2704 """
2705 if structured_array is not None:
2706 try:
2707 id = structured_array['id']
2708 description = structured_array['description']
2709 color = structured_array['color']
2710 connectivity = structured_array['connectivity']
2711 except (ValueError, TypeError):
2712 raise ValueError(
2713 'structured_array must be numpy.ndarray with dtype names "id", "color", "description", and "connectivity"')
2714 else:
2715 id = np.array(id)
2716 description = np.array(description)
2717 color = np.array(color)
2718 # Don't do this, need to keep as an object dtype
2719 # connectivity = np.array(connectivity)
2721 # Don't check shapes because we want them to be broadcastable, but id have to
2722 # be unique so we will use that for the shape
2723 tlarray = TracelineArray(id.shape)
2724 tlarray.id = id
2725 tlarray.description = description
2726 tlarray.color = color
2727 if tlarray.ndim == 0:
2728 np.ndarray.__getitem__(tlarray, 'connectivity')[()] = np.array(connectivity, dtype='uint32')
2729 else:
2730 connectivity = np.array(connectivity, dtype=object)
2731 for key, val in np.ndenumerate(id):
2732 tlarray.connectivity[key] = np.array(connectivity[key], dtype='uint32')
2734 return tlarray
2737from ..fem.sdynpy_exodus import Exodus, ExodusInMemory # noqa: E402
2738from .sdynpy_shape import shape_array # noqa: E402
2741class Geometry:
2742 """Container for nodes, coordinate systems, tracelines, and elements
2744 Geometry is the class that is most useful for working with the positioning
2745 and spatial visualization of test or analysis data. It contains functions
2746 to plot and animate 3D geometry"""
2748 def __init__(self, node=NodeArray((0,)), coordinate_system=CoordinateSystemArray((0,)),
2749 traceline=TracelineArray((0,)), element=ElementArray((0,))):
2750 """
2751 Initialize a geometry object with nodes, coordinate systems, tracelines,
2752 and elements.
2754 All input arguments will be flattened when passed to the Geometry, as
2755 the geometry does not support multi-dimensional object arrays.
2757 Parameters
2758 ----------
2759 node : NodeArray, optional
2760 The set of nodes in the geometry. The default is NodeArray((0,)),
2761 an empty NodeArray.
2762 coordinate_system : CoordinateSystemArray, optional
2763 The set of coordinate systems in the geometry. The default is
2764 CoordinateSystemArray((0,)), an empty CoordinateSystemArray
2765 traceline : TracelineArray, optional
2766 The set of tracelines defined in the geometry. The default is
2767 TracelineArray((0,)), an empty TracelineArray
2768 element : ElementArray, optional
2769 The set of elements defined in the geometry. The default is
2770 ElementArray((0,)), an empty element array.
2772 Returns
2773 -------
2774 None.
2776 """
2777 self.node = node.flatten()
2778 self.coordinate_system = coordinate_system.flatten()
2779 self.traceline = traceline.flatten()
2780 self.element = element.flatten()
2782 def __repr__(self):
2783 string_out = '\n'.join(val.capitalize(
2784 ) + '\n' + repr(self.__dict__[val]) for val in ['node', 'coordinate_system', 'traceline', 'element'])
2785 return string_out
2787 def reduce(self, node_list):
2788 """
2789 Reduce the geometry to only contain nodes in `node_list`
2791 Elements and tracelines will only be kept if all nodes in each element
2792 or traceline are in `node_list`. Coordinate systems will only be kept
2793 if they are required by a node in `node_list`
2795 Parameters
2796 ----------
2797 node_list : iterable
2798 An iterable of integer node id numbers.
2800 Returns
2801 -------
2802 Geometry
2803 A geometry only containing the nodes in `node_list`
2805 """
2806 node = self.node.reduce(node_list).copy()
2807 traceline = self.traceline.reduce(node_list).copy()
2808 element = self.element.reduce(node_list).copy()
2809 css = np.unique(np.concatenate((np.unique(node.def_cs), node.disp_cs)))
2810 coordinate_system = self.coordinate_system(css).copy()
2811 return Geometry(node, coordinate_system, traceline, element)
2813 def add_traceline(self, connectivity, id=None, description='', color=1):
2814 """
2815 Adds a traceline to the geometry
2817 Parameters
2818 ----------
2819 connectivity : iterable
2820 Iterable containing the node_ids to connect with the traceline
2821 id : TYPE, optional
2822 The id number of the new traceline. The default is None, which
2823 results in the new traceline having an id of one greater than the
2824 current maximum traceline it.
2825 description : str, optional
2826 A string description of the new traceline. The default is ''.
2827 color : int, optional
2828 An integer corresponding to the color of the new traceline.
2829 The default is 1.
2831 Returns
2832 -------
2833 None. Modifications are made in-place to the current geometry.
2835 """
2836 if id is None:
2837 if self.traceline.size > 0:
2838 id = np.max(self.traceline.id) + 1
2839 else:
2840 id = 1
2841 new_traceline = traceline_array((id,), description, color,
2842 [np.array(connectivity)])
2843 self.traceline = np.concatenate((self.traceline, new_traceline), axis=0)
2845 def add_element(self, elem_type, connectivity, id=None, color=1):
2846 if id is None:
2847 if self.element.size > 0:
2848 id = np.max(self.element.id) + 1
2849 else:
2850 id = 1
2851 new_element = element_array((id,), elem_type, color, [np.array(connectivity)])
2852 self.element = np.concatenate((self.element, new_element))
2854 @classmethod
2855 def from_unv(cls, unv_dict):
2856 """
2857 Create a geometry from universal file format data from readunv
2859 Parameters
2860 ----------
2861 unv_dict : dict
2862 Dictionary containing data from read_unv
2864 Returns
2865 -------
2866 Geometry
2867 Geometry object read from the unv data
2869 """
2870 try:
2871 coordinate_systems = CoordinateSystemArray.from_unv(unv_dict, combine=True)
2872 except KeyError:
2873 coordinate_systems = CoordinateSystemArray((0,))
2874 # Nodes
2875 try:
2876 nodes = NodeArray.from_unv(unv_dict, combine=True)
2877 cs_indices = coordinate_systems.id
2878 for node_index, node in nodes.ndenumerate():
2879 cs_index = np.where(cs_indices == node.def_cs)
2880 matching_cs = coordinate_systems[cs_index]
2881 if matching_cs.size == 0:
2882 print('Warning: Coordinate System {:} not found in universal file structure. Adding Global Cartesian with ID {:}'.format(
2883 node.def_cs, node.def_cs))
2884 coordinate_systems = np.concatenate(
2885 (coordinate_systems, coordinate_system_array((node.def_cs,))))
2886 cs_indices = coordinate_systems.id
2887 # cs_index = np.where(cs_indices == node.def_cs)
2888 # node.coordinate = local_coord(coordinate_systems[cs_index][0],node.coordinate)
2889 except KeyError:
2890 nodes = NodeArray((0,))
2891 try:
2892 elements = ElementArray.from_unv(unv_dict, combine=True)
2893 except KeyError:
2894 elements = ElementArray((0,))
2895 try:
2896 tracelines = TracelineArray.from_unv(unv_dict, combine=True)
2897 except KeyError:
2898 tracelines = TracelineArray((0,))
2900 return cls(nodes, coordinate_systems, tracelines, elements)
2902 from_uff = from_unv
2904 @staticmethod
2905 def write_excel_template(path_to_xlsx):
2906 """
2907 Writes an Excel File Template for Creating Geometry
2909 Parameters
2910 ----------
2911 path_to_xlsx : string
2912 Path to write xlsx Excel file
2914 Returns
2915 -------
2916 Nothing
2918 Notes
2919 -----
2920 See documentation for from_excel_template for instructions on filling
2921 out the template to create a geometry.
2922 """
2924 # Create a new workbook
2925 workbook = openpyxl.Workbook()
2927 # Add multiple sheets to the workbook
2928 sheet_names = ['Coordinate Systems', 'Nodes', 'Elements', 'Trace Lines']
2929 for sheet_name in sheet_names:
2930 workbook.create_sheet(sheet_name)
2932 # Remove the default sheet created at the start
2933 del workbook['Sheet']
2935 # Define column names for each sheet
2936 columns = {
2937 'Coordinate Systems': ['ID', 'Name', 'Color', 'Type', 'X Location', 'Y Location', 'Z Location', 'Rotation 1 Angle', 'Rotation 1 Axis', 'Rotation 2 Angle', 'Rotation 2 Axis', 'Rotation 3 Angle', 'Rotation 3 Axis'],
2938 'Nodes': ['ID', 'Color', 'X Location', 'Y Location', 'Z Location', 'Displacement CS', 'Definition CS'],
2939 'Elements': ['ID', 'Color', 'Type', 'Node 1', 'Node 2', 'Node 3', 'Node 4', 'Node 5'],
2940 'Trace Lines': ['ID', 'Description', 'Color', 'Node 1', 'Node 2', 'Node 3', 'Node 4', 'Node 5']
2941 }
2943 # Create data validations
2944 data_validation_ID = openpyxl.worksheet.datavalidation.DataValidation(type="whole", operator='greaterThan', formula1='-1',allow_blank=True,error='Must be interger greater than -1',errorStyle='stop',showErrorMessage=True)
2945 data_validation_Color = openpyxl.worksheet.datavalidation.DataValidation(type="list", formula1='"Black,Blue,Gray Blue,Light Blue,Cyan,Dark Olive,Dark Green,Green,Yellow,Golden Orange,Orange,Red,Magenta,Light Magenta,Pink,White"',allow_blank=True,error='Must be value from list',errorStyle='stop',showErrorMessage=True)
2946 data_validation_CSType = openpyxl.worksheet.datavalidation.DataValidation(type="list", formula1='"Cartesian,Cylindrical,Sypherical"',allow_blank=True,error='Must be value from list',errorStyle='stop',showErrorMessage=True)
2947 data_validation_Axis = openpyxl.worksheet.datavalidation.DataValidation(type="list", formula1='"X,Y,Z"',allow_blank=True,error='Must be value from list',errorStyle='stop',showErrorMessage=True)
2948 data_validation_ElemType = openpyxl.worksheet.datavalidation.DataValidation(type="list", formula1='"beam,triangle,quadrilateral,tetrahedral"',allow_blank=True,error='Must be value from list',errorStyle='stop',showErrorMessage=True)
2949 data_validation_Node = openpyxl.worksheet.datavalidation.DataValidation(type="whole", operator='greaterThan', formula1='-1',allow_blank=True,error='Must be value from list',errorStyle='stop',showErrorMessage=True)
2950 data_validation_Decimal = openpyxl.worksheet.datavalidation.DataValidation(type="decimal",allow_blank=True,error='Must be decimal',errorStyle='stop',showErrorMessage=True)
2952 # Specify Columns for data validation
2953 table_height = 10
2954 data_validation_ID_A = copy.deepcopy(data_validation_ID)
2955 data_validation_ID_A.add('A2:A'+str(table_height+1))
2956 data_validation_Color_B = copy.deepcopy(data_validation_Color)
2957 data_validation_Color_B.add('B2:B'+str(table_height+1))
2958 data_validation_Color_C = copy.deepcopy(data_validation_Color)
2959 data_validation_Color_C.add('C2:C'+str(table_height+1))
2960 data_validation_CSType.add('D2:D'+str(table_height+1))
2961 data_validation_Axis.add('I2:I'+str(table_height+1))
2962 data_validation_Axis.add('K2:K'+str(table_height+1))
2963 data_validation_Axis.add('M2:M'+str(table_height+1))
2964 data_validation_Decimal_Coordinate_System = data_validation_Decimal
2965 data_validation_Decimal_Coordinate_System.add('E2:H'+str(table_height+1))
2966 data_validation_Decimal_Coordinate_System.add('J2:J'+str(table_height+1))
2967 data_validation_Decimal_Coordinate_System.add('L2:L'+str(table_height+1))
2968 data_validation_Decimal_Nodes = copy.deepcopy(data_validation_Decimal)
2969 data_validation_Decimal_Nodes.add('C2:E'+str(table_height+1))
2970 data_validation_ID_Nodes = copy.deepcopy(data_validation_ID)
2971 data_validation_ID_Nodes.add('F2:G'+str(table_height+1))
2972 data_validation_ElemType.add('C2:C'+str(table_height+1))
2973 data_validation_Node.add('D2:H'+str(table_height+1))
2975 # Add data and tables to each sheet
2976 for sheet_name in sheet_names:
2977 sheet = workbook[sheet_name]
2979 # Define the column headers
2980 headers = columns[sheet_name]
2982 # Write the headers into the first row
2983 for col_num, header in enumerate(headers, start=1):
2984 sheet.cell(row=1, column=col_num, value=header)
2986 table_range = f'A1:{chr(65 + len(headers) - 1)}{table_height + 1}'
2988 table = openpyxl.worksheet.table.Table(displayName=sheet_name.replace(' ','_'), ref=table_range)
2990 # Apply table style
2991 style = openpyxl.worksheet.table.TableStyleInfo(name="TableStyleMedium9", showFirstColumn=False,showLastColumn=False, showRowStripes=True, showColumnStripes=False)
2992 table.tableStyleInfo = style
2994 # Add table to the worksheet
2995 sheet.add_table(table)
2997 # Apply Data Validations for "Coordinate Systems" Sheet
2998 sheet = workbook['Coordinate Systems']
2999 sheet.add_data_validation(data_validation_ID_A)
3000 sheet.add_data_validation(data_validation_Color_C)
3001 sheet.add_data_validation(data_validation_CSType)
3002 sheet.add_data_validation(data_validation_Decimal_Coordinate_System)
3003 sheet.add_data_validation(data_validation_Axis)
3005 # Apply Data Validations for "Nodes" Sheet
3006 sheet = workbook['Nodes']
3007 sheet.add_data_validation(data_validation_ID_A)
3008 sheet.add_data_validation(data_validation_ID_Nodes)
3009 sheet.add_data_validation(data_validation_Color_B)
3010 sheet.add_data_validation(data_validation_Decimal_Nodes)
3012 # Apply Data Validations for "Elements" Sheet
3013 sheet = workbook['Elements']
3014 sheet.add_data_validation(data_validation_ID_A)
3015 sheet.add_data_validation(data_validation_Color_B)
3016 sheet.add_data_validation(data_validation_ElemType)
3017 sheet.add_data_validation(data_validation_Node)
3019 # Apply Data Validations for "Trace Lines" Sheet
3020 sheet = workbook['Trace Lines']
3021 sheet.add_data_validation(data_validation_ID_A)
3022 sheet.add_data_validation(data_validation_Color_C)
3023 sheet.add_data_validation(data_validation_Node)
3025 # Save the workbook
3026 workbook.save(path_to_xlsx)
3028 @classmethod
3029 def from_excel_template(cls, path_to_xlsx):
3030 """
3031 Create a geometry from Excel file template
3033 Parameters
3034 ----------
3035 path_to_xlsx : string
3036 Path to xlsx Excel file containing geometry information
3038 Returns
3039 -------
3040 Geometry
3041 Geometry object created from the Excel file
3043 Notes
3044 -----
3045 To use this function, first save out an excel template file using the
3046 `write_excel_template` function. This will construct an excel
3047 workbook with four worksheets on which the different portions of the
3048 Geometry are defined.
3050 On the Coordinate Systems tab, users will define the various global and
3051 local coordinate systems in their geometry. Each geometry requires
3052 an ID number. A Name can optionally be given. The Color should be
3053 selectred from the dropdown list of colors in the Excel template. The
3054 type of the coordinate system should be selected from the dropdown list of
3055 coordinate system types.
3057 The origin of the coordinate system can be specified with the X Location,
3058 Y Location, Z Location columns. Then rotations of the coordinate system
3059 can be specified using rotations about axes. Up to three axes and angles
3060 can be specified to create arbitrary compound rotations. The rotation
3061 axes should be X, Y, or Z, and the rotation angles are in degrees.
3063 On the Nodes tab, all tabs must be filled out. ID numbers must be unique.
3064 Colors should be selected from the dropdown list in the Excel template. The
3065 position of the node is specified using the X Location, Y Location, Z
3066 Location columns. Each node has a displacement coordinate system in which
3067 its position is defined, and a a definition coordinate system in which
3068 its displacements are defined. These columns should consist of integers
3069 corresponding to ID numbers from the Coordinate Systems tab.
3071 On the Elements tab, elements connecting nodes are defined. The ID
3072 column must consist of unique integer identifiers. The Color tab should
3073 be selected from the dropdown list. The type can be selected from the dropdown list.
3074 If the type column
3075 is empty, an element type based on the number of connections given will
3076 be used. Defult element types for connection length of 2 is "Type 21 - Linear Beam",
3077 for a connection length of 3 is "Type 41 - Plane Stress Linear Triangle",
3078 and for a connection length of 4 is "Type 44 - Plane Stress Linear Quadrilateral".
3079 The columns Node 1 through Node 4 contain the nodes in each element. Only the required
3080 number of nodes must be filled out (e.g. a tri element would only contain
3081 3 nodes, so only columns Node 1, Node 2, and Node 3 would be filled).
3083 On the Trace Lines tab, lines connecting the nodes are defined. The
3084 ID column must consist of unique integer identifiers. The Description
3085 column contains a string description of the line. The Color should be
3086 selected from the dropdown list. The
3087 Node 1 through Node 5 columns should contain the nodes for each line.
3088 Only the number of nodes in the line must be filled out, so if a line
3089 only connects 5 nodes, only Node 1 though Node 5 must be filled. More columns can
3090 can be added for additional Nodes if longer tracelines are needed.
3093 """
3095 # % Read Geometry Information from Excel File
3096 df_coordinate_systems = pd.read_excel(io = path_to_xlsx,sheet_name='Coordinate Systems')
3097 df_nodes = pd.read_excel(io = path_to_xlsx,sheet_name='Nodes')
3098 df_elements = pd.read_excel(io = path_to_xlsx,sheet_name='Elements')
3099 df_trace_lines = pd.read_excel(io = path_to_xlsx,sheet_name='Trace Lines')
3101 color_number_dict = {
3102 'Black':1,
3103 'Blue':2,
3104 'Gray Blue':3,
3105 'Light Blue':4,
3106 'Cyan':5,
3107 'Dark Olive':6,
3108 'Dark Green':7,
3109 'Green':8,
3110 'Yellow':9,
3111 'Golden Orange':10,
3112 'Orange':11,
3113 'Red':12,
3114 'Magenta':13,
3115 'Light Magenta':14,
3116 'Pink':15,
3117 'White':16,
3118 }
3120 cs_type_dict = {
3121 'Cartesian':0,
3122 'Cylindrical':1,
3123 'Sypherical':2,
3124 }
3126 df_coordinate_systems = df_coordinate_systems.replace({'Color':color_number_dict})
3127 df_coordinate_systems = df_coordinate_systems.replace({'Type':cs_type_dict})
3128 df_nodes = df_nodes.replace({'Color':color_number_dict})
3129 df_elements = df_elements.replace({'Color':color_number_dict})
3130 df_trace_lines = df_trace_lines.replace({'Color':color_number_dict})
3132 # % Add Coordinate Systems
3133 # Start with an empty coordinate system
3134 coordinate_systems = coordinate_system_array([])
3135 for row_index, df_coordinate_system in df_coordinate_systems.iterrows():
3136 # Create Transformation Matrix
3137 T = np.eye(3)
3138 for rotation in range(1,(len(df_coordinate_system)-7)//2+1):
3139 rotation_angle = df_coordinate_system['Rotation ' + str(rotation) + ' Angle']
3140 rotation_axis = df_coordinate_system['Rotation ' + str(rotation) + ' Axis']
3141 rotation_axis_dict = {'X':0,'Y':1,'Z':2}
3142 if pd.notnull(rotation_angle) and pd.notnull(rotation_axis):
3143 T = R(rotation_axis_dict[rotation_axis],rotation_angle,degrees=True).T@T
3145 # Add Coordinate System to Geometry
3146 origin = np.array((df_coordinate_system['X Location'],df_coordinate_system['Y Location'],df_coordinate_system['Z Location']))[np.newaxis,:]
3147 coordinate_systems = np.concatenate((coordinate_systems,coordinate_system_array(id=df_coordinate_system['ID'], name=df_coordinate_system['Name'], color=df_coordinate_system['Color']-1, cs_type=df_coordinate_system['Type'], matrix = np.concatenate((T,origin),axis=0))),axis=None)
3149 # % Add Nodes
3150 # print(df_nodes['Color'])
3151 nodes = node_array(
3152 id=np.array(df_nodes['ID']),
3153 coordinate=np.concatenate(
3154 [np.array(df_nodes['X Location'])[:,np.newaxis],np.array(df_nodes['Y Location'])[:,np.newaxis],np.array(df_nodes['Z Location'])[:,np.newaxis]],axis=-1),
3155 color=np.array(df_nodes['Color'])-1,
3156 def_cs=np.array(df_nodes['Definition CS']),
3157 disp_cs=np.array(df_nodes['Displacement CS']))
3159 # % Add Elements
3160 elements = element_array([])
3161 element_connection_length_dict = {
3162 2:21, # 21 : 'Linear Beam'
3163 3:41, # 41 : 'Plane Stress Linear Triangle'
3164 4:44, # 44 : 'Plane Stress Linear Quadrilateral'
3165 }
3166 for row_index, df_element in df_elements.iterrows():
3167 connectivity = df_element[3:].values
3168 connectivity = connectivity[pd.notnull(connectivity)]
3169 if isinstance(df_element['Type'],str):
3170 element_type = _exodus_elem_type_map[df_element['Type'].lower()]
3171 elif isinstance(df_element['Type'],int):
3172 element_type = df_element['Type']
3173 elif pd.isnull(df_element['Type']):
3174 element_type = element_connection_length_dict[len(connectivity)]
3175 new_element = element_array(id=df_element['ID'], type=element_type, color=df_element['Color']-1,connectivity=connectivity)
3176 elements = np.concatenate((elements,new_element),axis=None)
3178 # % Add Tracelines
3179 tracelines = traceline_array([])
3180 for row_index, df_trace_line in df_trace_lines.iterrows():
3181 connectivity=df_trace_line[3:].values
3182 connectivity = connectivity[pd.notnull(connectivity)]
3183 new_traceline = traceline_array(connectivity=connectivity,id=df_trace_line['ID'],description=df_trace_line['Description'],color=df_trace_line['Color']-1)
3184 tracelines = np.concatenate((tracelines,new_traceline),axis=None)
3186 return cls(nodes, coordinate_systems, tracelines, elements)
3188 @classmethod
3189 def from_exodus(cls, exo: Exodus, blocks=None, local=False,
3190 preferred_local_orientation=np.array((0.0, 0.0, 1.0)),
3191 secondary_preferred_local_orientation=np.array((1.0, 0.0, 0.0)),
3192 local_nodes=None):
3193 """
3194 Generate a geometry from exodus file data
3196 Parameters
3197 ----------
3198 exo : Exodus or ExodusInMemory
3199 The exodus data from which geometry will be created.
3200 blocks : iterable, optional
3201 Iterable containing a set of block ids to import when creating the
3202 geometry. The default is None, which uses all blocks in the model.
3203 local : bool, optional
3204 Flag to specify whether or not to create local coordinate systems
3205 for each node in the model. This can be useful when creating
3206 instrumnetation positions from a finite element model, where the
3207 sensor will be oriented perpendicular to the surface it is mounted
3208 on. The default is False, which returns all data in a single
3209 global coordinate system.
3210 preferred_local_orientation : np.ndarray, optional
3211 A preferred direction for the local coordinate system. The first
3212 constraint is that the local Z+ axis is perpendicular to the
3213 surface. The coordinate system will then try to align itself as
3214 much as it can with this direction. The default is
3215 np.array((0.0,0.0,1.0)), which points the local coordinate systems
3216 along the global Z+ axis.
3217 secondary_preferred_local_orientation : np.ndarray, optional
3218 A secondary preferred direction is only used if the surface normal
3219 direction is parallel to the primary preferred_local_orientation.
3220 The default is np.array((1.0,0.0,0.0)), which points the local
3221 coordinate system along the local Z+ axis.
3222 local_nodes : np.ndarray, optional
3223 If specified, only create local coordinate systems at the specified
3224 nodes.
3226 Returns
3227 -------
3228 Geometry
3229 A geometry consisting of the finite element nodes and element
3230 connectivity.
3232 """
3233 if isinstance(exo, Exodus):
3234 exo = exo.load_into_memory(close=False, variables=[], timesteps=[], blocks=blocks)
3235 # Get nodes
3236 coordinates = exo.nodes.coordinates
3237 node_map = np.arange(
3238 coordinates.shape[0]) + 1 if exo.nodes.node_num_map is None else exo.nodes.node_num_map
3239 nodes = node_array(node_map, coordinates)
3240 tl_connectivity = []
3241 elem_connectivity = []
3242 tl_color = []
3243 elem_color = []
3244 elem_types = []
3245 if blocks is None:
3246 blocks = exo.blocks
3247 else:
3248 blocks = [block for block in exo.blocks if block.id in blocks]
3249 for i, block in enumerate(blocks):
3250 color_index = (i % 14) + 1
3251 try:
3252 connectivity = node_map[block.connectivity]
3253 except AttributeError:
3254 continue
3255 elem_type = block.elem_type
3256 if connectivity.shape[-1] <= 1:
3257 # Skip conmasses
3258 continue
3259 for element in connectivity:
3260 elem_connectivity.append(element)
3261 elem_color.append(color_index)
3262 elem_types.append(_exodus_elem_type_map[elem_type.lower()])
3263 tls = traceline_array(np.arange(len(tl_connectivity)) + 1,
3264 color=tl_color,
3265 connectivity=tl_connectivity)
3266 elems = element_array(np.arange(len(elem_connectivity)) + 1,
3267 elem_types, elem_color, elem_connectivity)
3268 if local:
3269 if local_nodes is not None:
3270 node_map_dict = {val: index for index, val in enumerate(node_map)}
3271 local_node_indices = np.array([node_map_dict[node] for node in local_nodes])
3272 local_node_index_map = {node_index: index for index, node_index in enumerate(local_node_indices)}
3273 preferred_local_orientation /= np.linalg.norm(preferred_local_orientation)
3274 secondary_preferred_local_orientation /= np.linalg.norm(
3275 secondary_preferred_local_orientation)
3276 # Go through each node, find the block that it's in, get its elements
3277 normal_sum = np.zeros(
3278 (coordinates.shape[0] if local_nodes is None else len(local_nodes),
3279 3))
3280 for block in blocks:
3281 if block.connectivity.shape[-1] <= 1:
3282 # Skip conmasses
3283 continue
3284 full_block_node_positions = coordinates[block.connectivity, :]
3285 normal_vectors = np.cross(full_block_node_positions[:, 2, :] - full_block_node_positions[:, 0, :],
3286 full_block_node_positions[:, 1, :] - full_block_node_positions[:, 0, :])
3287 normal_vectors /= np.linalg.norm(normal_vectors, axis=-1, keepdims=True)
3288 block_nodes = np.unique(block.connectivity)
3289 if local_nodes is not None:
3290 block_nodes = local_node_indices[np.isin(local_node_indices, block_nodes)]
3291 for index in block_nodes:
3292 node_elements = np.any(block.connectivity == index, axis=-1)
3293 node_normals = normal_vectors[node_elements]
3294 if local_nodes is None:
3295 normal_sum[index] += np.sum(node_normals, axis=0)
3296 else:
3297 normal_sum[local_node_index_map[index]] += np.sum(node_normals, axis=0)
3298 # Put small nodes to reasonable values to avoid divide by zero
3299 normal_sum[np.linalg.norm(normal_sum, axis=-1) < 1e-14] = preferred_local_orientation
3300 node_rotations = np.zeros((3,) + normal_sum.shape)
3301 node_rotations[2] = -normal_sum / np.linalg.norm(normal_sum, axis=-1, keepdims=True)
3302 # Find the nodes that are basically in the direction of the preferred direction
3303 nodes_in_preferred_direction = np.abs(
3304 np.einsum('ij,j->i', node_rotations[2], preferred_local_orientation)) > 0.99
3305 # Compute cross products if we are not parallel with the preferred orientation
3306 node_rotations[1, ~nodes_in_preferred_direction] = np.cross(
3307 node_rotations[2, ~nodes_in_preferred_direction], preferred_local_orientation)
3308 node_rotations[1, nodes_in_preferred_direction] = np.cross(
3309 node_rotations[2, nodes_in_preferred_direction], secondary_preferred_local_orientation)
3310 node_rotations[1] /= np.linalg.norm(node_rotations[1], axis=-1, keepdims=True)
3311 node_rotations[0] = np.cross(node_rotations[1], node_rotations[2])
3312 node_rotations[0] /= np.linalg.norm(node_rotations[0], axis=-1, keepdims=True)
3313 node_rotations = np.moveaxis(node_rotations, 0, 1)
3314 # Create the coordinate system matrix
3315 cs_matrix = np.concatenate((node_rotations, np.zeros(
3316 (node_rotations.shape[0], 1, node_rotations.shape[-1]))), axis=1)
3317 # Put a global coordinate system at the end
3318 cs_matrix = np.concatenate((cs_matrix, np.eye(4, 3)[np.newaxis]), axis=0)
3319 if local_nodes is None:
3320 cs_ids = np.concatenate((node_map, (node_map.max() + 1,)), axis=0)
3321 cs_descriptions = ['Node {:} Disp'.format(id) for id in node_map] + ['Global']
3322 css = coordinate_system_array(cs_ids, cs_descriptions, 1, 0, cs_matrix)
3323 nodes.def_cs = cs_ids[-1]
3324 nodes.disp_cs = cs_ids[:-1]
3325 else:
3326 cs_ids = np.concatenate((local_nodes, (node_map.max() + 1,)), axis=0)
3327 cs_descriptions = ['Node {:} Disp'.format(id) for id in local_nodes] + ['Global']
3328 css = coordinate_system_array(cs_ids, cs_descriptions, 1, 0, cs_matrix)
3329 nodes.def_cs = cs_ids[-1]
3330 nodes.disp_cs = cs_ids[-1]
3331 nodes.disp_cs[local_node_indices] = cs_ids[:-1]
3332 else:
3333 css = coordinate_system_array((1,))
3334 return Geometry(nodes, css, tls, elems)
3336 @classmethod
3337 def camera_visualization(cls, K, RT, image_size, size=1, colors=1):
3338 """
3339 Create a geometry used to visualize a camera specified with K and RT
3341 Parameters
3342 ----------
3343 K : ndarray
3344 A (...,3,3) matrix containing the intrinsic parameters of the cameras
3345 RT : ndarray
3346 A (...,3,4) matrix containing the extrinsic parameters of the cameras
3347 image_size : ndarray
3348 A (...,2) matrix containing the width,height of each camera
3349 size : float
3350 The distance that the rays will project from the pinhole of the camera
3351 colors : int or ndarray
3352 The colors assigned to each camera
3354 Returns
3355 -------
3356 geometry : cls
3357 A geoemtry containing the cameras that can be used for visualization
3359 """
3360 if K.ndim == 2:
3361 K = K[np.newaxis, ...]
3362 if K.ndim > 3:
3363 K = K.reshape(-1, 3, 3)
3364 if RT.ndim == 2:
3365 RT = RT[np.newaxis, ...]
3366 if RT.ndim > 3:
3367 RT = RT.reshape(-1, 3, 4)
3368 if image_size.ndim == 1:
3369 image_size = image_size[np.newaxis, ...]
3370 if image_size.ndim > 2:
3371 image_size = image_size.reshape(-1, 2)
3372 if np.array(colors).ndim == 0:
3373 colors = colors * np.ones(K.shape[0])
3374 if np.array(colors).ndim > 1:
3375 colors = np.array(colors).flatten()
3377 corners = np.array([(0, 0),
3378 (1, 0),
3379 (1, 1),
3380 (0, 1)])
3382 geometry = cls(coordinate_system=coordinate_system_array((1,)))
3384 for i, (thisK, thisRT, this_size, this_color) in enumerate(zip(K, RT, image_size, colors)):
3385 pixel_corners = corners * this_size
3387 # Compute the positions
3388 image_corners = point_on_pixel(thisK, thisRT, pixel_corners, size)
3389 pinhole = -thisRT[:3, :3].T @ thisRT[..., :3, 3:]
3390 # Create a geometry
3391 pinhole_node = node_array(((i + 1) * 10,), pinhole.T, disp_cs=i + 2, color=this_color)
3392 pixel_nodes = node_array((i + 1) * 10 + np.arange(4) + 1,
3393 image_corners, disp_cs=i + 2, color=this_color)
3394 cs = coordinate_system_array((i + 2,), color=this_color)
3395 cs.matrix[..., :3, :3] = thisRT[:3, :3]
3396 geometry.node = np.concatenate((geometry.node, pinhole_node, pixel_nodes))
3397 geometry.coordinate_system = np.concatenate((geometry.coordinate_system, cs))
3398 geometry.add_traceline(np.array([0, 1, 2, 3, 4, 1]) + 10 * (i + 1), color=this_color)
3399 geometry.add_traceline(np.array([0, 2]) + 10 * (i + 1), color=this_color)
3400 geometry.add_traceline(np.array([0, 3]) + 10 * (i + 1), color=this_color)
3401 geometry.add_traceline(np.array([0, 4]) + 10 * (i + 1), color=this_color)
3402 return geometry
3404 def compress_ids(self, compress_nodes=True, compress_elements=True,
3405 compress_tracelines=True, compress_cs=True,
3406 return_maps=False):
3407 """
3408 Compresses ID numbers to make node, element, etc. contiguous
3410 Parameters
3411 ----------
3412 compress_nodes : bool, optional
3413 If True, compress the node ids. The default is True.
3414 compress_elements : bool, optional
3415 If True, compress the element ids. The default is True.
3416 compress_tracelines : bool, optional
3417 If True, compress the traceline ids. The default is True.
3418 compress_cs : bool, optional
3419 If True, compress the coordinate system ids. The default is True.
3420 return_maps : bool, optional
3421 If True, return id maps for nodes, elements, tracelines, and
3422 coordinate systems. Maps will be equal to `None` if no compression
3423 was done on that field. The default is False.
3425 Returns
3426 -------
3427 mapped_geometry : Geometry
3428 Geometry with contiguous id numbers for the selected fields.
3429 node_map : id_map or `None`
3430 Mapping from the old set of nodes to the new set of nodes, only
3431 returned if `return_maps` is True. If `compress_nodes` is False,
3432 this will be `None`.
3433 traceline_map : id_map or `None`
3434 Mapping from the old set of tracelines to the new set of tracelines, only
3435 returned if `return_maps` is True. If `compress_tracelines` is False,
3436 this will be `None`.
3437 element_map : id_map or `None`
3438 Mapping from the old set of elements to the new set of elements, only
3439 returned if `return_maps` is True. If `compress_elements` is False,
3440 this will be `None`.
3441 cs_map : id_map or `None`
3442 Mapping from the old set of coordinate systems to the new set of
3443 coordinate systems, only
3444 returned if `return_maps` is True. If `compress_cs` is False,
3445 this will be `None`.
3447 """
3448 if compress_nodes:
3449 to_nodes = np.arange(self.node.size)+1
3450 from_nodes = self.node.id.flatten()
3451 node_map = id_map(from_nodes, to_nodes)
3452 else:
3453 node_map = None
3454 if compress_cs:
3455 to_cs = np.arange(self.coordinate_system.size)+1
3456 from_cs = self.coordinate_system.id.flatten()
3457 cs_map = id_map(from_cs, to_cs)
3458 else:
3459 cs_map = None
3460 if compress_elements:
3461 to_elements = np.arange(self.element.size)+1
3462 from_elements = self.coordinate_system.id.flatten()
3463 element_map = id_map(from_elements, to_elements)
3464 else:
3465 element_map = None
3466 if compress_tracelines:
3467 to_tls = np.arange(self.traceline.size)+1
3468 from_tls = self.coordinate_system.id.flatten()
3469 traceline_map = id_map(from_tls, to_tls)
3470 else:
3471 traceline_map = None
3472 mapped_geometry = self.map_ids(node_map, traceline_map, element_map, cs_map)
3473 if return_maps:
3474 return mapped_geometry, node_map, traceline_map, element_map, cs_map
3475 else:
3476 return mapped_geometry
3478 @classmethod
3479 def from_imat_struct(cls, imat_fem_struct):
3480 """
3481 Constructs a Geometry from an imat_fem class saved to a Matlab structure
3483 In IMAT, a structure can be created from an `imat_fem` by using the get()
3484 function. This can then be saved to a .mat file and loaded using
3485 `scipy.io.loadmat`. The output from loadmat can be passed into this function
3487 Parameters
3488 ----------
3489 imat_fem_struct : np.ndarray
3490 structure from loadmat containing data from an imat_fem
3492 Returns
3493 -------
3494 Geometry
3495 Geometry constructed from the data in the imat structure
3497 """
3498 node_struct = imat_fem_struct['node'][0, 0]
3499 cs_array = node_struct['cs'][0, 0].reshape(-1, 3)
3500 nodes = node_array(id=node_struct['id'][0, 0].flatten(),
3501 coordinate=node_struct['coord'][0, 0].reshape(-1, 3),
3502 color=node_struct['color'][0, 0].flatten(),
3503 def_cs=cs_array[:, 0], disp_cs=cs_array[:, 1])
3504 cs_struct = imat_fem_struct['cs'][0, 0]
3505 cs_ids = cs_struct['id'][0, 0].flatten()
3506 if cs_ids.size > 0:
3507 cs_names = np.concatenate(cs_struct['name'][0, 0].flatten()).flatten()
3508 else:
3509 cs_names = []
3510 css = coordinate_system_array(id=cs_ids,
3511 name=cs_names,
3512 color=cs_struct['color'][0, 0].flatten(),
3513 cs_type=cs_struct['type'][0, 0].flatten(),
3514 matrix=cs_struct['matrix'][0, 0].reshape(4, 3, -1).transpose(2, 0, 1))
3515 elem_struct = imat_fem_struct['elem'][0, 0]
3516 elems = element_array(id=elem_struct['id'][0, 0].flatten(),
3517 type=elem_struct['type'][0, 0].flatten(),
3518 color=elem_struct['color'][0, 0].flatten(),
3519 connectivity=[val.flatten() for val in elem_struct['conn'][0, 0].flatten()])
3520 tl_struct = imat_fem_struct['tl'][0, 0]
3521 tl_ids = tl_struct['id'][0, 0].flatten()
3522 if tl_ids.size > 0:
3523 tl_names = np.concatenate(tl_struct['desc'][0, 0].flatten()).flatten()
3524 else:
3525 tl_names = []
3526 tls = traceline_array(id=tl_ids,
3527 description=tl_names,
3528 color=tl_struct['color'][0, 0].flatten(),
3529 connectivity=[val.flatten() for val in tl_struct['conn'][0, 0].flatten()])
3530 return cls(nodes, css, tls, elems)
3532 def global_node_coordinate(self, node_ids=None):
3533 """
3534 Position of the Geometry's nodes in the global coordinate system
3536 Parameters
3537 ----------
3538 node_ids : np.ndarray
3539 An array of node id numbers to keep in the output coordinate. If
3540 not specified, all nodes will be returned in the order they are in
3541 the model.
3544 Returns
3545 -------
3546 np.ndarray
3547 numpy array with shape (n,3) where n is the number of nodes in the
3548 geometry. This array contains the 3d position of each node in the
3549 global coordinate system.
3551 """
3552 if node_ids is None:
3553 node = self.node
3554 else:
3555 node = self.node(node_ids)
3556 cs_array = self.coordinate_system(node.def_cs)
3557 return node.global_coordinate(cs_array)
3559 def node_by_global_position(self, global_position_array):
3560 """
3561 Select node by closest position
3563 Parameters
3564 ----------
3565 global_position_array : np.ndarray
3566 A (...,3) shape array containing positions of nodes to keep
3568 Returns
3569 -------
3570 NodeArray
3571 NodArray containing nodes that were closest to the positions in
3572 position_array.
3574 """
3575 node_differences = np.linalg.norm(
3576 self.global_node_coordinate() - global_position_array[(Ellipsis,) + (np.newaxis,) * self.node.ndim + (slice(None),)], axis=-1)
3577 node_differences = node_differences.reshape(*global_position_array.shape[:-1], -1)
3578 node_indices = np.argmin(
3579 node_differences,
3580 axis=-1)
3581 return self.node.flatten()[node_indices]
3583 def node_by_global_grid(self,grid_spacing, x_min = None, y_min = None,
3584 z_min = None, x_max = None, y_max = None,
3585 z_max = None):
3586 """
3587 Selects nodes in a grid
3589 Parameters
3590 ----------
3591 grid_spacing : float
3592 Approximate grid spacing between selected nodes
3593 x_min : float, optional
3594 Minimum X-value of the grid.
3595 The default is the minimum node x-coordinate.
3596 y_min : float, optional
3597 Minimun Y-value of the grid.
3598 The default is the minimum node y-coordinate.
3599 z_min : float, optional
3600 Minimum Z-value of the grid.
3601 The default is the minimum node z-coordinate.
3602 x_max : float, optional
3603 Maximum X-value of the grid.
3604 The default is the maximum node x-coordinate.
3605 y_max : float, optional
3606 Maximum Y-value of the grid.
3607 The default is the maximum node y-coordinate.
3608 z_max : float, optional
3609 Maximum Z-value of the grid.
3610 The default is the maximum node z-coordinate.
3612 Returns
3613 -------
3614 NodeArray
3615 NodeArray containing the closest nodes to the grid points
3616 """
3617 coordinate = self.global_node_coordinate()
3618 if x_min is None:
3619 x_min = coordinate[:, 0].min()
3620 if y_min is None:
3621 y_min = coordinate[:, 1].min()
3622 if z_min is None:
3623 z_min = coordinate[:, 2].min()
3624 if x_max is None:
3625 x_max = coordinate[:, 0].max()
3626 if y_max is None:
3627 y_max = coordinate[:, 1].max()
3628 if z_max is None:
3629 z_max = coordinate[:, 2].max()
3630 min_coords = np.array((x_min, y_min, z_min))
3631 max_coords = np.array((x_max, y_max, z_max))
3632 range_coords = max_coords - min_coords
3633 num_grids = np.ceil(range_coords / grid_spacing).astype(int)
3634 # Now we will create the grid in each dimension using linspace
3635 grid_arrays = [np.linspace(min_coord, max_coord, num_coord)
3636 for min_coord, max_coord, num_coord
3637 in zip(min_coords, max_coords, num_grids)]
3638 # We will then use meshgrid to assemble the array grid. We will flatten the
3639 # point dimension and transpose so the array has shape (n_points x 3)
3640 grid_points = np.array(np.meshgrid(*grid_arrays, indexing='ij')).reshape(3, -1).T
3641 # Now we will select nodes that are closest to the points in the grid
3642 candidate_nodes = self.node_by_global_position(grid_points)
3643 # Reduce to only nodes that are within one grid spacing of their requested point
3644 # this reduces the nodes that are selected by points far from the model
3645 candidate_nodes = candidate_nodes[
3646 np.linalg.norm(self.global_node_coordinate(candidate_nodes.id) - grid_points, axis=-1) < grid_spacing]
3647 # Remove any duplicates
3648 candidate_node_ids = np.unique(candidate_nodes.id)
3649 return candidate_nodes(candidate_node_ids)
3651 def global_deflection(self, coordinate_array):
3652 """
3653 Direction of local deflection in the global coordinate system
3655 Parameters
3656 ----------
3657 coordinate_array : CoordinateArray
3658 A list of coordinates for which the global deformation direction
3659 will be computed
3661 Returns
3662 -------
3663 global_deflections : np.ndarray
3664 numpy array with shape (n,3) where n is the number of coordinates
3665 in the specified `coordinate_array`. This array contains the 3d
3666 direction of motion for each coordinate in the global coordinate
3667 system.
3669 """
3670 local_deformations = coordinate_array.local_direction()
3671 ordered_nodes = self.node(coordinate_array.node)
3672 coordinate_systems = self.coordinate_system(ordered_nodes.disp_cs)
3673 points = global_coord(self.coordinate_system(
3674 ordered_nodes.def_cs), ordered_nodes.coordinate)
3675 global_deflections = global_deflection(coordinate_systems, local_deformations, points)
3676 return global_deflections
3678 @staticmethod
3679 def overlay_geometries(geometries, color_override=None, return_node_id_offset=False):
3680 """
3681 Combines several geometries, offsetting the id numbers to avoid conflicts
3683 Parameters
3684 ----------
3685 geometries : iterable
3686 An iterable of geometry objects that will be combined into a single
3687 geometry
3688 color_override : iterable, optional
3689 An iterble of integers specifying colors, which will override the
3690 existing geometry colors. This should have the same length as the
3691 `geometries` input. The default is None, which keeps the original
3692 geometry colors.
3693 return_node_id_offset : bool, optional
3694 Specifies whether or not the applied node offset should be returned,
3695 which is useful if the new node numbers are to be mapped to the old
3696 node numbers. The default is False, which only returns the combined
3697 geometry.
3699 Returns
3700 -------
3701 Geometry
3702 A geometry consisting of a combination of the specified geometries
3703 node_offset
3704 An integer specifying the node id offset applied to avoid conflicts
3706 """
3707 node_ids = np.concatenate([geometry.node.id for geometry in geometries])
3708 cs_ids = np.concatenate([geometry.coordinate_system.id for geometry in geometries])
3709 tl_ids = np.concatenate([geometry.traceline.id for geometry in geometries])
3710 elem_ids = np.concatenate([geometry.element.id for geometry in geometries])
3711 try:
3712 max_node_id = np.max(node_ids)
3713 except ValueError:
3714 max_node_id = 1
3715 try:
3716 max_cs_id = np.max(cs_ids)
3717 except ValueError:
3718 max_cs_id = 1
3719 try:
3720 max_tl_id = np.max(tl_ids)
3721 except ValueError:
3722 max_tl_id = 1
3723 try:
3724 max_elem_id = np.max(elem_ids)
3725 except ValueError:
3726 max_elem_id = 1
3727 node_length = int(np.floor(np.log10(max_node_id)) + 1)
3728 tl_length = int(np.floor(np.log10(max_tl_id)) + 1)
3729 elem_length = int(np.floor(np.log10(max_elem_id)) + 1)
3730 cs_length = int(np.floor(np.log10(max_cs_id)) + 1)
3731 new_geometries = [geometry.modify_ids(
3732 node_change=10**(node_length) * (i + 1),
3733 traceline_change=10**(tl_length) * (i + 1),
3734 element_change=10**(elem_length) * (i + 1),
3735 coordinate_system_change=10**(cs_length) * (i + 1))
3736 for i, geometry in enumerate(geometries)]
3737 if color_override is not None:
3738 for i, color in enumerate(color_override):
3739 new_geometries[i].node.color = color
3740 new_geometries[i].traceline.color = color
3741 new_geometries[i].element.color = color
3742 new_geometries[i].coordinate_system.color = color
3744 new_geometry = new_geometries[0]
3745 for i in range(1, len(new_geometries)):
3746 new_geometry = new_geometry + new_geometries[i]
3748 if return_node_id_offset:
3749 return new_geometry, 10**(node_length)
3750 else:
3751 return new_geometry
3753 def map_coordinate(self,dof_list,to_geometry,node_map=None,
3754 plot_maps = False):
3755 """
3756 Map degrees of freedom from one geometry to the closest aligned on a
3757 second geometry.
3759 Parameters
3760 ----------
3761 dof_list : CoordinateArray
3762 An array of degrees of freedom on the geometry to map.
3763 to_geometry : Geometry
3764 The geometry to which the degrees of freedom will be mapped.
3765 node_map : None, str, or id_map, optional
3766 If specified, it will map the node ids from the present geometry to
3767 the mapping geometry. If the string 'position' is specified, then
3768 the closest node by position on the mapping geometry to the present
3769 geometry will be chosen. To explicitly map between nodes, pass a
3770 sdpy.id_map. If not specified, it is assumed that the node numbers
3771 in the present geometry map to the nodes in the mapping geometry.
3772 plot_maps : bool
3773 If True, it will plot the CoordinateArray on the geometries to
3774 visualize the two mapped degrees of freedom.
3776 Raises
3777 ------
3778 ValueError
3779 If an invalid mapping is specified.
3781 Returns
3782 -------
3783 output_dof_list : CoordinateArray
3784 A CoordinateArray of the same size and shape of dof_list with each
3785 entry in dof_list corresponding to the same position in output_dof_list.
3787 """
3788 dof_list = dof_list.copy()
3789 if node_map == 'position':
3790 from_nodes = np.unique(dof_list.node)
3791 to_nodes = to_geometry.node_by_global_position(self.global_node_coordinate(from_nodes)).id
3792 to_geometry = to_geometry.reduce(to_nodes)
3793 node_map = id_map(from_nodes,to_nodes)
3794 self = self.reduce(from_nodes).map_ids(node_id_map = node_map)
3795 dof_list.node = node_map(dof_list.node)
3796 elif isinstance(node_map,id_map):
3797 self = self.reduce(node_map.from_ids).map_ids(node_id_map = node_map)
3798 to_geometry = to_geometry.reduce(node_map.to_nodes)
3799 dof_list.node = node_map(dof_list.node)
3800 elif node_map is None:
3801 self = self.reduce(np.unique(dof_list.node))
3802 to_geometry = to_geometry.reduce(np.unique(dof_list.node))
3803 else:
3804 raise ValueError('Invalid node_map. Must be an id_map, "position", or None.')
3806 output_dof_list = dof_list.copy()
3807 global_directions = self.global_deflection(dof_list)
3808 to_directions = {id_num:to_geometry.global_deflection(coordinate_array(id_num,[1,2,3]))
3809 for id_num in to_geometry.node.id}
3810 for index,value in dof_list.ndenumerate():
3811 dirs = to_directions[value.node]
3812 dots = np.einsum('ij,j',dirs,global_directions[index])
3813 direction = np.argmax(abs(dots))
3814 sign = np.sign(dots[direction])
3815 output_dof_list[index].direction = sign*(direction+1)
3816 if plot_maps:
3817 plotter = self.plot_coordinate(dof_list[index])
3818 to_geometry.plot_coordinate(output_dof_list[index],plot_kwargs={'plotter':plotter})
3820 return output_dof_list
3822 def convert_elements_to_tracelines(self,keep_ids = False, start_id = None,
3823 in_place = False):
3824 """
3825 Transforms all elements in the geometry into tracelines
3827 Parameters
3828 ----------
3829 keep_ids : bool, optional
3830 If True, the element ids will be used as the new traceline ids.
3831 The default is False.
3832 start_id : int, optional
3833 The starting id number for the tracelines converted from elements.
3834 Only used if keep_ids is False. If not specified, each subsequent
3835 element will use the current highest traceline id number + 1.
3836 in_place : bool, optional
3837 If True, modify the geometry in place. If False, a copy will be
3838 created and returned. The default is False.
3840 Raises
3841 ------
3842 ValueError
3843 If conflicts are expected to occur between the newly created
3844 traceline IDs and the current existing traceline IDs.
3845 NotImplementedError
3846 If solid element types are converted.
3848 Returns
3849 -------
3850 geometry : Geometry
3851 Only returned if in_place is False.
3853 """
3854 if in_place:
3855 geometry = self
3856 else:
3857 geometry = self.copy()
3858 if keep_ids:
3859 if np.any(np.isin(geometry.element.id,geometry.traceline.id)):
3860 raise ValueError('Keeping element IDs would result in a conflict with existing traceline IDs')
3861 if start_id is not None:
3862 if np.any(np.isin(np.arange(geometry.element.size)+start_id,geometry.traceline.id)):
3863 raise ValueError('Setting new traceline ID starting at {:} would result in a conflict with existing traceline IDs'.format(start_id))
3864 if not keep_ids and start_id is None:
3865 if len(geometry.traceline) > 0:
3866 start_id = np.max(geometry.traceline.id) + 1
3867 else:
3868 start_id = 1
3869 for idx, element in enumerate(geometry.element.flatten()):
3870 if keep_ids:
3871 new_id = element.id
3872 else:
3873 new_id = start_id + idx
3874 if element.type <= 96:
3875 conn = element.connectivity
3876 conn = np.concatenate((conn,[conn[0]]))
3877 else:
3878 raise NotImplementedError('Solid Elements are not implemented yet!')
3879 geometry.add_traceline(conn,new_id,color=element.color)
3880 geometry.element = ElementArray((0,))
3881 if not in_place:
3882 return geometry
3884 def split_tracelines_into_segments(self, in_place = False):
3885 """
3886 Splits long tracelines into many length-2 tracelines
3888 Parameters
3889 ----------
3890 in_place : bool, optional
3891 If True, the current geometry will be modified in place. The
3892 default is False.
3894 Returns
3895 -------
3896 geometry : Geometry
3897 Only returned if in_place is False.
3899 """
3900 temp_geometry = Geometry()
3901 for traceline in self.traceline:
3902 connectivity = np.array((traceline.connectivity[:-1],traceline.connectivity[1:])).T
3903 connectivity = connectivity[np.all(connectivity != 0,axis=-1)]
3904 for conn in connectivity:
3905 temp_geometry.add_traceline(conn, traceline.id, traceline.description, traceline.color)
3906 if not in_place:
3907 geometry = self.copy()
3908 geometry.traceline = temp_geometry.traceline
3909 return geometry
3910 else:
3911 self.traceline = temp_geometry.traceline
3913 def remove_duplicate_tracelines(self, in_place = False):
3914 """
3915 Removes tracelines that are effectively equivalent to each other.
3917 If two tracelines have the same connectivity (or reversed connectivity),
3918 then the first will be kept and the second will be discarded. This will
3919 not take into account other fields such as description or color.
3921 Parameters
3922 ----------
3923 in_place : bool, optional
3924 If True, the current geometry will be modified in place. The
3925 default is False.
3927 Returns
3928 -------
3929 geometry : Geometry
3930 Only returned if in_place is False.
3932 """
3933 conn_set = set([])
3934 keep_array = []
3935 for traceline in self.traceline:
3936 conn = tuple(traceline.connectivity)
3937 if conn in conn_set or conn[::-1] in conn_set:
3938 keep_array.append(False)
3939 else:
3940 keep_array.append(True)
3941 conn_set.add(conn)
3942 if not in_place:
3943 geometry = self.copy()
3944 geometry.traceline = geometry.traceline[keep_array]
3945 return geometry
3946 else:
3947 self.traceline = self.traceline[keep_array]
3949 def plot(self, node_size: int = 5, line_width: int = 1, opacity=1.0, view_up=None, view_from=None, plotter=None,
3950 show_edges=False, label_nodes = False, label_tracelines = False,
3951 label_elements = False, label_font_size = 16,
3952 plot_individual_items = False):
3953 """
3954 Plots the geometry in an interactive 3D window.
3956 Parameters
3957 ----------
3958 node_size : int, optional
3959 Size to display the nodes in pixels. Set to 0 to not display nodes.
3960 The default is 5.
3961 line_width : int, optional
3962 Width to display tracelines and element edges in pixels. Set to 0
3963 to not show tracelines or edges. The default is 1.
3964 opacity : float, optional
3965 A float between 0 and 1 to specify the transparency of the geometry.
3966 Set to 1 for completely opaque, and 0 for completely transparent
3967 (invisible). The default is 1.0, no transparency.
3968 view_up : np.ndarray, optional
3969 Set the "up" direction in the plot by passing in a size-3 numpy
3970 array. The default is None.
3971 view_from : np.ndarray, optional
3972 Specify the direction from which the geometry is viewed. The
3973 default is None.
3974 plotter : BackgroundPlotter, optional
3975 A plotter can be specified to plot the geometry in an existing
3976 plot. The default is None, which creates a new window and plot.
3977 show_edges : bool, optional
3978 Specify whether or not to draw edges on elements. The default is False.
3979 label_nodes : bool or iteratble, optional
3980 Specify whether or not to label the nodes in the geometry. If True,
3981 all nodes will be plotted. If label_nodes is an array or list, it
3982 should contain the node numbers to label. The default is False,
3983 which does not label nodes.
3984 label_tracelines : bool or iteratble, optional
3985 Specify whether or not to label the tracelines in the geometry. If True,
3986 all tracelines will be plotted. If label_tracelines is an array or list, it
3987 should contain the traceline numbers to label. The default is False,
3988 which does not label tracelines.
3989 label_elements : bool or iteratble, optional
3990 Specify whether or not to label the elements in the geometry. If True,
3991 all elements will be plotted. If label_elements is an array or list, it
3992 should contain the element numbers to label. The default is False,
3993 which does not label elements.
3994 label_font_size : int, optional
3995 Specify the font size of the node labels. The default is 16.
3996 plot_individual_items : bool, optional
3997 If True, then each item will be plotted one at a time. This will
3998 make the plotting much slower, but it is required for certain
3999 features (like exporting to 3D PDF content) where having all data
4000 in one mesh is not feasible. Note that volume elements will NOT
4001 be plotted in this case.
4003 Raises
4004 ------
4005 KeyError
4006 If referenced id numbers are not found in the corresponding object,
4007 for example if a traceline references node 11 but there is no node
4008 11 in the NodeArray
4009 ValueError
4010 If an invalid or unknown element type is used
4012 Returns
4013 -------
4014 plotter : BackgroundPlotter
4015 A reference to the plotter object that the geometry was plotted in
4016 face_mesh : TYPE
4017 A reference to the mesh used to plot surface elements
4018 point_mesh : TYPE
4019 A reference to the mesh used to plot nodes and tracelines
4020 solid_mesh : TYPE
4021 A reference to the mesh used to plot volume elements
4023 """
4024 if IGNORE_PLOTS:
4025 return None
4026 # Get part information
4027 nodes = self.node.flatten()
4028 css = self.coordinate_system.flatten()
4029 elems = self.element.flatten()
4030 tls = self.traceline.flatten()
4032 num_labels = ((1 if (label_nodes is True or (label_nodes is not False and len(label_nodes) > 0)) else 0)
4033 + (1 if (label_tracelines is True or (label_tracelines is not False and len(label_tracelines) > 0)) else 0)
4034 + (1 if (label_elements is True or (label_elements is not False and len(label_elements) > 0)) else 0))
4036 coordinate_system = CoordinateSystemArray(nodes.shape)
4037 cs_map = {cs['id']:idx for idx,cs in np.ndenumerate(self.coordinate_system)}
4038 for key, node in nodes.ndenumerate():
4039 coordinate_system[key] = css[cs_map[node.def_cs]]
4040 global_node_positions = global_coord(coordinate_system, nodes.coordinate)
4041 node_index_dict = {node.id[()]: index[0] for index, node in nodes.ndenumerate()}
4042 node_index_map = np.vectorize(node_index_dict.__getitem__)
4043 tl_index_dict = {tl.id[()]: index[0] for index, tl in tls.ndenumerate()}
4044 tl_index_map = np.vectorize(tl_index_dict.__getitem__)
4045 elem_index_dict = {elem.id[()]: index[0] for index, elem in elems.ndenumerate()}
4046 elem_index_map = np.vectorize(elem_index_dict.__getitem__)
4048 # Now go through and get the element and line connectivity
4049 if not plot_individual_items:
4050 face_element_connectivity = []
4051 face_element_colors = []
4052 solid_element_connectivity = []
4053 solid_element_colors = []
4054 solid_element_types = []
4055 node_colors = []
4056 line_connectivity = []
4057 line_colors = []
4058 for index, node in nodes.ndenumerate():
4059 # element_connectivity.append(1)
4060 # element_connectivity.append(index[0])
4061 # element_colors.append(node.color)
4062 node_colors.append(node.color)
4063 for index, element in elems.ndenumerate():
4064 # Check which type of element it is
4065 if element.type in _beam_elem_types: # Beamlike element, use a line
4066 line_connectivity.append(len(element.connectivity))
4067 try:
4068 line_connectivity.extend(node_index_map(element.connectivity))
4069 except KeyError:
4070 raise KeyError(
4071 'Element {:} contains a node id not found in the node array'.format(element.id))
4072 line_colors.append(element.color)
4073 elif element.type in _face_element_types:
4074 face_element_connectivity.append(len(element.connectivity))
4075 try:
4076 face_element_connectivity.extend(node_index_map(element.connectivity))
4077 except KeyError:
4078 raise KeyError(
4079 'Element {:} contains a node id not found in the node array'.format(element.id))
4080 face_element_colors.append(element.color)
4081 elif element.type in _solid_element_types:
4082 try:
4083 solid_element_types.append(_vtk_element_map[element.type])
4084 except KeyError:
4085 raise ValueError('Do not know equivalent VTK element type for {:}: {:}'.format(
4086 element.type, _element_types[element.type]))
4087 solid_element_connectivity.append(len(element.connectivity))
4088 try:
4089 if solid_element_types[-1] in _vtk_connectivity_reorder:
4090 solid_element_connectivity.extend(node_index_map(
4091 element.connectivity[..., _vtk_connectivity_reorder[solid_element_types[-1]]]))
4092 else:
4093 solid_element_connectivity.extend(node_index_map(element.connectivity))
4094 except KeyError:
4095 raise KeyError(
4096 'Element {:} contains a node id not found in the node array'.format(element.id))
4097 solid_element_colors.append(element.color)
4098 else:
4099 raise ValueError('Unknown element type {:}'.format(element.type))
4101 for index, tl in tls.ndenumerate():
4102 for conn_group in split_list(tl.connectivity, 0):
4103 if len(conn_group) == 0:
4104 continue
4105 line_connectivity.append(len(conn_group))
4106 try:
4107 line_connectivity.extend(node_index_map(conn_group))
4108 except KeyError:
4109 raise KeyError(
4110 'Traceline {:} contains a node id not found in the node array'.format(tl.id))
4111 line_colors.append(tl.color)
4112 else:
4113 face_element_connectivity = []
4114 face_element_colors = []
4115 node_colors = []
4116 line_connectivity = []
4117 line_colors = []
4118 for index, node in nodes.ndenumerate():
4119 # element_connectivity.append(1)
4120 # element_connectivity.append(index[0])
4121 # element_colors.append(node.color)
4122 node_colors.append(node.color)
4123 for index, element in elems.ndenumerate():
4124 # Check which type of element it is
4125 if element.type in _beam_elem_types: # Beamlike element, use a line
4126 try:
4127 line_connectivity.append(node_index_map(element.connectivity))
4128 except KeyError:
4129 raise KeyError(
4130 'Element {:} contains a node id not found in the node array'.format(element.id))
4131 line_colors.append(element.color)
4132 elif element.type in _face_element_types:
4133 try:
4134 if len(element.connectivity) == 3:
4135 face_element_connectivity.append(node_index_map(element.connectivity))
4136 face_element_colors.append(element.color)
4137 else:
4138 face_element_connectivity.append(node_index_map(element.connectivity[[0,1,3]]))
4139 face_element_colors.append(element.color)
4140 face_element_connectivity.append(node_index_map(element.connectivity[[1,2,3]]))
4141 face_element_colors.append(element.color)
4142 except KeyError:
4143 raise KeyError(
4144 'Element {:} contains a node id not found in the node array'.format(element.id))
4145 elif element.type in _solid_element_types:
4146 warnings.warn('Solid Elements are currently not supported and will be skipped!')
4147 else:
4148 raise ValueError('Unknown element type {:}'.format(element.type))
4150 for index, tl in tls.ndenumerate():
4151 for conn_group in split_list(tl.connectivity, 0):
4152 if len(conn_group) == 0:
4153 continue
4154 try:
4155 mapped_conn_group = node_index_map(conn_group)
4156 except KeyError:
4157 raise KeyError(
4158 'Traceline {:} contains a node id not found in the node array'.format(tl.id))
4159 for indices in zip(mapped_conn_group[:-1],mapped_conn_group[1:]):
4160 line_connectivity.append(indices)
4161 line_colors.append(tl.color)
4163 # Now we start to plot
4164 if plotter is None:
4165 plotter = GeometryPlotter(editor=False)
4167 if not plot_individual_items:
4168 # Need to split up between point mesh, face/line mesh, and solid mesh
4169 if len(face_element_connectivity) == 0 and len(line_connectivity) == 0:
4170 face_mesh = None
4171 else:
4172 face_mesh = pv.PolyData(global_node_positions,
4173 faces=None if len(
4174 face_element_connectivity) == 0 else face_element_connectivity,
4175 lines=None if len(line_connectivity) == 0 else line_connectivity)
4176 face_mesh.cell_data['color'] = line_colors + face_element_colors
4178 plotter.add_mesh(face_mesh, scalars='color', cmap=colormap, clim=[0, 15],
4179 show_edges=show_edges, # True if line_width > 0 else False,
4180 show_scalar_bar=False, line_width=line_width,
4181 opacity=opacity)
4182 if len(solid_element_connectivity) == 0:
4183 solid_mesh = None
4184 else:
4185 solid_mesh = pv.UnstructuredGrid(np.array(solid_element_connectivity),
4186 np.array(solid_element_types, dtype='uint8'),
4187 np.array(global_node_positions))
4188 solid_mesh.cell_data['color'] = solid_element_colors
4190 plotter.add_mesh(solid_mesh, scalars='color', cmap=colormap, clim=[0, 15],
4191 show_edges=show_edges, # True if line_width > 0 else False,
4192 show_scalar_bar=False, line_width=line_width,
4193 opacity=opacity)
4195 if node_size > 0:
4196 point_mesh = pv.PolyData(global_node_positions)
4197 point_mesh.cell_data['color'] = node_colors
4198 plotter.add_mesh(point_mesh, scalars=node_colors, cmap=colormap, clim=[0, 15],
4199 show_edges=show_edges, # True if line_width > 0 else False,
4200 show_scalar_bar=False, point_size=node_size,
4201 opacity=opacity)
4202 else:
4203 point_mesh = None
4205 if (label_nodes is True or (label_nodes is not False and len(label_nodes) > 0)):
4206 try:
4207 label_node_ids = [id_num for id_num in label_nodes]
4208 label_node_indices = node_index_map(label_node_ids)
4209 except TypeError:
4210 label_node_ids = [id_num for id_num in self.node.id]
4211 label_node_indices = node_index_map(label_node_ids)
4213 node_label_mesh = pv.PolyData(global_node_positions[label_node_indices])
4214 node_label_mesh['NODE_ID'] = [
4215 ('N' if num_labels > 1 else '') + str(val) for val in label_node_ids]
4216 plotter.add_point_labels(node_label_mesh, 'NODE_ID', tolerance=0.0, shape=None,
4217 show_points=False, always_visible=True, font_size=label_font_size)
4219 if (label_tracelines is True or (label_tracelines is not False and len(label_tracelines) > 0)):
4220 try:
4221 label_traceline_ids = [id_num for id_num in label_tracelines]
4222 label_traceline_indices = tl_index_map(label_traceline_ids)
4223 except TypeError:
4224 label_traceline_ids = [id_num for id_num in self.traceline.id]
4225 label_traceline_indices = tl_index_map(label_traceline_ids)
4227 # Need to get the point to plot for each traceline
4228 traceline_locations = []
4229 for idx in label_traceline_indices:
4230 tl = tls[idx]
4231 connectivity = tl.connectivity[tl.connectivity != 0]
4232 # Find the middle of the traceline
4233 center_nodes = [connectivity[len(connectivity)//2-1],connectivity[len(connectivity)//2]]
4234 position = np.mean(global_node_positions[node_index_map(center_nodes)],axis=0)
4235 traceline_locations.append(position)
4237 tl_label_mesh = pv.PolyData(np.array(traceline_locations))
4238 tl_label_mesh['TL_ID'] = [
4239 ('L' if num_labels > 1 else '') + str(val) for val in label_traceline_ids]
4240 plotter.add_point_labels(tl_label_mesh, 'TL_ID', tolerance=0.0, shape=None,
4241 show_points=False, always_visible=True, font_size=label_font_size)
4243 if (label_elements is True or (label_elements is not False and len(label_elements) > 0)):
4244 try:
4245 label_element_ids = [id_num for id_num in label_elements]
4246 label_element_indices = elem_index_map(label_element_ids)
4247 except TypeError:
4248 label_element_ids = [id_num for id_num in self.element.id]
4249 label_element_indices = elem_index_map(label_element_ids)
4251 # Need to get the point to plot for each traceline
4252 element_locations = []
4253 for idx in label_element_indices:
4254 elem = elems[idx]
4255 connectivity = elem.connectivity
4256 # Find the middle of the traceline
4257 position = np.mean(global_node_positions[node_index_map(connectivity)],axis=0)
4258 element_locations.append(position)
4260 elem_label_mesh = pv.PolyData(np.array(element_locations))
4261 elem_label_mesh['ELEM_ID'] = [
4262 ('E' if num_labels > 1 else '') + str(val) for val in label_element_ids]
4263 plotter.add_point_labels(elem_label_mesh, 'ELEM_ID', tolerance=0.0, shape=None,
4264 show_points=False, always_visible=True, font_size=label_font_size)
4266 else:
4267 face_map = []
4268 face_mesh = []
4269 for conn,color in zip(face_element_connectivity,face_element_colors):
4270 node_indices,inverse = np.unique(conn,return_inverse=True)
4271 node_positions = global_node_positions[node_indices][inverse]
4272 node_connectivity = np.arange(node_positions.shape[0])
4273 mesh = pv.PolyData(node_positions,faces=[len(node_connectivity)]+[c for c in node_connectivity])
4274 mesh.cell_data['color'] = color
4275 mesh.SetObjectName('Elem: {:} {:} {:}'.format(*self.node.id[node_indices[inverse]]))
4276 plotter.add_mesh(mesh,scalars='color',cmap=colormap, clim = [0,15],
4277 show_edges=show_edges, show_scalar_bar = False,
4278 line_width=line_width,opacity=opacity,label='Elem: {:} {:} {:}'.format(*self.node.id[node_indices[inverse]]))
4279 face_map.append(self.node.id[node_indices[inverse]])
4280 face_mesh.append(mesh)
4282 line_map = []
4283 line_mesh = []
4284 for conn,color in zip(line_connectivity,line_colors):
4285 node_indices,inverse = np.unique(conn,return_inverse=True)
4286 node_positions = global_node_positions[node_indices][inverse]
4287 node_connectivity = np.arange(node_positions.shape[0])
4288 mesh = pv.PolyData(node_positions,lines=[len(node_connectivity)]+[c for c in node_connectivity])
4289 mesh.cell_data['color'] = color
4290 mesh.SetObjectName('Line: {:} {:}'.format(*self.node.id[node_indices[inverse]]))
4291 plotter.add_mesh(mesh,scalars='color',cmap=colormap, clim = [0,15],
4292 show_edges=show_edges, show_scalar_bar = False,
4293 line_width=line_width,opacity=opacity,label='Line: {:} {:}'.format(*self.node.id[node_indices[inverse]]))
4294 line_map.append(self.node.id[node_indices[inverse]])
4295 line_mesh.append(mesh)
4297 point_map = []
4298 point_mesh = []
4299 for node,position,color in zip(self.node.id,global_node_positions,self.node.color):
4300 if node_size > 0:
4301 mesh = pv.PolyData(position)
4302 mesh.cell_data['color'] = color
4303 plotter.add_mesh(mesh, scalars = color, cmap=colormap, clim=[0,15],
4304 show_edges=show_edges, show_scalar_bar=False, point_size=node_size,
4305 opacity=opacity)
4306 point_map.append(node)
4307 point_mesh.append(mesh)
4309 if label_nodes:
4310 try:
4311 label_node_ids = [id_num for id_num in label_nodes]
4312 label_node_indices = node_index_map(label_node_ids)
4313 except TypeError:
4314 label_node_ids = [id_num for id_num in self.node.id]
4315 label_node_indices = node_index_map(label_node_ids)
4317 for index,label in zip(label_node_indices,label_node_ids):
4318 node_label_mesh = pv.PolyData(global_node_positions[index])
4319 node_label_mesh['NODE_ID'] = [str(label)]
4320 plotter.add_point_labels(node_label_mesh, 'NODE_ID', tolerance=0.0, shape=None,
4321 show_points=False, always_visible=True, font_size=label_font_size)
4323 if view_from is not None:
4324 focus = plotter.camera.focal_point
4325 distance = plotter.camera.distance
4326 plotter.camera.position = np.array(
4327 focus) + distance * np.array(view_from) / np.linalg.norm(view_from)
4328 plotter.camera.focal_point = focus
4329 if view_up is not None:
4330 plotter.camera.up = view_up
4331 plotter.show()
4332 plotter.render()
4333 QApplication.processEvents()
4334 if plot_individual_items:
4335 return plotter, (face_mesh,face_map), (line_mesh,line_map), (point_mesh,point_map)
4336 else:
4337 return plotter, face_mesh, point_mesh, solid_mesh
4339 def global_geometry_table(self,coordinate_array_or_nodelist = None):
4340 if coordinate_array_or_nodelist is None:
4341 coord = coordinate_array(self.node.id,[1,2,3],force_broadcast=True)
4342 elif isinstance(coordinate_array_or_nodelist,CoordinateArray):
4343 coord = coordinate_array_or_nodelist.flatten()
4344 else:
4345 coord = coordinate_array(coordinate_array_or_nodelist,[1,2,3],force_broadcast=True)
4346 global_positions = self.global_node_coordinate(coord.node)
4347 global_directions = self.global_deflection(coord)
4348 coord_string = coord.string_array()
4349 df_dict = {}
4350 df_dict['Degree of Freedom'] = coord_string
4351 for label,column in zip('XYZ',global_positions.T):
4352 df_dict[f'Coordinate {label}'] = column
4353 for label,column in zip('XYZ',global_directions.T):
4354 df_dict[f'Direction {label}'] = column
4355 df = pd.DataFrame(df_dict)
4356 return df
4358 def plot_coordinate(self, coordinates: CoordinateArray = None,
4359 arrow_scale=0.1,
4360 arrow_scale_type='bbox', label_dofs=False,
4361 label_font_size=16, opacity=1.0,
4362 arrow_ends_on_node=False, plot_kwargs={},
4363 plot_individual_items = False):
4364 """
4365 Plots coordinate arrows on the geometry
4367 Parameters
4368 ----------
4369 coordinates : CoordinateArray, optional
4370 Coordinates to draw on the geometry. If no coordinates are specified,
4371 all translation degrees of freedom at each node will be plotted.
4372 arrow_scale : float | np.ndarray, optional
4373 Size of the arrows in proportion to the length of the diagonal of
4374 the bounding box of the Geometry if `arrow_scale_type` is 'bbox',
4375 otherwise the raw length of the arrow. If an ndarray is provided,
4376 it should match the shape of `coordinates`. The default is 0.1.
4377 arrow_scale_type : str, optional
4378 Specifies how to compute the size of the arrows. If 'bbox', then
4379 the arrow is scaled based on the size of the geometry. Otherwise,
4380 the arrow size is the specified length. The default is 'bbox'.
4381 label_dofs : bool, optional
4382 Specify whether or not to label the coordinates with strings.
4383 The default is False.
4384 label_font_size : int, optional
4385 Specifies the font size for the node labels.
4386 Default is 16.
4387 opacity : float, optional
4388 A float between 0 and 1 to specify the transparency of the geometry.
4389 Set to 1 for completely opaque, and 0 for completely transparent
4390 (invisible). The default is 1.0, no transparency.
4391 arrow_ends_on_node : bool, optional
4392 If True, arrow tip ends at the node, otherwise the arrow begins at
4393 node.
4394 Defualt is False
4395 plot_kwargs : dict, optional
4396 Any additional keywords that should be passed to the Geometry.plot
4397 function. The default is {}.
4398 plot_individual_items : bool, optional
4399 If True, then each item will be plotted one at a time. This will
4400 make the plotting much slower, but it is required for certain
4401 features (like exporting to 3D PDF content) where having all data
4402 in one mesh is not feasible. Note that volume elements will NOT
4403 be plotted in this case.
4405 Returns
4406 -------
4407 plotter : BackgroundPlotter
4408 A reference to the window the geometry was plotted in.
4410 """
4411 if IGNORE_PLOTS:
4412 return None
4413 plot_kwargs = plot_kwargs.copy()
4414 plot_kwargs['plot_individual_items'] = plot_individual_items
4415 plotter, face_mesh, point_mesh, solid_mesh = self.plot(opacity=opacity, **plot_kwargs)
4416 # Add deflection directions to each node
4417 # Get local deformation for each coordinate direction
4419 if coordinates is None:
4420 coordinates = from_nodelist(self.node.id)
4422 def build_coord_mesh(coordinates, geom, scale_factors=None):
4423 nodes = coordinates.node
4424 indices = np.isin(coordinates.node, geom.node.id)
4425 coordinates = coordinates[indices]
4426 nodes = nodes[indices]
4427 local_deformations = coordinates.local_direction()
4428 coordinate_systems = geom.coordinate_system(geom.node(nodes).disp_cs)
4429 points = global_coord(geom.coordinate_system(
4430 geom.node(nodes).def_cs), geom.node(nodes).coordinate)
4431 global_deflections = global_deflection(coordinate_systems, local_deformations, points)
4432 # Now add the point array to the mesh
4433 coord_mesh = pv.PolyData(points)
4434 coord_mesh.point_data['Coordinates'] = global_deflections
4435 coord_mesh.point_data['Direction'] = abs(coordinates.direction)
4436 coord_mesh.point_data['Node ID'] = nodes
4437 if scale_factors is not None:
4438 coord_mesh.point_data['Scale'] = scale_factors
4439 return coord_mesh, points, global_deflections
4441 coordinates = coordinates.flatten()
4443 if isinstance(arrow_scale, int | float):
4444 scale = False
4445 arrow_scale_array = np.ones(len(coordinates))
4446 global_arrow_scale = arrow_scale
4447 elif isinstance(arrow_scale, np.ndarray):
4448 scale = 'Scale'
4449 arrow_scale_array = arrow_scale
4450 global_arrow_scale = 1
4452 if arrow_scale_type == 'bbox':
4453 node_coords = self.global_node_coordinate()
4454 bbox_diagonal = np.linalg.norm(
4455 np.max(node_coords, axis=0) - np.min(node_coords, axis=0))
4456 arrow_factor = bbox_diagonal * global_arrow_scale
4457 else:
4458 arrow_factor = global_arrow_scale
4460 # Create Straight Arrows and Labels
4461 coordinates_straight = coordinates[np.abs(coordinates.direction) < 4]
4462 if arrow_scale_array is not None:
4463 straight_scaling = arrow_scale_array[np.abs(coordinates.direction) < 4]
4464 else:
4465 straight_scaling = None
4466 if coordinates_straight.size:
4467 if plot_individual_items:
4468 for coordinate_straight in coordinates_straight:
4469 [coord_mesh_straight, points_straight, global_deflections_straight] = build_coord_mesh(
4470 coordinate_straight[np.newaxis], self, straight_scaling)
4472 if arrow_ends_on_node:
4473 arrow_start = (-1.0, 0.0, 0.0)
4474 else:
4475 arrow_start = (0.0, 0.0, 0.0)
4476 coord_arrows_stright = coord_mesh_straight.glyph(
4477 orient='Coordinates', scale=scale, factor=arrow_factor, geom=pv.Arrow(start=arrow_start))
4478 color = ((1.0,0.0,0.0) if abs(coordinate_straight.direction) == 1 else
4479 (0.0,1.0,0.0) if abs(coordinate_straight.direction) == 2 else
4480 (0.0,0.0,1.0))
4481 plotter.add_mesh(coord_arrows_stright, color = color, show_scalar_bar=False)
4482 else:
4483 [coord_mesh_straight, points_straight, global_deflections_straight] = build_coord_mesh(
4484 coordinates_straight, self, straight_scaling)
4486 if arrow_ends_on_node:
4487 arrow_start = (-1.0, 0.0, 0.0)
4488 else:
4489 arrow_start = (0.0, 0.0, 0.0)
4490 coord_arrows_stright = coord_mesh_straight.glyph(
4491 orient='Coordinates', scale=scale, factor=arrow_factor, geom=pv.Arrow(start=arrow_start))
4492 plotter.add_mesh(coord_arrows_stright, scalars='Direction',
4493 cmap=coord_colormap, clim=[1, 6], show_scalar_bar=False)
4495 if label_dofs:
4496 if arrow_ends_on_node:
4497 dof_label_mesh = pv.PolyData(
4498 points_straight - global_deflections_straight * arrow_factor * straight_scaling[:,np.newaxis])
4499 else:
4500 dof_label_mesh = pv.PolyData(
4501 points_straight + global_deflections_straight * arrow_factor * straight_scaling[:,np.newaxis])
4502 dof_label_mesh['DOF'] = [val for val in coordinates_straight.string_array()]
4503 plotter.add_point_labels(dof_label_mesh, 'DOF', tolerance=0.0, shape=None,
4504 show_points=False, always_visible=True, font_size=label_font_size)
4506 # Create Rotational Arrows and Labels
4507 coordinates_rotation = coordinates[np.abs(coordinates.direction) > 3]
4508 if arrow_scale_array is not None:
4509 rotation_scaling = arrow_scale_array[np.abs(coordinates.direction) > 3]
4510 else:
4511 rotation_scaling = None
4512 if coordinates_rotation.size:
4513 if plot_individual_items:
4514 for coordinate_rotation in coordinates_rotation:
4515 [coord_mesh_rotation, points_rotation, global_deflections_rotation] = build_coord_mesh(
4516 coordinate_rotation[np.newaxis], self, rotation_scaling)
4518 arrow_index = np.arange(4)
4519 head_location_angles = 1 / 8 * np.pi + 1 / 2 * np.pi * arrow_index
4520 tail_location_angles = 3 / 8 * np.pi + 1 / 2 * np.pi * arrow_index
4522 r = 1
4523 arrow_head_locations = np.array([np.zeros(
4524 head_location_angles.size), r * np.sin(head_location_angles), r * np.cos(head_location_angles)]).T
4525 arrow_tail_locations = np.array([np.zeros(
4526 tail_location_angles.size), r * np.sin(tail_location_angles), r * np.cos(tail_location_angles)]).T
4527 cone_vectors = np.array([np.zeros(head_location_angles.size), np.sin(
4528 head_location_angles - np.pi / 2), np.cos(head_location_angles - np.pi / 2)]).T
4530 arc = pv.merge([pv.CircularArc(pointa=start, pointb=stop, center=(0, 0, 0)).tube(
4531 radius=0.05) for start, stop in zip(arrow_head_locations, arrow_tail_locations)])
4532 cone = pv.merge([pv.Cone(center=cone_center, direction=cone_vector, height=0.3, radius=0.1, resolution=20)
4533 for cone_center, cone_vector in zip(arrow_head_locations, cone_vectors)])
4535 curved_arrow = pv.merge([arc, cone])
4536 coord_arrows_rotation = coord_mesh_rotation.glyph(
4537 orient='Coordinates', scale=scale, factor=arrow_factor, geom=curved_arrow)
4538 color = ((1.0,0.0,0.0) if abs(coordinate_straight.direction) == 1 else
4539 (0.0,1.0,0.0) if abs(coordinate_straight.direction) == 2 else
4540 (0.0,0.0,1.0))
4541 plotter.add_mesh(coord_arrows_rotation, color=color, show_scalar_bar=False)
4542 else:
4543 [coord_mesh_rotation, points_rotation, global_deflections_rotation] = build_coord_mesh(
4544 coordinates_rotation, self, rotation_scaling)
4546 arrow_index = np.arange(4)
4547 head_location_angles = 1 / 8 * np.pi + 1 / 2 * np.pi * arrow_index
4548 tail_location_angles = 3 / 8 * np.pi + 1 / 2 * np.pi * arrow_index
4550 r = 1
4551 arrow_head_locations = np.array([np.zeros(
4552 head_location_angles.size), r * np.sin(head_location_angles), r * np.cos(head_location_angles)]).T
4553 arrow_tail_locations = np.array([np.zeros(
4554 tail_location_angles.size), r * np.sin(tail_location_angles), r * np.cos(tail_location_angles)]).T
4555 cone_vectors = np.array([np.zeros(head_location_angles.size), np.sin(
4556 head_location_angles - np.pi / 2), np.cos(head_location_angles - np.pi / 2)]).T
4558 arc = pv.merge([pv.CircularArc(pointa=start, pointb=stop, center=(0, 0, 0)).tube(
4559 radius=0.05) for start, stop in zip(arrow_head_locations, arrow_tail_locations)])
4560 cone = pv.merge([pv.Cone(center=cone_center, direction=cone_vector, height=0.3, radius=0.1, resolution=20)
4561 for cone_center, cone_vector in zip(arrow_head_locations, cone_vectors)])
4563 curved_arrow = pv.merge([arc, cone])
4564 coord_arrows_rotation = coord_mesh_rotation.glyph(
4565 orient='Coordinates', scale=scale, factor=arrow_factor, geom=curved_arrow)
4566 # Make Colors of Rotations same as Straight
4567 coord_arrows_rotation['Direction'] = coord_arrows_rotation['Direction'] - 3
4569 plotter.add_mesh(coord_arrows_rotation, scalars='Direction',
4570 cmap=coord_colormap, clim=[1, 6], show_scalar_bar=False)
4572 if label_dofs:
4573 nodes = coordinates_rotation.node
4574 indices = np.isin(coordinates_rotation.node, self.node.id)
4575 nodes = nodes[indices]
4576 local_deformations = coordinates_rotation.local_direction()
4577 local_deformation_transformation = np.array(([arrow_head_locations[0, 0], arrow_head_locations[0, 1], arrow_head_locations[0, 2]],
4578 [arrow_head_locations[0, 2], arrow_head_locations[0,
4579 0], arrow_head_locations[0, 1]],
4580 [arrow_head_locations[0, 1], arrow_head_locations[0, 2], arrow_head_locations[0, 0]]))
4581 local_deformations_new = np.transpose(
4582 local_deformation_transformation.T @ local_deformations.T)
4583 coordinate_systems = self.coordinate_system(self.node(nodes).disp_cs)
4584 points = global_coord(self.coordinate_system(
4585 self.node(nodes).def_cs), self.node(nodes).coordinate)
4586 global_deflections_rotation_new = global_deflection(
4587 coordinate_systems, local_deformations_new, points)
4589 dof_label_mesh = pv.PolyData(
4590 points_rotation + global_deflections_rotation_new * arrow_factor * rotation_scaling[:,np.newaxis])
4592 dof_label_mesh['DOF'] = [val for val in coordinates_rotation.string_array()]
4593 plotter.add_point_labels(dof_label_mesh, 'DOF', tolerance=0.0, shape=None,
4594 show_points=False, always_visible=True, font_size=label_font_size)
4595 return plotter
4597 def plot_shape(self, shape, plot_kwargs={}, background_plotter_kwargs={'editor': False}, undeformed_opacity=0.25, deformed_opacity=1.0, starting_scale=1.0):
4598 """
4599 Plot mode shapes on the geometry
4601 Parameters
4602 ----------
4603 shape : ShapeArray
4604 The set of shapes to plot
4605 plot_kwargs : dict, optional
4606 Any additional keywords that should be passed to the Geometry.plot
4607 function. The default is {}.
4608 background_plotter_kwargs : dict, optional
4609 Any additional arguments that should be passed to the
4610 BackgroundPlotter initializer. The default is {'editor':False}.
4611 undeformed_opacity : float, optional
4612 A float between 0 and 1 to specify the transparency of the undeformed
4613 geometry. Set to 1 for completely opaque, and 0 for completely
4614 transparent (invisible). The default is 0.25.
4615 deformed_opacity : float, optional
4616 A float between 0 and 1 to specify the transparency of the deformed
4617 geometry. Set to 1 for completely opaque, and 0 for completely
4618 transparent (invisible). The default is 1.0.
4619 starting_scale : float, optional
4620 The starting scale factor of the animation. The default is 1.0.
4622 Returns
4623 -------
4624 ShapePlotter
4625 A reference to the ShapePlotter class that is created to plot the
4626 animated shapes
4628 """
4629 if IGNORE_PLOTS:
4630 return None
4631 return ShapePlotter(self, shape, plot_kwargs, background_plotter_kwargs,
4632 undeformed_opacity, deformed_opacity=deformed_opacity, starting_scale=starting_scale)
4634 def plot_deflection_shape(self, deflection_shape_data, plot_kwargs={},
4635 background_plotter_kwargs={'editor': False},
4636 undeformed_opacity=0.25, deformed_opacity=1.0,
4637 starting_scale=1.0):
4638 """
4639 Plot deflection shapes shapes on the geometry
4641 Parameters
4642 ----------
4643 deflection_shape_data : NDDataArray
4644 Data array containing the deflection shapes to plot
4645 plot_kwargs : dict, optional
4646 Any additional keywords that should be passed to the Geometry.plot
4647 function. The default is {}.
4648 background_plotter_kwargs : dict, optional
4649 Any additional arguments that should be passed to the
4650 BackgroundPlotter initializer. The default is {'editor':False}.
4651 undeformed_opacity : float, optional
4652 A float between 0 and 1 to specify the transparency of the undeformed
4653 geometry. Set to 1 for completely opaque, and 0 for completely
4654 transparent (invisible). The default is 0.25.
4655 deformed_opacity : float, optional
4656 A float between 0 and 1 to specify the transparency of the deformed
4657 geometry. Set to 1 for completely opaque, and 0 for completely
4658 transparent (invisible). The default is 1.0.
4659 starting_scale : float, optional
4660 The starting scale factor of the animation. The default is 1.0.
4662 Returns
4663 -------
4664 deflection_shape_data
4665 A reference to the deflection_shape_data object that is created to
4666 plot the animated shapes
4668 """
4669 if IGNORE_PLOTS:
4670 return None
4671 return DeflectionShapePlotter(self, deflection_shape_data, plot_kwargs, background_plotter_kwargs,
4672 undeformed_opacity, deformed_opacity=deformed_opacity, starting_scale=starting_scale)
4674 def plot_transient(self, displacement_data, displacement_scale=1.0,
4675 frames_per_second=20,
4676 undeformed_opacity=0.0, deformed_opacity=1.0,
4677 plot_kwargs={},
4678 transformation_shapes=None,
4679 num_curves=50,
4680 show: bool = True,
4681 app=None,
4682 window_size=None,
4683 off_screen=None,
4684 allow_quit_keypress=True,
4685 toolbar=True,
4686 menu_bar=True,
4687 editor=False,
4688 update_app_icon=None,
4689 **kwargs):
4690 """
4691 Create a TransientPlotter object to plot displacements over time
4693 Parameters
4694 ----------
4695 displacement_data : TimeHistoryArray
4696 Transient displacement data that will be applied
4697 displacement_scale : float, optional
4698 Scale factor applied to displacements. The default is 1.0.
4699 frames_per_second : float, optional
4700 Number of time steps to plot per second while the displacement is
4701 animating. Default is 20.
4702 undeformed_opacity : float, optional
4703 Opacity of the undeformed geometry. The default is 0.0, or
4704 completely transparent.
4705 deformed_opacity : float, optional
4706 Opacity of the deformed geometry. The default is 1.0, or completely
4707 opaque.
4708 plot_kwargs : dict, optional
4709 Keyword arguments passed to the Geometry.plot function
4710 transformation_shapes : ShapeArray
4711 Shape matrix that will be used to expand the data. Must be the
4712 same size as the `displacement_data`
4713 num_curves : int, optional
4714 Maximum number of curves to plot on the time selector. Default is
4715 50.
4716 show : bool, optional
4717 Show the plotting window. If ``False``, show this window by
4718 running ``show()``. The default is True.
4719 app : QApplication, optional
4720 Creates a `QApplication` if left as `None`. The default is None.
4721 window_size : list of int, optional
4722 Window size in pixels. Defaults to ``[1024, 768]``
4723 off_screen : TYPE, optional
4724 Renders off screen when True. Useful for automated
4725 screenshots or debug testing. The default is None.
4726 allow_quit_keypress : bool, optional
4727 Allow user to exit by pressing ``"q"``. The default is True.
4728 toolbar : bool, optional
4729 If True, display the default camera toolbar. Defaults to True.
4730 menu_bar : bool, optional
4731 If True, display the default main menu. Defaults to True.
4732 editor : TYPE, optional
4733 If True, display the VTK object editor. Defaults to False.
4734 update_app_icon : bool, optional
4735 If True, update_app_icon will be called automatically to update the
4736 Qt app icon based on the current rendering output. If None, the
4737 logo of PyVista will be used. If False, no icon will be set.
4738 Defaults to None. The default is None.
4739 title : str, optional
4740 Title of plotting window.
4741 multi_samples : int, optional
4742 The number of multi-samples used to mitigate aliasing. 4 is a
4743 good default but 8 will have better results with a potential
4744 impact on performance.
4745 line_smoothing : bool, optional
4746 If True, enable line smothing
4747 point_smoothing : bool, optional
4748 If True, enable point smothing
4749 polygon_smoothing : bool, optional
4750 If True, enable polygon smothing
4751 auto_update : float, bool, optional
4752 Automatic update rate in seconds. Useful for automatically
4753 updating the render window when actors are change without
4754 being automatically ``Modified``. If set to ``True``, update
4755 rate will be 1 second.
4757 Returns
4758 -------
4759 TransientPlotter
4760 """
4761 if IGNORE_PLOTS:
4762 return None
4763 return TransientPlotter(self, displacement_data, displacement_scale,
4764 frames_per_second,
4765 undeformed_opacity, deformed_opacity,
4766 plot_kwargs,
4767 transformation_shapes,
4768 num_curves,
4769 show,
4770 app,
4771 window_size,
4772 off_screen,
4773 allow_quit_keypress,
4774 toolbar,
4775 menu_bar,
4776 editor,
4777 update_app_icon,
4778 **kwargs)
4780 def __add__(self, geometry):
4781 if not isinstance(self, geometry.__class__):
4782 return NotImplemented
4783 final_fields = {}
4784 for field in ['node', 'traceline', 'element', 'coordinate_system']:
4785 # Look through and see which are in common
4786 common_ids = np.intersect1d(getattr(self, field).id, getattr(geometry, field).id)
4787 if common_ids.size > 0:
4788 equal_ids = getattr(self, field)(common_ids) == getattr(geometry, field)(common_ids)
4789 if not all(equal_ids):
4790 raise ValueError('Both geometries contain {:} with ID {:} but they are not equivalent'.format(
4791 field, common_ids[~equal_ids]))
4792 self_ids = np.concatenate(
4793 (np.setdiff1d(getattr(self, field).id, getattr(geometry, field).id), common_ids))
4794 geometry_ids = np.setdiff1d(getattr(geometry, field).id, getattr(self, field).id)
4795 final_fields[field] = np.concatenate(
4796 (getattr(self, field)(self_ids), getattr(geometry, field)(geometry_ids)))
4797 return Geometry(**final_fields)
4799 def modify_ids(self, node_change=0, traceline_change=0, element_change=0, coordinate_system_change=0):
4800 """
4801 Shifts the id numbers in the geometry
4803 Parameters
4804 ----------
4805 node_change : int, optional
4806 The amount to shift the node ids. The default is 0.
4807 traceline_change : int, optional
4808 The amount to shift the traceline ids.. The default is 0.
4809 element_change : int, optional
4810 The amount to shift the element ids. The default is 0.
4811 coordinate_system_change : int, optional
4812 The amount to shift the coordinate system ids. The default is 0.
4814 Returns
4815 -------
4816 geom_out : Geometry
4817 A copy of the original geometry with id numbers modified
4819 """
4820 geom_out = self.copy()
4821 geom_out.node.id += node_change
4822 geom_out.traceline.id += traceline_change
4823 geom_out.element.id += element_change
4824 geom_out.coordinate_system.id += coordinate_system_change
4825 geom_out.node.def_cs += coordinate_system_change
4826 geom_out.node.disp_cs += coordinate_system_change
4827 geom_out.traceline.connectivity += node_change
4828 geom_out.element.connectivity += node_change
4829 # Set values of tracelines back to zero if they just equal node_change,
4830 # because that means they were zero and should represent picking up the
4831 # pen
4832 for connectivity_array in geom_out.traceline.connectivity:
4833 connectivity_array[connectivity_array == node_change] = 0
4834 return geom_out
4836 def map_ids(self, node_id_map=None, traceline_id_map=None, element_id_map=None, coordinate_system_id_map=None):
4837 """
4838 Maps id numbers from an original set of ids to a new set of ids.
4840 This function accepts id_map classes defining "from" and "to" ids
4841 Existing ids found in the "from" set are transformed to the corresponding
4842 id in the "to" set.
4844 Parameters
4845 ----------
4846 node_id_map : id_map, optional
4847 An id_map defining the mapping applied to node ids.
4848 The default is None, which results in the ids being unchanged
4849 traceline_id_map : id_map, optional
4850 An id_map defining the mapping applied to traceline ids.
4851 The default is None, which results in the ids being unchanged
4852 element_id_map : id_map, optional
4853 An id_map defining the mapping applied to element ids.
4854 The default is None, which results in the ids being unchanged
4855 coordinate_system_id_map : id_map, optional
4856 An id_map defining the mapping applied to coordinate system ids.
4857 The default is None, which results in the ids being unchanged
4859 Returns
4860 -------
4861 geom_out : Geometry
4862 A copy of the original geometry with id numbers modified
4864 """
4865 geom_out = self.copy()
4866 if node_id_map is not None:
4867 geom_out.node.id = node_id_map(self.node.id)
4868 for key, val in self.traceline.ndenumerate():
4869 geom_out.traceline.connectivity[key] = node_id_map(self.traceline.connectivity[key])
4870 for key, val in self.element.ndenumerate():
4871 geom_out.element.connectivity[key] = node_id_map(self.element.connectivity[key])
4872 if coordinate_system_id_map is not None:
4873 geom_out.coordinate_system.id = coordinate_system_id_map(self.coordinate_system.id)
4874 geom_out.node.disp_cs = coordinate_system_id_map(self.node.disp_cs)
4875 geom_out.node.def_cs = coordinate_system_id_map(self.node.def_cs)
4876 if traceline_id_map is not None:
4877 geom_out.traceline.id = traceline_id_map(self.traceline.id)
4878 if element_id_map is not None:
4879 geom_out.element.id = element_id_map(self.element.id)
4880 return geom_out
4882 def response_kinematic_transformation(self, response_coordinate,
4883 virtual_point_node_number,
4884 virtual_point_location=[0,0,0]):
4885 """
4886 Creates a kinematic response transformation from the rigid body body
4887 shapes for the associated geometry and response coordinates.
4889 Parameters
4890 ----------
4891 response_coordinate : CoordinateArray
4892 The "physical" response coordinates for the kinematic transformation.
4893 virtual_point_node_number : int
4894 The numeric ID for the virtual point in the kinematic transformation.
4895 virtual_point_location : list or ndarray
4896 The [X, Y, Z] coordinates that defined the location of the virtual
4897 point. The default is is [0, 0, 0].
4899 Returns
4900 -------
4901 transformation : Matrix
4902 The kinematic transformation as a matrix object. It is organized with
4903 the virtual point CoordinateArray on the rows and the physical
4904 response CoordinateArray on the columns.
4906 Notes
4907 -----
4908 The transformation automatically handles polarity differences in the geometry
4909 and response_coordinate.
4911 References
4912 ----------
4913 .. [1] M. Van der Seijs, D. van den Bosch, D. Rixen, and D. Klerk, "An improved
4914 methodology for the virtual point transformation of measured frequency
4915 response functions in dynamic substructuring," in Proceedings of the 4th
4916 International Conference on Computational Methods in Structural Dynamics
4917 and Earthquake Engineering, Kos Island, 2013, pp. 4334-4347,
4918 doi: 10.7712/120113.4816.C1539.
4919 """
4920 shapes = self.rigid_body_shapes(response_coordinate, cg = virtual_point_location)
4921 transformation_array = np.linalg.pinv(shapes.shape_matrix.T)
4923 virtual_point_coordinate = coordinate_array(node=virtual_point_node_number, direction=[1,2,3,4,5,6])
4925 return matrix(transformation_array, virtual_point_coordinate, response_coordinate)
4927 def force_kinematic_transformation(self, force_coordinate,
4928 virtual_point_node_number,
4929 virtual_point_location=[0,0,0]):
4930 """
4931 Creates a kinematic force transformation from the rigid body body
4932 shapes for the associated geometry and force coordinates.
4934 Parameters
4935 ----------
4936 force_coordinate : CoordinateArray
4937 The "physical" force coordinates for the kinematic transformation.
4938 virtual_point_node_number : int
4939 The numeric ID for the virtual point in the kinematic transformation.
4940 virtual_point_location : list or ndarray
4941 The [X, Y, Z] coordinates that defined the location of the virtual
4942 point. The default is is [0, 0, 0].
4944 Returns
4945 -------
4946 transformation : Matrix
4947 The kinematic transformation as a matrix object. It is organized with
4948 the virtual point CoordinateArray on the rows and the physical
4949 force CoordinateArray on the columns.
4951 Notes
4952 -----
4953 The transformation automatically handles polarity differences in the geometry
4954 and force_coordinate.
4956 References
4957 ----------
4958 .. [1] M. Van der Seijs, D. van den Bosch, D. Rixen, and D. Klerk, "An improved
4959 methodology for the virtual point transformation of measured frequency
4960 response functions in dynamic substructuring," in Proceedings of the 4th
4961 International Conference on Computational Methods in Structural Dynamics
4962 and Earthquake Engineering, Kos Island, 2013, pp. 4334-4347,
4963 doi: 10.7712/120113.4816.C1539.
4964 """
4965 shapes = self.rigid_body_shapes(force_coordinate, cg = virtual_point_location)
4966 virtual_point_coordinate = coordinate_array(node=virtual_point_node_number, direction=[1,2,3,4,5,6])
4968 return matrix(shapes.shape_matrix, virtual_point_coordinate, force_coordinate)
4970 def copy(self):
4971 """
4972 Return's a copy of the current Geometry
4974 Changes to the copy will not also be applied to the original geometry
4976 Returns
4977 -------
4978 Geometry
4979 A copy of the current Geometry
4981 """
4982 return Geometry(self.node.copy(), self.coordinate_system.copy(), self.traceline.copy(), self.element.copy())
4984 def save(self, filename):
4985 """Saves the geometry to a numpy .npz file
4987 The .npz file will have fields 'node', 'coordinate_system', 'element',
4988 and 'traceline', with each field storing the respective portion of the
4989 geometry
4991 Parameters
4992 ----------
4993 filename : str
4994 Filename to save the geometry to. If the filename doesn't end with
4995 .npz, it will be appended.
4997 Returns
4998 -------
4999 None.
5001 """
5002 np.savez(filename,
5003 node=self.node.view(np.ndarray),
5004 coordinate_system=self.coordinate_system.view(np.ndarray),
5005 element=self.element.view(np.ndarray),
5006 traceline=self.traceline.view(np.ndarray))
5008 def savemat(self, filename):
5009 """Saves the geometry to a .mat file
5011 The .mat file will have fields 'node', 'coordinate_system', 'element',
5012 and 'traceline', with each field storing the respective portion of the
5013 geometry
5015 Parameters
5016 ----------
5017 filename : str
5018 Filename to save the geometry to.
5020 Returns
5021 -------
5022 None.
5023 """
5025 mat_dict = {key: getattr(self, key).assemble_mat_dict() for key in
5026 ['node', 'coordinate_system', 'element', 'traceline']}
5027 scipy_savemat(filename,mat_dict)
5029 def rigid_body_shapes(self, coordinates, mass=1, inertia=np.eye(3), cg=np.zeros(3), principal_axes=False):
5030 """
5031 Creates a set of shapes corresponding to the rigid body motions
5033 Rigid body translation and rotation shapes are computed analytically
5034 from the Geoemtry
5036 Parameters
5037 ----------
5038 coordinates : CoordinateArray
5039 coordinates at which to compute deformations
5040 mass : float, optional
5041 The mass of the geometry used to scale the rigid body translation
5042 shapes. The default is 1.
5043 inertia : np.ndarray, optional
5044 A 3x3 array consisting of the mass moments of inertia, used to
5045 scale the rotation shapes. The default is np.eye(3).
5046 cg : np.ndarray, optional
5047 The center of gravity of the Geometry about which the rotations
5048 occur. The default is np.zeros(3).
5049 principal_axes : bool, optional
5050 If True, compute the principal axes of the test article and
5051 perform rotations about those axes. The default is False.
5053 Returns
5054 -------
5055 output_shape : ShapeArray
5056 A set of rigid body shapes for the current geometry
5058 """
5059 if principal_axes:
5060 rotation_scale, rotation_directions = np.linalg.eig(inertia)
5061 rotation_directions /= np.linalg.norm(rotation_directions, axis=0, keepdims=True)
5062 else:
5063 rotation_scale = np.diag(inertia)
5064 rotation_directions = np.eye(3)
5065 M = np.diag(np.concatenate((mass * np.ones(3), rotation_scale)))
5066 phi = np.sqrt(np.linalg.inv(M))
5067 # Construct the full coordinate array for all degrees of freedom for each node
5068 full_coordinates = coordinate_array(np.unique(coordinates.node)[:, np.newaxis],
5069 ['X+', 'Y+', 'Z+']).flatten()
5070 # Assume local coordinates for now, transform later
5071 translation_shape_matrix = full_coordinates.local_direction().T/np.sqrt(mass)
5072 rotation_shape_matrix = np.zeros(translation_shape_matrix.shape)
5073 for i in range(3):
5074 for j in range(len(full_coordinates)):
5075 direction = full_coordinates[j].local_direction()
5076 coord = self.global_node_coordinate(full_coordinates[j].node)
5077 rotation_shape_matrix[i, j] = phi[i + 3, i + 3] * \
5078 np.dot(direction, np.cross(rotation_directions[:, i], coord - cg))
5079 # Concatenate shapes together
5080 shape_matrix = np.concatenate((translation_shape_matrix, rotation_shape_matrix), axis=0)
5081 full_shapes = shape_array(full_coordinates, shape_matrix)
5082 # Now we need to transform to the local coordinate systems
5083 # First set up a geometry with just the global coordinate system in it
5084 global_geometry = Geometry(node=self.node.copy(),
5085 coordinate_system=coordinate_system_array((1,)))
5086 global_geometry.node.disp_cs = 1
5087 global_geometry.node.def_cs = 1
5088 transformed_shape = full_shapes.transform_coordinate_system(global_geometry, self)
5089 # Now turn it into our requested coordinates
5090 reduced_shape_matrix = transformed_shape[coordinates]
5091 output_shape = shape_array(coordinates, reduced_shape_matrix)
5092 return output_shape
5094 def coordinate_transformation_matrix(self, to_geometry : 'Geometry',
5095 nodes : np.ndarray = None,
5096 rotations : bool = False):
5097 """
5098 Creates a transformation matrix that transforms one geometry to a new
5099 geometry with different coordinate systems defined.
5101 Parameters
5102 ----------
5103 to_geometry : Geometry
5104 Geometry object with coordinate systems that the transformation
5105 matrix will transform into
5106 nodes : np.ndarray, optional
5107 An array of node id numbers to include in the transformation.
5108 The default is to include all nodes in the geometry.
5109 rotations : bool, optional
5110 If True, create degrees of freedom for rotations as well as
5111 translations.
5114 Returns
5115 -------
5116 Matrix
5117 A Matrix object that transforms data or shapes currently
5118 represented by this Geometry to the Geometry specified in the
5119 `to_geometry` argument.
5120 """
5121 if nodes is None:
5122 nodes = np.intersect1d(self.node.id,to_geometry.node.id)
5123 else:
5124 nodes = np.ravel(nodes)
5125 if rotations:
5126 dofs = coordinate_array(nodes[:,np.newaxis,np.newaxis],[[1, 2, 3],
5127 [4, 5, 6]])
5128 else:
5129 dofs = coordinate_array(nodes[:,np.newaxis],[1, 2, 3])
5130 transform_from_original = self.global_deflection(dofs)
5131 transform_to_new = to_geometry.global_deflection(dofs)
5132 # We need to transform this matrix into block-diagonal form.
5133 transform_from_original = transform_from_original.reshape(-1,3,3)
5134 transform_to_new = transform_to_new.reshape(-1,3,3)
5135 transform_from_original = block_diag(*transform_from_original)
5136 transform_to_new = block_diag(*transform_to_new)
5137 transform = transform_to_new @ np.linalg.pinv(transform_from_original)
5138 transformation_matrix = matrix(transform,dofs.flatten(),dofs.flatten())
5139 return transformation_matrix
5141 @classmethod
5142 def load(cls, filename):
5143 """
5144 Loads a geometry from a numpy .npz or .unv file
5146 The .npz file must have fields 'node', 'coordinate_system', 'element',
5147 and 'traceline', with each field storing the respective portion of the
5148 geometry
5150 The .unv file will need to have the proper datasets specified to define
5151 a geometry
5153 Parameters
5154 ----------
5155 filename : str
5156 Filename to load the geometry from.
5158 Raises
5159 ------
5160 AttributeError
5161 Raised if the calling class does not have a `from_unv` method
5162 defined
5164 Returns
5165 -------
5166 Geometry
5167 Geometry constructed from the data in the loaded file
5169 """
5170 if filename[-4:].lower() in ['.unv', '.uff']:
5171 try:
5172 from ..fileio.sdynpy_uff import readunv
5173 unv_dict = readunv(filename)
5174 return cls.from_unv(unv_dict)
5175 except AttributeError:
5176 raise AttributeError('Class {:} has no from_unv attribute defined'.format(cls))
5177 else:
5178 with np.load(filename, allow_pickle=True) as data:
5179 return Geometry(node=data['node'].view(NodeArray),
5180 coordinate_system=data['coordinate_system'].view(
5181 CoordinateSystemArray),
5182 traceline=data['traceline'].view(TracelineArray),
5183 element=data['element'].view(ElementArray))
5185 def write_to_unv(self, filename, write_nodes=True, write_coordinate_systems=True,
5186 write_tracelines=True, write_elements=True,
5187 dataset_2412_kwargs={},
5188 dataset_2420_kwargs={}):
5189 """
5190 Write the geometry to a unversal file format file
5192 Parameters
5193 ----------
5194 filename : str
5195 Filename to which the geometry will be written. If None, a
5196 unv data dictionary will be returned instead, similar to that
5197 obtained from the readunv function in sdynpy
5198 write_nodes : bool, optional
5199 If True, write the geometry's nodes to dataset 2411 in the output
5200 file. The default is True.
5201 write_coordinate_systems : True, optional
5202 If True, write the geometry's coordinate systems to dataset 2420
5203 in the output file. The default is True.
5204 write_tracelines : bool, optional
5205 If True, write the geometry's tracelines to dataset 82 in the
5206 output file. The default is True.
5207 write_elements : TYPE, optional
5208 If True, write the geometry's elements to dataset 2412 in the
5209 output file. The default is True.
5210 dataset_2412_kwargs : dict, optional
5211 Allows users to specify additional element parameters not stored
5212 by the Geometry. The default is {}.
5213 dataset_2420_kwargs : dict, optional
5214 Allows users to specify additional coordinate system parameters not
5215 stored by the Geometry. The default is {}.
5217 Returns
5218 -------
5219 unv_dict : dict
5220 Dictionary containing unv information, similar to that obtained from
5221 readunv. Only returned if filename is None.
5223 """
5224 from ..fileio.sdynpy_uff import dataset_82, dataset_2411, dataset_2412, dataset_2420
5225 # Write Coordinate systems
5226 if write_coordinate_systems:
5227 if self.coordinate_system.size > 0:
5228 cs_unv = dataset_2420.Sdynpy_UFF_Dataset_2420(
5229 cs_labels=self.coordinate_system.id.flatten(),
5230 cs_types=self.coordinate_system.cs_type.flatten(),
5231 cs_colors=self.coordinate_system.color.flatten(),
5232 cs_names=[cs.name for index, cs in self.coordinate_system.ndenumerate()],
5233 cs_matrices=self.coordinate_system.matrix.reshape(-1, 4, 3),
5234 **dataset_2420_kwargs)
5235 else:
5236 cs_unv = None
5237 if write_nodes:
5238 if self.node.size > 0:
5239 node_unv = dataset_2411.Sdynpy_UFF_Dataset_2411(
5240 self.node.id.flatten(),
5241 self.node.def_cs.flatten(),
5242 self.node.disp_cs.flatten(),
5243 self.node.color.flatten(),
5244 self.node.coordinate.reshape(-1, 3))
5245 else:
5246 node_unv = None
5247 if write_elements:
5248 if self.element.size > 0:
5249 if 'physical_property_table_numbers' not in dataset_2412_kwargs:
5250 dataset_2412_kwargs['physical_property_table_numbers'] = [1] * self.element.size
5251 if 'material_property_table_numbers' not in dataset_2412_kwargs:
5252 dataset_2412_kwargs['material_property_table_numbers'] = [1] * self.element.size
5253 if 'beam_orientations' not in dataset_2412_kwargs:
5254 dataset_2412_kwargs['beam_aft_cross_section_numbers'] = [
5255 None] * self.element.size
5256 if 'beam_aft_cross_section_numbers' not in dataset_2412_kwargs:
5257 dataset_2412_kwargs['beam_orientations'] = [None] * self.element.size
5258 if 'beam_fore_cross_section_numbers' not in dataset_2412_kwargs:
5259 dataset_2412_kwargs['beam_fore_cross_section_numbers'] = [
5260 None] * self.element.size
5261 elem_unv = dataset_2412.Sdynpy_UFF_Dataset_2412(
5262 element_labels=self.element.id.flatten(),
5263 fe_descriptor_ids=self.element.type.flatten(),
5264 colors=self.element.color.flatten(),
5265 connectivities=[elem.connectivity for index,
5266 elem in self.element.ndenumerate()],
5267 **dataset_2412_kwargs)
5268 else:
5269 elem_unv = None
5270 if write_tracelines:
5271 tl_unv = []
5272 for index, tl in self.traceline.ndenumerate():
5273 tl_unv.append(
5274 dataset_82.Sdynpy_UFF_Dataset_82(
5275 tl.id, tl.color, tl.description, tl.connectivity,
5276 ))
5277 if filename is None:
5278 unv_dict = {}
5279 unv_dict[2420] = [cs_unv]
5280 unv_dict[2411] = [node_unv]
5281 unv_dict[2412] = [elem_unv]
5282 unv_dict[82] = tl_unv
5283 return unv_dict
5284 else:
5285 with open(filename, 'w') as f:
5286 if cs_unv is not None:
5287 f.write(' -1\n')
5288 f.write(' 2420\n')
5289 f.write(cs_unv.write_string())
5290 f.write(' -1\n')
5291 if node_unv is not None:
5292 f.write(' -1\n')
5293 f.write(' 2411\n')
5294 f.write(node_unv.write_string())
5295 f.write(' -1\n')
5296 if elem_unv is not None:
5297 f.write(' -1\n')
5298 f.write(' 2412\n')
5299 f.write(elem_unv.write_string())
5300 f.write(' -1\n')
5301 for dataset in tl_unv:
5302 f.write(' -1\n')
5303 f.write(' 82\n')
5304 f.write(dataset.write_string())
5305 f.write(' -1\n')
5308from_excel_template = Geometry.from_excel_template
5309write_excel_template = Geometry.write_excel_template
5310from_imat_struct = Geometry.from_imat_struct
5311from_exodus = Geometry.from_exodus
5312from_unv = Geometry.from_unv
5313from_uff = Geometry.from_uff
5314load = Geometry.load
5317class id_map:
5318 """Class defining mapping between two sets of id numbers"""
5320 def __init__(self, from_ids, to_ids):
5321 """
5322 Initializes the id map
5324 Parameters
5325 ----------
5326 from_ids : np.ndarray
5327 Id numbers to map from
5328 to_ids : np.ndarray
5329 Id numbers to map to
5331 Returns
5332 -------
5333 None.
5335 """
5336 self.from_ids = from_ids
5337 self.to_ids = to_ids
5338 map_dict = {f: t for f, t in zip(from_ids, to_ids)}
5339 map_dict[0] = 0
5340 self.mapper = np.vectorize(map_dict.__getitem__, otypes=['uint64'])
5342 def __call__(self, val):
5343 """
5344 Map id numbers
5346 Parameters
5347 ----------
5348 val : np.ndarray
5349 A set of id numbers to map. These must all be contained in the
5350 id_map's from_ids array.
5352 Returns
5353 -------
5354 np.ndarray
5355 The set of id numbers from the id_map's to_ids array corresponding
5356 to those passed in to the function from the from_ids array.
5358 """
5359 return self.mapper(val)
5361 def inverse(self):
5362 """
5363 Returns an inverse map, swapping the from and to ids.
5365 If duplicate id entries exist, they will be overwritten.
5367 Returns
5368 -------
5369 id_map
5370 An id_map object with swapped from and to ids.
5372 """
5373 return id_map(self.to_ids,self.from_ids)