Coverage for src/sdynpy/fem/sdynpy_exodus.py: 8%
1773 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-18 21:32 +0000
1# -*- coding: utf-8 -*-
2"""
3Import data from and save data to Exodus files.
4"""
5"""
6Copyright 2022 National Technology & Engineering Solutions of Sandia,
7LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
8Government retains certain rights in this software.
10This program is free software: you can redistribute it and/or modify
11it under the terms of the GNU General Public License as published by
12the Free Software Foundation, either version 3 of the License, or
13(at your option) any later version.
15This program is distributed in the hope that it will be useful,
16but WITHOUT ANY WARRANTY; without even the implied warranty of
17MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18GNU General Public License for more details.
20You should have received a copy of the GNU General Public License
21along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""
24import netCDF4
25import numpy as np
26from types import SimpleNamespace
27import datetime
28from importlib.metadata import version
29import copy
30from ..core.sdynpy_geometry import global_coord, local_coord, global_deflection, _exodus_elem_type_map
31from ..core.sdynpy_shape import ShapeArray
32from ..core.sdynpy_data import TimeHistoryArray
34_inverse_exodus_elem_type_map = {val: key for key, val in _exodus_elem_type_map.items()}
35_inverse_exodus_elem_type_map[61] = "tri3"
36_inverse_exodus_elem_type_map[91] = "tri3"
37_inverse_exodus_elem_type_map[41] = "tri3"
38_inverse_exodus_elem_type_map[64] = "quad4"
39_inverse_exodus_elem_type_map[94] = "quad4"
40_inverse_exodus_elem_type_map[44] = "quad4"
42__version__ = '0.99'
45def face_connectivity(reduced_element_type, keep_midside_nodes):
46 if reduced_element_type.lower() == 'tetra10' and keep_midside_nodes:
47 return (('tri6', np.array([[0, 1, 3, 4, 8, 7],
48 [1, 2, 3, 5, 9, 8],
49 [2, 0, 3, 6, 7, 9],
50 [0, 2, 1, 6, 5, 4]])),)
51 elif reduced_element_type.lower() == 'tetra4' or (reduced_element_type == 'tetra10' and not keep_midside_nodes):
52 return (('tri3', np.array([[0, 1, 3],
53 [1, 2, 3],
54 [2, 0, 3],
55 [0, 2, 1]])),)
56 elif reduced_element_type.lower() == 'hex20' and keep_midside_nodes:
57 return (('quad8', np.array([[0, 1, 5, 4, 8, 13, 16, 12],
58 [1, 2, 6, 5, 9, 14, 17, 13],
59 [2, 3, 7, 6, 10, 15, 18, 14],
60 [3, 0, 4, 7, 11, 12, 19, 15],
61 [3, 2, 1, 0, 11, 8, 9, 10],
62 [4, 5, 6, 7, 16, 17, 18, 19]])),)
63 elif reduced_element_type.lower() == 'hex8' or reduced_element_type.lower() == 'hex' or (reduced_element_type == 'hex20' and not keep_midside_nodes):
64 return (('quad4', np.array([[0, 1, 5, 4],
65 [1, 2, 6, 5],
66 [2, 3, 7, 6],
67 [3, 0, 4, 7],
68 [3, 2, 1, 0],
69 [4, 5, 6, 7]])),)
70 else:
71 return ((None, None),)
74def mesh_triangulation_array(element_type):
75 if element_type == 'quad4':
76 triangles = np.array([[0, 1, 2],
77 [0, 2, 3]])
78 elif element_type == 'quad8':
79 triangles = np.array([[0, 4, 7],
80 [4, 1, 5],
81 [5, 2, 6],
82 [7, 6, 3],
83 [4, 6, 7],
84 [4, 5, 6]])
85 else:
86 triangles = None
87 return triangles
90class ExodusError(Exception):
91 pass
94class Exodus:
95 '''Read or write exodus files.
97 This class creates functionality to read or write exodus files using the
98 netCDF4 python module.
100 Parameters
101 ----------
102 filename : str
103 The path string to the file that will be opened
104 mode : str
105 Mode with which the file is opened, 'r' - Read, 'w' - Write,
106 'a' - Append. Default is 'r'.
107 title : str
108 Title of the exodus file, only required if mode='w'.
109 num_dims : int
110 Number of dimensions in the exodus file, only required if mode='w'.
111 num_nodes : int
112 Number of nodes in the exodus file, only required if mode='w'.
113 num_elem : int
114 Number of elements in the exodus file, only required if mode='w'.
115 num_blocks : int
116 Number of blocks in the exodus file, only required if mode='w'.
117 num_node_sets : int
118 Number of node sets in the exodus file, only required if mode = 'w'.
119 num_side_sets : int
120 Number of side sets in the exodus file, only required if mode='w'.
121 '''
123 def __init__(self, filename, mode='r',
124 title=None, num_dims=None, num_nodes=None, num_elem=None,
125 num_blocks=None, num_node_sets=None, num_side_sets=None,
126 clobber=False):
127 self.filename = filename
128 self.mode_character = mode
129 if mode in ['r', 'a']:
130 # Here we read in a file that already exists
131 self._ncdf_handle = netCDF4.Dataset(filename, mode)
132 else:
133 # Here we have to create one first
134 self._ncdf_handle = netCDF4.Dataset(filename, mode, clobber=clobber,
135 format='NETCDF3_64BIT_OFFSET')
136 # Assign attributes
137 self._ncdf_handle.api_version = 5.22
138 self._ncdf_handle.version = 5.22
139 self._ncdf_handle.floating_point_word_size = 8
140 self._ncdf_handle.file_size = 1
141 self._ncdf_handle.int64_status = 0
142 if title is None:
143 title = ''
144 self._ncdf_handle.title = title
146 # Assign Dimensions
147 self._ncdf_handle.createDimension('len_string', 33)
148 self._ncdf_handle.createDimension('len_name', 33)
149 self._ncdf_handle.createDimension('len_line', 81)
150 self._ncdf_handle.createDimension('four', 4)
151 if num_dims is None:
152 raise ValueError("num_dims must be assigned for mode 'w'")
153 self._ncdf_handle.createDimension('num_dim', num_dims)
154 if num_nodes is None:
155 raise ValueError("num_nodes must be assigned for mode 'w'")
156 self._ncdf_handle.createDimension('num_nodes', num_nodes)
157 if num_elem is None:
158 raise ValueError("num_elem must be assigned for mode 'w'")
159 self._ncdf_handle.createDimension('num_elem', num_elem)
160 if num_blocks is None:
161 raise ValueError("num_blocks must be assigned for mode 'w'")
162 self._ncdf_handle.createDimension('num_el_blk', num_blocks)
163 if num_node_sets is None:
164 raise ValueError("num_node_sets must be assigned for mode 'w'")
165 if num_node_sets > 0:
166 self._ncdf_handle.createDimension('num_node_sets', num_node_sets)
167 if num_side_sets is None:
168 raise ValueError("num_side_sets must be assigned for mode 'w'")
169 if num_side_sets > 0:
170 self._ncdf_handle.createDimension('num_side_sets', num_side_sets)
172 # Create convenience variables
173 if self._ncdf_handle.floating_point_word_size == 8:
174 self.floating_point_type = 'f8'
175 else:
176 self.floating_point_type = 'f4'
178 # Initialize time
179 # The time dimension is expandable in the Exodus format so we will
180 # set the size of the dimension to none, which means unlimited.
181 self._ncdf_handle.createDimension('time_step', None)
182 # Initialize the time variable
183 self._ncdf_handle.createVariable('time_whole', self.floating_point_type,
184 ('time_step',))
186 @property
187 def title(self):
188 '''Get the title of the exodus file'''
189 return self._ncdf_handle.title
191 def get_qa_records(self):
192 '''Get the quality assurance records in the exodus file
194 Returns
195 -------
196 qa_records : tuple
197 Returns a nested tuple of strings where the first index is the
198 record number, and the second is the line in the record.
200 '''
201 try:
202 raw_records = self._ncdf_handle.variables['qa_records'][:]
203 except KeyError:
204 raise ExodusError('QA Records are not yet defined!')
205 qa_records = tuple(tuple(''.join(value.decode() for value in line)
206 for line in record) for record in raw_records.data)
207 return qa_records
209 def put_qa_records(self, records):
210 '''Puts the quality assurance records in the exodus file
212 Parameters
213 ----------
214 qa_records : sequence of sequence of strings
215 A nested sequence (list/tuple/etc.) containing the quality assurance
216 records.
219 Notes
220 -----
221 Each index in qa_records should consist of a length-4 tuple of strings:
223 1. Analysis Code Name
224 2. Analysis Code Version
225 3. Analysis Date
226 4. Analysis Time
227 '''
228 if 'num_qa_rec' in self._ncdf_handle.dimensions:
229 raise ExodusError('QA Records have already been put to the exodus file')
230 # Make sure that it is the correct shape
231 num_qa_rec = len(records)
232 len_string = self._ncdf_handle.dimensions['len_string'].size
233 array = np.zeros((num_qa_rec, 4, len_string), dtype='|S1')
234 for i, record in enumerate(records):
235 if len(record) > 4:
236 raise ValueError('QA Records can only have 4 entries per record!')
237 for j, line in enumerate(record):
238 if len(line) > len_string:
239 print('Warning:When setting QA records, entry "{}" (record {} entry {}) will be truncated'.format(
240 line, i, j))
241 line = line[:len_string]
242 listed_string = [val for i, val in enumerate(line)]
243 array[i, j, :len(listed_string)] = listed_string
244 # Create the dimension
245 self._ncdf_handle.createDimension('num_qa_rec', num_qa_rec)
246 # Create the variable
247 self._ncdf_handle.createVariable('qa_records', 'S1', ('num_qa_rec', 'four', 'len_string'))
248 # Assign to it
249 self._ncdf_handle.variables['qa_records'][:] = array[:]
251 def get_info_records(self):
252 '''Get the information records in the exodus file
254 Returns
255 -------
256 info_records : tuple
257 Returns a tuple of strings where the index is the
258 record number.
259 '''
260 try:
261 raw_records = self._ncdf_handle.variables['info_records'][:]
262 except KeyError:
263 raise ExodusError('Information Records are not yet defined!')
264 info_records = tuple(''.join(value.decode() for value in line) for line in raw_records.data)
265 return info_records
267 def put_info_records(self, records):
268 '''Puts the information records in the exodus file
270 Parameters
271 ----------
272 info_records : sequence of strings
273 A sequence (list/tuple/etc.) containing the information records.
274 '''
275 if 'num_info' in self._ncdf_handle.dimensions:
276 raise ExodusError('Information Records have already been put to the exodus file')
277 # Make sure that it is the correct shape
278 num_info = len(records)
279 len_line = self._ncdf_handle.dimensions['len_line'].size
280 array = np.zeros((num_info, len_line), dtype='|S1')
281 for i, line in enumerate(records):
282 if len(line) > len_line:
283 print('Warning:When setting info records, entry "{}" (record {}) will be truncated'.format(line, i))
284 line = line[:len_line]
285 listed_string = [val for i, val in enumerate(line)]
286 array[i, :len(listed_string)] = listed_string
287 # Create the dimension
288 self._ncdf_handle.createDimension('num_info', num_info)
289 # Create the variable
290 self._ncdf_handle.createVariable('info_records', 'S1', ('num_info', 'len_line'))
291 # Assign to it
292 self._ncdf_handle.variables['info_records'][:] = array[:]
294 @property
295 def num_times(self):
296 return self._ncdf_handle.dimensions['time_step'].size
298 def get_times(self, indices=None):
299 '''Gets the time values from the exodus file
301 Returns
302 -------
303 time_array : np.masked_array
304 A masked_array containing the time values.
305 '''
306 if indices is None:
307 indices = slice(None)
308 elif np.size(indices) == 0:
309 return np.ma.zeros((0,))
310 return self._ncdf_handle.variables['time_whole'][indices]
312 def set_time(self, step, value):
313 '''Sets the time value of a given step in the exodus file
315 Parameters
316 ----------
317 step : int
318 The index in the time vector to set to value
319 value : float
320 A real number to set the value of the specified index in the time
321 vector to.
323 Notes
324 -----
325 If step is not a valid index for the time vector, the time vector will
326 be expanded so that it is.
327 '''
328 self._ncdf_handle.variables['time_whole'][step] = value
330 def set_times(self, values):
331 '''Sets the time vector for the exodus file
333 Parameters
334 ----------
335 values : array-like
336 A 1-dimensional array that has the time step values as it's entries
338 Notes
339 -----
340 If the vector is longer than the current time step vector, it will be
341 expanded. If the vector is shorter than the current time step vector,
342 only the first len(values) entries of the time vector will be assigned.
343 '''
344 self._ncdf_handle.variables['time_whole'][:] = values
346 @property
347 def num_dimensions(self):
348 return self._ncdf_handle.dimensions['num_dim'].size
350 r'''
351 _ _ _
352 | \ | | ___ __| | ___ ___
353 | \| |/ _ \ / _` |/ _ \/ __|
354 | |\ | (_) | (_| | __/\__ \
355 |_| \_|\___/ \__,_|\___||___/
356 '''
358 @property
359 def num_nodes(self):
360 return self._ncdf_handle.dimensions['num_nodes'].size
362 def get_coord_names(self):
363 '''Retrieve the coordinate names in the exodus file.
365 Returns
366 -------
367 coord_names : tuple
368 Returns a tuple of strings containing the coordinate names in the
369 exodus file.
370 '''
371 try:
372 raw_records = self._ncdf_handle.variables['coor_names']
373 except KeyError:
374 raise ExodusError('Coordinate Names are not yet defined!')
375 coord_names = tuple(''.join(value.decode() for value in line if not isinstance(
376 value, np.ma.core.MaskedConstant)) for line in raw_records)
378 return coord_names
380 def put_coord_names(self, coord_names):
381 '''Puts the coordinate names into the exodus file
383 Parameters
384 ----------
385 coord_names : sequence of strings
386 A sequence (list/tuple/etc.) containing the coordinate names.
387 '''
388 if 'coor_names' in self._ncdf_handle.variables:
389 raise ExodusError('Coordinate Names have already been put to the exodus file')
390 # Make sure that it is the correct shape
391 num_dim = self.num_dimensions # self._ncdf_handle.dimensions['num_dim'].size
392 len_name = self._ncdf_handle.dimensions['len_name'].size
393 array = np.zeros((num_dim, len_name), dtype='|S1')
394 for i, line in enumerate(coord_names):
395 if len(line) > len_name:
396 print('Warning:When setting info records, entry "{}" (record {}) will be truncated'.format(line, i))
397 line = line[:len_name]
398 listed_string = [val for i, val in enumerate(line)]
399 # Create the coordinate variables
400 array[i, :len(listed_string)] = listed_string
401 # Create the variable
402 self._ncdf_handle.createVariable('coor_names', 'S1', ('num_dim', 'len_name'))
403 # Assign to it
404 self._ncdf_handle.variables['coor_names'][:] = array[:]
406 def get_coords(self):
407 '''Retrieve the coordinates of the nodes in the exodus file.
409 Returns
410 -------
411 coords : np.array
412 Returns a 2D array with size num_dims x num_nodes
413 '''
414 # TODO Add error checking
415 coord_names = ('coordx', 'coordy', 'coordz')[:self.num_dimensions]
416 raw_list = [self._ncdf_handle.variables[name][:] for name in coord_names]
417 coords = np.array(raw_list)
418 return coords
420 def get_coord(self, index):
421 '''Retrieve the coordinates of the specified node in the exodus file.
423 Parameters
424 ----------
425 index : int
426 The global node index (not node number) of the node of which the
427 coordinates are desired
429 Returns
430 -------
431 coords : np.array
432 Returns an array with size num_dims
433 '''
434 # TODO Add Error Checking
435 coord_names = ('coordx', 'coordy', 'coordz')[:self.num_dimensions]
436 raw_list = [self._ncdf_handle.variables[name][index] for name in coord_names]
437 coords = np.array(raw_list)
438 return coords
440 def put_coords(self, coords):
441 '''Puts the coordinate values into the exodus file
443 Parameters
444 ----------
445 coords : np.ndarray
446 A 2d array containing coordinate values.
447 '''
448 coord_names = ('coordx', 'coordy', 'coordz')[:self._ncdf_handle.dimensions['num_dim'].size]
449 if coord_names[0] in self._ncdf_handle.variables:
450 raise ExodusError('Coordinates have already been put to the exodus file')
451 coords = np.array(coords)
452 if coords.shape != (self.num_dimensions, self.num_nodes):
453 raise ExodusError('coords.shape should be (self.num_dimensions,self.num_nodes)')
454 for i, name in enumerate(coord_names):
455 self._ncdf_handle.createVariable(name, self.floating_point_type, ('num_nodes',))
456 self._ncdf_handle.variables[name][:] = coords[i]
458 def get_node_num_map(self):
459 '''Retrieve the list of local node IDs from the exodus file.
461 Returns
462 -------
463 node_num_map : np.array
464 Returns a 1D array with size num_nodes, denoting the node number
465 for the node in each index
467 Notes
468 -----
469 If there is no node_num_map in the exodus file, this function simply
470 returns an array from 1 to self.num_nodes
471 '''
472 if 'node_num_map' in self._ncdf_handle.variables:
473 return self._ncdf_handle.variables['node_num_map'][:]
474 else:
475 return np.ma.MaskedArray(np.arange(self.num_nodes) + 1)
477 def put_node_num_map(self, node_num_map):
478 '''Puts a list of local node IDs into the exodus file.
480 Parameters
481 ----------
482 node_num_map : np.array
483 A 1D array with size num_nodes, denoting the node number
484 for the node in each index
486 '''
487 if 'node_num_map' in self._ncdf_handle.variables:
488 raise ExodusError('node_num_map has already been put to the exodus file.')
489 # Create the variable
490 self._ncdf_handle.createVariable('node_num_map', 'int32', ('num_nodes',))
491 self._ncdf_handle.variables['node_num_map'][:] = node_num_map
493 def get_node_variable_names(self):
494 '''Gets a tuple of nodal variable names from the exodus file
496 Returns
497 -------
498 node_var_names : tuple of strings
499 Returns a tuple containing the names of the nodal variables in the
500 exodus file.
501 '''
502 try:
503 raw_records = self._ncdf_handle.variables['name_nod_var']
504 except KeyError:
505 raise ExodusError('Node Variable Names are not defined!')
506 node_var_names = tuple(''.join(value.decode() for value in line if not isinstance(
507 value, np.ma.core.MaskedConstant)) for line in raw_records)
508 return node_var_names
510 def put_node_variable_names(self, node_var_names):
511 '''Puts the specified variable names in the exodus file
513 Parameters
514 ----------
515 node_var_names : tuple of strings
516 A tuple containing the names of the nodal variables in the model
517 '''
518 if 'name_nod_var' in self._ncdf_handle.variables:
519 raise ExodusError('Nodal Variable Names have already been put to the exodus file')
520 # Make sure that it is the correct shape
521 num_names = len(node_var_names)
522 string_length = self._ncdf_handle.dimensions['len_name'].size
523 array = np.zeros((num_names, string_length), dtype='|S1')
524 for i, line in enumerate(node_var_names):
525 if len(line) > string_length:
526 print('Warning:When setting nodal variable names, entry "{}" (record {}) will be truncated'.format(
527 line, i))
528 line = line[:string_length]
529 listed_string = [val for i, val in enumerate(line)]
530 # Create the coordinate variables
531 array[i, :len(listed_string)] = listed_string
532 self._ncdf_handle.createVariable('vals_nod_var{:d}'.format(
533 i + 1), self.floating_point_type, ('time_step', 'num_nodes'))
534 # Create the dimension
535 self._ncdf_handle.createDimension('num_nod_var', num_names)
536 # Create the variable
537 self._ncdf_handle.createVariable('name_nod_var', 'S1', ('num_nod_var', 'len_name'))
538 # Assign to it
539 self._ncdf_handle.variables['name_nod_var'][:] = array[:]
541 @property
542 def num_node_variables(self):
543 try:
544 return self._ncdf_handle.dimensions['num_nod_var'].size
545 except KeyError:
546 raise ExodusError('Number of Node Variables is not defined!')
548 def get_node_variable_values(self, name_or_index, step=None):
549 '''Gets the node variable values for the specified timestep
551 Parameters
552 ----------
553 name_or_index : str or int
554 Name or Index of the nodal variable that is desired. If
555 type(name_or_index) == str, then it is assumed to be the name. If
556 type(name_or_index) == int, then it is assumed to be the index.
557 step : int
558 Time step at which to recover the nodal variable
560 Returns
561 -------
562 node_variable_values : maskedarray
563 A 1d array consisting of the variable values for each node at the
564 specified time step.
565 '''
566 if isinstance(name_or_index, (int, np.integer)):
567 index = name_or_index
568 elif isinstance(name_or_index, (str, np.character)):
569 try:
570 index = self.get_node_variable_names().index(name_or_index)
571 except ValueError:
572 raise ExodusError('Name {} not found in self.get_node_variable_names(). Options are {}'.format(
573 name_or_index, self.get_node_variable_names()))
574 vals_nod_var_name = 'vals_nod_var{:d}'.format(index + 1)
575 if step is not None:
576 if step >= self.num_times:
577 raise ExodusError('Invalid Time Step')
578 return self._ncdf_handle.variables[vals_nod_var_name][step, :]
579 else:
580 return self._ncdf_handle.variables[vals_nod_var_name][:]
582 def get_node_variable_value(self, name_or_index, node_index, step=None):
583 '''Gets a node variable value for the specified timestep
585 Parameters
586 ----------
587 name_or_index : str or int
588 Name or Index of the nodal variable that is desired. If
589 type(name_or_index) == str, then it is assumed to be the name. If
590 type(name_or_index) == int, then it is assumed to be the index.
591 node_index : int
592 Node index at which to recover the nodal variable
593 step : int
594 Time step at which to recover the nodal variable
595 Returns
596 -------
597 node_variable_value : float
598 The variable values for the specified node at the specified time step.
599 '''
600 if isinstance(name_or_index, (int, np.integer)):
601 index = name_or_index
602 elif isinstance(name_or_index, (str, np.character)):
603 try:
604 index = self.get_node_variable_names().index(name_or_index)
605 except ValueError:
606 raise ExodusError('Name {} not found in self.get_node_variable_names(). Options are {}'.format(
607 name_or_index, self.get_node_variable_names()))
608 vals_nod_var_name = 'vals_nod_var{:d}'.format(index + 1)
609 if step is not None:
610 if step >= self.num_times:
611 raise ExodusError('Invalid Time Step')
612 return self._ncdf_handle.variables[vals_nod_var_name][step, node_index]
613 else:
614 return self._ncdf_handle.variables[vals_nod_var_name][:, node_index]
616 def set_node_variable_values(self, name_or_index, step, values):
617 '''Sets the node variable values for the specified timestep
619 Parameters
620 ----------
621 name_or_index : str or int
622 Name or Index of the nodal variable that is desired. If
623 type(name_or_index) == str, then it is assumed to be the name. If
624 type(name_or_index) == int, then it is assumed to be the index.
625 step : int
626 Time step at which to recover the nodal variable
627 values : array-like
628 A 1d array consisting of the variable values for each node at the
629 specified time step.
631 Notes
632 -----
633 If step is not a valid index for the time vector, the time vector will
634 be expanded so that it is.
635 '''
636 if isinstance(name_or_index, (int, np.integer)):
637 index = name_or_index
638 elif isinstance(name_or_index, (str, np.character)):
639 try:
640 index = self.get_node_variable_names().index(name_or_index)
641 except ValueError:
642 raise ValueError('Name {} not found in self.get_node_variable_names(). Options are {}'.format(
643 name_or_index, self.get_node_variable_names()))
644 vals_nod_var_name = 'vals_nod_var{:d}'.format(index + 1)
645 self._ncdf_handle.variables[vals_nod_var_name][step, :] = values
647 def set_node_variable_value(self, name_or_index, node_index, step, value):
648 '''Sets the node variable values for the specified timestep
650 Parameters
651 ----------
652 name_or_index : str or int
653 Name or Index of the nodal variable that is desired. If
654 type(name_or_index) == str, then it is assumed to be the name. If
655 type(name_or_index) == int, then it is assumed to be the index.
656 node_index : int
657 Node index at which to recover the nodal variable
658 step : int
659 Time step at which to recover the nodal variable
660 value : float
661 The variable value for the specified node at the specified time step.
663 Notes
664 -----
665 If step is not a valid index for the time vector, the time vector will
666 be expanded so that it is.
667 '''
668 if isinstance(name_or_index, (int, np.integer)):
669 index = name_or_index
670 elif isinstance(name_or_index, (str, np.character)):
671 try:
672 index = self.get_node_variable_names().index(name_or_index)
673 except ValueError:
674 raise ValueError('Name {} not found in self.get_node_variable_names(). Options are {}'.format(
675 name_or_index, self.get_node_variable_names()))
676 vals_nod_var_name = 'vals_nod_var{:d}'.format(index + 1)
677 self._ncdf_handle.variables[vals_nod_var_name][step, node_index] = value
679 def get_displacements(self, displacement_name='Disp', capital_coordinates=True):
680 return np.array([self.get_node_variable_values(displacement_name + (val.upper() if capital_coordinates else val.lower()))
681 for val in 'xyz'])
683 r'''
684 _____ _ _
685 | ____| | ___ _ __ ___ ___ _ __ | |_ ___
686 | _| | |/ _ \ '_ ` _ \ / _ \ '_ \| __/ __|
687 | |___| | __/ | | | | | __/ | | | |_\__ \
688 |_____|_|\___|_| |_| |_|\___|_| |_|\__|___/
689 '''
691 @property
692 def num_elems(self):
693 if 'num_elem' in self._ncdf_handle.dimensions:
694 return self._ncdf_handle.dimensions['num_elem'].size
695 else:
696 return 0
698 def get_elem_num_map(self):
699 '''Retrieve the list of local element IDs from the exodus file.
701 Returns
702 -------
703 elem_num_map : np.array
704 Returns a 1D array with size num_elems, denoting the element number
705 for the element in each index
707 Notes
708 -----
709 If there is no elem_num_map in the exodus file, this function simply
710 returns an array from 1 to self.num_elems
711 '''
712 if 'elem_num_map' in self._ncdf_handle.variables:
713 return self._ncdf_handle.variables['elem_num_map'][:]
714 else:
715 return np.ma.MaskedArray(np.arange(self.num_elems) + 1)
717 def put_elem_num_map(self, elem_num_map):
718 '''Puts a list of local element IDs into the exodus file.
720 Parameters
721 ----------
722 elem_num_map : np.array
723 A 1D array with size num_elems, denoting the element number
724 for the element in each index
726 '''
727 if 'elem_num_map' in self._ncdf_handle.variables:
728 raise ExodusError('elem_num_map has already been put to the exodus file.')
729 # Create the variable
730 self._ncdf_handle.createVariable('elem_num_map', 'int32', ('num_elem',))
731 self._ncdf_handle.variables['elem_num_map'][:] = elem_num_map
733 @property
734 def num_blks(self):
735 return self._ncdf_handle.dimensions['num_el_blk'].size
737 def get_elem_blk_ids(self):
738 ''' Gets a list of the element block ID numbers in the exodus file
740 Returns
741 -------
742 block_ids : tuple
743 A tuple containing the block ID numbers
744 '''
745 try:
746 return tuple(self._ncdf_handle.variables['eb_prop1'][:])
747 except KeyError:
748 raise ExodusError('Element Block IDs are not defined!')
750 def put_elem_blk_ids(self, block_ids):
751 '''Puts a list of the element block ID numbers in the exodus file
753 Parameters
754 ----------
755 block_ids : array-like
756 A sequency of integers specifying the block id numbers
757 '''
758 if 'eb_prop1' in self._ncdf_handle.variables:
759 raise ExodusError('Element Block IDs have already been put to the exodus file.')
760 block_ids = np.array(block_ids).flatten()
761 if len(block_ids) != self._ncdf_handle.dimensions['num_el_blk'].size:
762 raise ExodusError('Invalid Number of Block IDs Specified')
763 self._ncdf_handle.createVariable('eb_prop1', 'int32', ('num_el_blk',))
764 self._ncdf_handle.variables['eb_prop1'][:] = block_ids
765 self._ncdf_handle.variables['eb_prop1'].setncattr('name', 'ID')
766 # Create block names too
767 block_names = ['block_{:d}'.format(id) for id in block_ids]
768 num_blk, len_name = (self._ncdf_handle.dimensions['num_el_blk'].size,
769 self._ncdf_handle.dimensions['len_name'].size)
770 array = np.zeros((num_blk, len_name), dtype='|S1')
771 for i, line in enumerate(block_names):
772 if len(line) > len_name:
773 print('Warning:When setting block names, entry "{}" (record {}) will be truncated'.format(line, i))
774 line = line[:len_name]
775 listed_string = [val for i, val in enumerate(line)]
776 # Create the coordinate variables
777 array[i, :len(listed_string)] = listed_string
778 # Create the variable
779 self._ncdf_handle.createVariable('eb_names', 'S1', ('num_el_blk', 'len_name'))
780 # Assign to it
781 self._ncdf_handle.variables['eb_names'][:] = array[:]
782 # Create eb_status as well
783 self._ncdf_handle.createVariable('eb_status', 'int32', ('num_el_blk',))
784 self._ncdf_handle.variables['eb_status'][:] = 0
786 def get_elem_blk_info(self, id):
787 ''' Gets the element block information for the specified element block ID
789 Parameters
790 ----------
791 id : int
792 Element Block ID number
794 Returns
795 -------
796 element_type : str
797 Name of the Element Type
798 elements_in_block : int
799 Number of elements in the block
800 nodes_per_element : int
801 Number of nodes per element in the block
802 attributes_per_element : int
803 Number of attributes per element
804 '''
805 block_ids = self.get_elem_blk_ids()
806 try:
807 block_index = block_ids.index(id) + 1
808 except ValueError:
809 raise ExodusError('Invalid Block ID')
810 try:
811 num_el_in_blk = self._ncdf_handle.dimensions['num_el_in_blk{:d}'.format(
812 block_index)].size
813 num_nod_per_el = self._ncdf_handle.dimensions['num_nod_per_el{:d}'.format(
814 block_index)].size
815 element_type = self._ncdf_handle.variables['connect{:d}'.format(
816 block_index)].getncattr('elem_type')
817 if 'num_att_in_blk{:d}'.format(block_index) in self._ncdf_handle.dimensions:
818 attr_per_element = self._ncdf_handle.dimensions['num_att_in_blk{:d}'.format(
819 block_index)].size
820 else:
821 attr_per_element = 0
822 except KeyError:
823 raise ExodusError('Block {:d} not initialized correctly')
824 return element_type, num_el_in_blk, num_nod_per_el, attr_per_element
826 def put_elem_blk_info(self, id, elem_type, num_elements, num_nodes_per_element, num_attrs_per_elem=0):
827 '''Puts the element block information for an element block ID into the exodus file
829 Parameters
830 ----------
831 id : int
832 The block ID (not index) that the information is being specified for
833 elem_type : str
834 The element type ('SHELL4','HEX8', etc.) in the block
835 num_elements : int
836 The number of elements in the block
837 num_nodes_per_element : int
838 The number of nodes per element in the block
839 num_attrs_per_element : int
840 The number of attributes per element in the block
841 '''
842 block_ids = self.get_elem_blk_ids()
843 try:
844 block_index = block_ids.index(id) + 1
845 except ValueError:
846 raise ExodusError('Invalid Block ID')
847 num_el_in_blk_name = 'num_el_in_blk{:d}'.format(block_index)
848 num_nod_per_el_name = 'num_nod_per_el{:d}'.format(block_index)
849 if num_el_in_blk_name in self._ncdf_handle.dimensions:
850 raise ExodusError(
851 'Information for block {:d} has already been put to the exodus file.'.format(id))
852 connect_name = 'connect{:d}'.format(block_index)
853 # Create the dimensions
854 self._ncdf_handle.createDimension(num_el_in_blk_name, num_elements)
855 self._ncdf_handle.createDimension(num_nod_per_el_name, num_nodes_per_element)
856 # Create the connectivity matrix
857 self._ncdf_handle.createVariable(
858 connect_name, 'int32', (num_el_in_blk_name, num_nod_per_el_name))
859 self._ncdf_handle.variables[connect_name].setncattr('elem_type', elem_type)
860 # Create attributes if necessary
861 if num_attrs_per_elem > 0:
862 print('Warning:Element attribute functionality has not been thoroughly checked out, use with caution.')
863 num_att_in_blk_name = 'num_att_in_blk{:d}'.format(block_index)
864 attrib_name = 'attrib{:d}'.format(block_index)
865 self._ncdf_handle.createDimension(num_att_in_blk_name, num_attrs_per_elem)
866 self._ncdf_handle.createVariable(
867 attrib_name, self.floating_point_type, (num_el_in_blk_name, num_att_in_blk_name))
868 # Set the eb_status to 1 when it is defined (I'm not sure this is
869 # actually what it should be, just every exodus file I've looked at has
870 # the eb_status as all ones...)
871 self._ncdf_handle.variables['eb_status'][block_index - 1] = 1
873 def get_elem_connectivity(self, id):
874 '''Gets the element connectivity matrix for a given block
876 Parameters
877 ----------
878 id : int
879 The block id for which the element connectivity matrix is desired.
881 Returns
882 -------
883 connectivity : np.masked_array
884 The 2d connectivity matrix for the block with dimensions num_elem x
885 num_nodes_per_element. The returned value is converted from the
886 1-based Exodus indexing to 0-based Python/NumPy indexing.
887 '''
888 block_ids = self.get_elem_blk_ids()
889 try:
890 block_index = block_ids.index(id) + 1
891 except ValueError:
892 raise ExodusError('Invalid Block ID')
893 conn_name = 'connect{:d}'.format(block_index)
894 if not self.num_elems == 0:
895 try:
896 return self._ncdf_handle.variables[conn_name][:] - 1
897 except KeyError:
898 raise ExodusError('Element Block {:d} has not been defined yet')
900 def set_elem_connectivity(self, id, connectivity):
901 '''Sets the element connectivity matrix for a given block
903 Parameters
904 ----------
905 id : int
906 The block id for which the element connectivity is being assigned
907 connectivity : 2D array-like
908 A 2D array of dimension num_elements x num_nodes_per_element
909 defining the element connectivity. Note that the connectivity
910 matrix should be 0-based (first node index is zero) per Python
911 conventions. It will be converted to one-based when it is written
912 to the Exodus file.
914 '''
915 block_ids = self.get_elem_blk_ids()
916 try:
917 block_index = block_ids.index(id) + 1
918 except ValueError:
919 raise ExodusError('Invalid Block ID')
920 conn_name = 'connect{:d}'.format(block_index)
921 if conn_name not in self._ncdf_handle.variables:
922 raise ExodusError('Element Block {:d} has not been defined yet')
923 connectivity = np.array(connectivity)
924 if self._ncdf_handle.variables[conn_name][:].shape != connectivity.shape:
925 raise ExodusError('Shape of connectivity matrix should be {} (recieved {})'.format(self._ncdf_handle.variables[conn_name][:].shape,
926 connectivity.shape))
927 self._ncdf_handle.variables[conn_name][:] = connectivity + 1
929 def num_attr(self, id):
930 '''Gets the number of attributes per element for an element block
932 Parameters
933 ----------
934 d : int
935 The block id for which the number of attributes is desired
937 Returns
938 -------
939 attr_per_element : int
940 The number of attributes per element in a block
941 '''
942 block_ids = self.get_elem_blk_ids()
943 try:
944 block_index = block_ids.index(id) + 1
945 except ValueError:
946 raise ExodusError('Invalid Block ID')
947 if 'num_att_in_blk{:d}'.format(block_index) in self._ncdf_handle.dimensions:
948 attr_per_element = self._ncdf_handle.dimensions['num_att_in_blk{:d}'.format(
949 block_index)].size
950 else:
951 attr_per_element = 0
952 return attr_per_element
954 def num_elems_in_blk(self, id):
955 '''Gets the number of elements in an element block
957 Parameters
958 ----------
959 d : int
960 The block id for which the number of attributes is desired
962 Returns
963 -------
964 elem_per_element : int
965 The number of elements in the block
966 '''
967 block_ids = self.get_elem_blk_ids()
968 try:
969 block_index = block_ids.index(id) + 1
970 except ValueError:
971 raise ExodusError('Invalid Block ID')
972 try:
973 return self._ncdf_handle.dimensions['num_el_in_blk{:d}'.format(block_index)].size
974 except KeyError:
975 raise ExodusError('Block {:d} not initialized correctly'.format(id))
977 def num_nodes_per_elem(self, id):
978 '''Gets the number of nodes per element in an element block
980 Parameters
981 ----------
982 d : int
983 The block id for which the number of attributes is desired
985 Returns
986 -------
987 nodes_per_elem : int
988 The number of nodes per element in the block
989 '''
990 block_ids = self.get_elem_blk_ids()
991 try:
992 block_index = block_ids.index(id) + 1
993 except ValueError:
994 raise ExodusError('Invalid Block ID')
995 try:
996 return self._ncdf_handle.dimensions['num_nod_per_el{:d}'.format(block_index)].size
997 except KeyError:
998 raise ExodusError('Block {:d} not initialized correctly'.format(id))
1000 def get_elem_attr(self, id):
1001 '''Gets the element attributes for a given block
1003 Parameters
1004 ----------
1005 id : int
1006 The block id for which the element connectivity matrix is desired.
1008 Returns
1009 -------
1010 attributes : np.masked_array
1011 The 2d attribute matrix for the block with dimensions num_elem x
1012 num_attributes
1013 '''
1014 block_ids = self.get_elem_blk_ids()
1015 try:
1016 block_index = block_ids.index(id) + 1
1017 except ValueError:
1018 raise ExodusError('Invalid Block ID')
1019 try:
1020 return self._ncdf_handle.variables['attrib{:d}'.format(block_index)][:]
1021 except KeyError:
1022 raise ExodusError('No attributes defined for Block {:d}'.format(id))
1024 def set_elem_attr(self, id, attributes):
1025 '''Sets the element attributes for a given block
1027 Parameters
1028 ----------
1029 id : int
1030 The block id for which the element connectivity matrix is desired.
1031 attributes : 2D array-like
1032 The 2d attribute matrix for the block with dimensions num_elem x
1033 num_attributes
1034 '''
1035 block_ids = self.get_elem_blk_ids()
1036 try:
1037 block_index = block_ids.index(id) + 1
1038 except ValueError:
1039 raise ExodusError('Invalid Block ID')
1040 attr_name = 'attrib{:d}'.format(block_index)
1041 if attr_name not in self._ncdf_handle.variables:
1042 raise ExodusError('Element Block {:d} has no defined attributes')
1043 attributes = np.array(attributes)
1044 if self._ncdf_handle.variables[attr_name][:].shape != attributes.shape:
1045 raise ExodusError('Shape of attributes matrix should be {} (recieved {})'.format(self._ncdf_handle.variables[attr_name][:].shape,
1046 attributes.shape))
1047 self._ncdf_handle.variables[attr_name][:] = attributes
1049 def get_elem_type(self, id):
1050 '''Gets the element type for a given block
1052 Parameters
1053 ----------
1054 id : int
1055 The block id for which the element connectivity matrix is desired.
1057 Returns
1058 -------
1059 type : str
1060 The element type of the block
1061 '''
1062 block_ids = self.get_elem_blk_ids()
1063 try:
1064 block_index = block_ids.index(id) + 1
1065 except ValueError:
1066 raise ExodusError('Invalid Block ID')
1067 try:
1068 return self._ncdf_handle.variables['connect{:d}'.format(block_index)].getncattr('elem_type')
1069 except KeyError:
1070 raise ExodusError('Element Block {:d} has not been defined yet')
1072 def get_elem_variable_names(self):
1073 '''Gets a tuple of element variable names from the exodus file
1075 Returns
1076 -------
1077 elem_var_names : tuple of strings
1078 Returns a tuple containing the names of the element variables in the
1079 exodus file.
1080 '''
1081 try:
1082 raw_records = self._ncdf_handle.variables['name_elem_var']
1083 except KeyError:
1084 raise ExodusError('Element Variable Names are not defined!')
1085 elem_var_names = tuple(''.join(value.decode() for value in line if not isinstance(
1086 value, np.ma.core.MaskedConstant)) for line in raw_records)
1087 return elem_var_names
1089 def put_elem_variable_names(self, elem_var_names, elem_var_table=None):
1090 '''Puts the specified variable names in the exodus file
1092 Parameters
1093 ----------
1094 elem_var_names : tuple of strings
1095 A tuple containing the names of the element variables in the exodus
1096 file
1097 elem_var_table : 2d array-like
1098 A 2d array of shape num_el_blk x num_elem_var defining which
1099 element variables are defined for which element blocks
1100 elem_var_table[i,j] should be True if the j-th element variable is
1101 defined for the i-th block (index, not id number) in the model.
1102 If not specified, it is assumed that all variables are defined for
1103 all blocks.
1104 '''
1105 if 'name_elem_var' in self._ncdf_handle.variables:
1106 raise ExodusError('Element Variable Names have already been put to the exodus file')
1107 # Make sure that it is the correct shape
1108 num_names = len(elem_var_names)
1109 if elem_var_table is None:
1110 elem_var_table = np.ones((self.num_blks, num_names), dtype=bool)
1111 elem_var_table = np.array(elem_var_table)
1112 if elem_var_table.shape != (self.num_blks, num_names):
1113 raise ExodusError('Shape of elem_var_table matrix should be {} (recieved {})'.format((self.num_blks, num_names),
1114 elem_var_table.shape))
1115 string_length = self._ncdf_handle.dimensions['len_name'].size
1116 array = np.zeros((num_names, string_length), dtype='|S1')
1117 for i, line in enumerate(elem_var_names):
1118 if len(line) > string_length:
1119 print('Warning:When setting nodal variable names, entry "{}" (record {}) will be truncated'.format(
1120 line, i))
1121 line = line[:string_length]
1122 listed_string = [val for i, val in enumerate(line)]
1123 # Create the coordinate variables
1124 array[i, :len(listed_string)] = listed_string
1125 # Create the element variables if they exist
1126 for j, blkid in enumerate(self.get_elem_blk_ids()):
1127 if elem_var_table[j, i]:
1128 self._ncdf_handle.createVariable('vals_elem_var{:d}eb{:d}'.format(
1129 i + 1, j + 1), self.floating_point_type, ('time_step', 'num_el_in_blk{:d}'.format(j + 1)))
1130 # Create the dimension
1131 self._ncdf_handle.createDimension('num_elem_var', num_names)
1132 # Create the variable
1133 self._ncdf_handle.createVariable('name_elem_var', 'S1', ('num_elem_var', 'len_name'))
1134 self._ncdf_handle.createVariable('elem_var_tab', 'int32', ('num_el_blk', 'num_elem_var'))
1135 # Assign to it
1136 self._ncdf_handle.variables['name_elem_var'][:] = array[:]
1137 self._ncdf_handle.variables['elem_var_tab'][:] = elem_var_table
1139 def get_elem_variable_table(self):
1140 '''Gets the element variable table
1142 Gets the element variable table showing which elements are defined for
1143 which blocks.
1145 Returns
1146 -------
1147 elem_var_table : np.masked_array
1148 A 2D array with dimension num_blocks x num_element_variables. If
1149 the jth element variable is defined for the ith block,
1150 elem_var_table[i,j] == True
1151 '''
1152 try:
1153 return self._ncdf_handle.variables['elem_var_tab'][:]
1154 except KeyError:
1155 raise ExodusError('Element Variable Table has not been defined in the exodus file')
1157 @property
1158 def num_elem_variables(self):
1159 try:
1160 return self._ncdf_handle.dimensions['num_elem_var'].size
1161 except KeyError:
1162 raise ExodusError('Number of element variables is not defined')
1164 def get_elem_variable_values(self, block_id, name_or_index, step=None):
1165 '''Gets a block's element variable values for the specified timestep
1167 Parameters
1168 ----------
1169 block_id : int
1170 Block id number for the block from which element variable values
1171 are desired.
1172 name_or_index : str or int
1173 Name or Index of the element variable that is desired. If
1174 type(name_or_index) == str, then it is assumed to be the name. If
1175 type(name_or_index) == int, then it is assumed to be the index.
1176 step : int
1177 Time step at which to recover the element variable
1179 Returns
1180 -------
1181 elem_variable_values : maskedarray
1182 A 1d array consisting of the variable values for each element in
1183 the specified block at the specified time step.
1184 '''
1185 block_ids = self.get_elem_blk_ids()
1186 try:
1187 block_index = block_ids.index(block_id)
1188 except ValueError:
1189 raise ExodusError('Invalid Block ID')
1190 var_names = self.get_elem_variable_names()
1191 if isinstance(name_or_index, (int, np.integer)):
1192 var_index = name_or_index
1193 var_name = var_names[var_index]
1194 elif isinstance(name_or_index, (str, np.character)):
1195 try:
1196 var_name = name_or_index
1197 var_index = var_names.index(var_name)
1198 except ValueError:
1199 raise ExodusError(
1200 'Name {} not found in self.get_elem_variable_names(). Options are {}'.format(var_name, var_names))
1201 if not self.get_elem_variable_table()[block_index, var_index]:
1202 raise ExodusError('Variable {} is not defined for block {}'.format(var_name, block_id))
1203 variable_value_name = 'vals_elem_var{:d}eb{:d}'.format(var_index + 1, block_index + 1)
1204 if step is not None:
1205 if step >= self.num_times:
1206 raise ExodusError('Invalid Timestep')
1207 return self._ncdf_handle.variables[variable_value_name][step, :]
1208 else:
1209 return self._ncdf_handle.variables[variable_value_name][:, :]
1211 def get_elem_variable_value(self, block_id, name_or_index, element_index, step=None):
1212 '''Gets an element's variable value for the specified timestep
1214 Parameters
1215 ----------
1216 block_id : int
1217 Block id number for the block from which element variable value
1218 is desired.
1219 name_or_index : str or int
1220 Name or Index of the element variable that is desired. If
1221 type(name_or_index) == str, then it is assumed to be the name. If
1222 type(name_or_index) == int, then it is assumed to be the index.
1223 element_index : int
1224 element index at which to recover the nodal variable
1225 step : int
1226 Time step at which to recover the nodal variable
1227 Returns
1228 -------
1229 elem_variable_value : float
1230 The variable values for the specified element at the specified time
1231 step.
1232 '''
1233 block_ids = self.get_elem_blk_ids()
1234 try:
1235 block_index = block_ids.index(block_id)
1236 except ValueError:
1237 raise ExodusError('Invalid Block ID')
1238 var_names = self.get_elem_variable_names()
1239 if isinstance(name_or_index, (int, np.integer)):
1240 var_index = name_or_index
1241 var_name = var_names[var_index]
1242 elif isinstance(name_or_index, (str, np.character)):
1243 try:
1244 var_name = name_or_index
1245 var_index = var_names.index(var_name)
1246 except ValueError:
1247 raise ExodusError(
1248 'Name {} not found in self.get_elem_variable_names(). Options are {}'.format(var_name, var_names))
1249 if not self.get_elem_variable_table()[block_index, var_index]:
1250 raise ExodusError('Variable {} is not defined for block {}'.format(var_name, block_id))
1251 variable_value_name = 'vals_elem_var{:d}eb{:d}'.format(var_index + 1, block_index + 1)
1252 if step is not None:
1253 if step >= self.num_times:
1254 raise ExodusError('Invalid Timestep')
1255 return self._ncdf_handle.variables[variable_value_name][step, element_index]
1256 else:
1257 return self._ncdf_handle.variables[variable_value_name][:, element_index]
1259 def set_elem_variable_values(self, block_id, name_or_index, step, values):
1260 '''Sets a block's element variable values for the specified timestep
1262 Parameters
1263 ----------
1264 block_id : int
1265 Block id number for the block from which element variable values
1266 are desired.
1267 name_or_index : str or int
1268 Name or Index of the element variable that is desired. If
1269 type(name_or_index) == str, then it is assumed to be the name. If
1270 type(name_or_index) == int, then it is assumed to be the index.
1271 step : int
1272 Time step at which to recover the element variable
1273 values : array-like
1274 A 1d array consisting of the variable values for each element in
1275 the specified block at the specified time step.
1277 Notes
1278 -----
1279 If step is not a valid index for the time vector, the time vector will
1280 be expanded so that it is.
1281 '''
1282 block_ids = self.get_elem_blk_ids()
1283 try:
1284 block_index = block_ids.index(block_id)
1285 except ValueError:
1286 raise ExodusError('Invalid Block ID')
1287 var_names = self.get_elem_variable_names()
1288 if isinstance(name_or_index, (int, np.integer)):
1289 var_index = name_or_index
1290 var_name = var_names[var_index]
1291 elif isinstance(name_or_index, (str, np.character)):
1292 try:
1293 var_name = name_or_index
1294 var_index = var_names.index(var_name)
1295 except ValueError:
1296 raise ExodusError(
1297 'Name {} not found in self.get_elem_variable_names(). Options are {}'.format(var_name, var_names))
1298 if not self.get_elem_variable_table()[block_index, var_index]:
1299 raise ExodusError('Variable {} is not defined for block {}'.format(var_name, block_id))
1300 variable_value_name = 'vals_elem_var{:d}eb{:d}'.format(var_index + 1, block_index + 1)
1301# values = np.array(values).flatten()
1302# if values.shape != (self.num_elems_in_blk(block_id),):
1303# raise ExodusError('values should have {:d} elements'.format(self.num_elems_in_blk(block_id)))
1304 self._ncdf_handle.variables[variable_value_name][step, :] = values
1306 def set_elem_variable_value(self, block_id, name_or_index, element_index,
1307 step, value):
1308 '''Sets an element variable value for the specified timestep
1310 Parameters
1311 ----------
1312 block_id : int
1313 Block id number for the block from which element variable values
1314 are desired.
1315 name_or_index : str or int
1316 Name or Index of the element variable that is desired. If
1317 type(name_or_index) == str, then it is assumed to be the name. If
1318 type(name_or_index) == int, then it is assumed to be the index.
1319 step : int
1320 Time step at which to recover the element variable
1321 value : float
1322 The variable values for the specified element in the specified
1323 block at the specified time step.
1325 Notes
1326 -----
1327 If step is not a valid index for the time vector, the time vector will
1328 be expanded so that it is.
1329 '''
1330 block_ids = self.get_elem_blk_ids()
1331 try:
1332 block_index = block_ids.index(block_id)
1333 except ValueError:
1334 raise ExodusError('Invalid Block ID')
1335 var_names = self.get_elem_variable_names()
1336 if isinstance(name_or_index, (int, np.integer)):
1337 var_index = name_or_index
1338 var_name = var_names[var_index]
1339 elif isinstance(name_or_index, (str, np.character)):
1340 try:
1341 var_name = name_or_index
1342 var_index = var_names.index(var_name)
1343 except ValueError:
1344 raise ExodusError(
1345 'Name {} not found in self.get_elem_variable_names(). Options are {}'.format(var_name, var_names))
1346 if not self.get_elem_variable_table()[block_index, var_index]:
1347 raise ExodusError('Variable {} is not defined for block {}'.format(var_name, block_id))
1348 variable_value_name = 'vals_elem_var{:d}eb{:d}'.format(var_index + 1, block_index + 1)
1349 self._ncdf_handle.variables[variable_value_name][step, element_index] = value
1351 def get_element_property_names(self):
1352 raise NotImplementedError
1354 def get_element_property_value(self):
1355 raise NotImplementedError
1357 def put_element_property_names(self):
1358 raise NotImplementedError
1360 def put_element_property_value(self):
1361 raise NotImplementedError
1362 r'''
1363 _ _ _ _
1364 | \ | | ___ __| | ___ ___ ___| |_ ___
1365 | \| |/ _ \ / _` |/ _ \/ __|/ _ \ __/ __|
1366 | |\ | (_) | (_| | __/\__ \ __/ |_\__ \
1367 |_| \_|\___/ \__,_|\___||___/\___|\__|___/
1368 '''
1370 @property
1371 def num_node_sets(self):
1372 return self._ncdf_handle.dimensions['num_node_sets'].size
1374 def get_node_set_names(self):
1375 '''Retrieve the node set names in the exodus file.
1377 Returns
1378 -------
1379 ns_names : tuple
1380 Returns a tuple of strings containing the node set names in the
1381 exodus file.
1382 '''
1383 try:
1384 raw_records = self._ncdf_handle.variables['ns_names']
1385 except KeyError:
1386 raise ExodusError('Node set names are not yet defined!')
1387 ns_names = tuple(''.join(value.decode() for value in line if not isinstance(
1388 value, np.ma.core.MaskedConstant)) for line in raw_records)
1389 return ns_names
1391 def put_node_set_names(self, ns_names):
1392 '''Puts the node set names into the exodus file
1394 Parameters
1395 ----------
1396 ns_names : sequence of strings
1397 A sequence (list/tuple/etc.) containing the node set names.
1398 '''
1399 if 'ns_names' in self._ncdf_handle.variables:
1400 raise ExodusError('Node set names have already been put to the exodus file')
1401 # Make sure that it is the correct shape
1402 num_ns = self.num_node_sets # self._ncdf_handle.dimensions['num_dim'].size
1403 if num_ns != len(ns_names):
1404 raise ExodusError('Length of ns_names should be same as self.num_node_sets')
1405 len_name = self._ncdf_handle.dimensions['len_name'].size
1406 array = np.zeros((num_ns, len_name), dtype='|S1')
1407 for i, line in enumerate(ns_names):
1408 if len(line) > len_name:
1409 print('Warning:When setting info records, entry "{}" (record {}) will be truncated'.format(line, i))
1410 line = line[:len_name]
1411 listed_string = [val for i, val in enumerate(line)]
1412 # Create the coordinate variables
1413 array[i, :len(listed_string)] = listed_string
1414 # Create the variable
1415 self._ncdf_handle.createVariable('ns_names', 'S1', ('num_node_sets', 'len_name'))
1416 # Assign to it
1417 self._ncdf_handle.variables['ns_names'][:] = array[:]
1419 def get_node_set_ids(self):
1420 ''' Gets a list of the node set ID numbers in the exodus file
1422 Returns
1423 -------
1424 ns_ids : tuple
1425 A tuple containing the node set ID numbers
1426 '''
1427 try:
1428 return tuple(self._ncdf_handle.variables['ns_prop1'][:])
1429 except KeyError:
1430 raise ExodusError('Node set IDs are not defined!')
1432 def put_node_set_ids(self, ns_ids):
1433 '''Puts a list of the node set ID numbers in the exodus file
1435 Parameters
1436 ----------
1437 ns_ids : array-like
1438 A sequency of integers specifying the node set id numbers
1439 '''
1440 if 'ns_prop1' in self._ncdf_handle.variables:
1441 raise ExodusError('Node set IDs have already been put to the exodus file.')
1442 ns_ids = np.array(ns_ids).flatten()
1443 if len(ns_ids) != self._ncdf_handle.dimensions['num_node_sets'].size:
1444 raise ExodusError('Invalid number of Node set IDs specified')
1445 self._ncdf_handle.createVariable('ns_prop1', 'int32', ('num_node_sets',))
1446 self._ncdf_handle.variables['ns_prop1'][:] = ns_ids
1447 self._ncdf_handle.variables['ns_prop1'].setncattr('name', 'ID')
1448 # Create ns_status as well
1449 self._ncdf_handle.createVariable('ns_status', 'int32', ('num_node_sets',))
1450 self._ncdf_handle.variables['ns_status'][:] = 0
1452 def get_node_set_num_nodes(self, id):
1453 '''Get the number of nodes in the specified node set
1455 Parameters
1456 ----------
1457 id : int
1458 Node set ID (not index)
1460 Returns
1461 -------
1462 num_nodes : int
1463 The number of nodes in the node set.
1464 '''
1465 ns_ids = self.get_node_set_ids()
1466 try:
1467 ns_index = ns_ids.index(id) + 1
1468 except ValueError:
1469 raise ExodusError('Invalid node set ID')
1470 try:
1471 return self._ncdf_handle.dimensions['num_nod_ns{:d}'.format(ns_index)].size
1472 except KeyError:
1473 raise ExodusError('Node set {:d} not initialized correctly'.format(id))
1475 def get_node_set_dist_factors(self, id):
1476 '''Get the distribution factors of the specified node set
1478 Parameters
1479 ----------
1480 id : int
1481 Node set ID (not index)
1483 Returns
1484 -------
1485 dist_factor : np.array
1486 The distribution factors in the node set
1487 '''
1488 ns_ids = self.get_node_set_ids()
1489 try:
1490 ns_index = ns_ids.index(id) + 1
1491 except ValueError:
1492 raise ExodusError('Invalid node set ID')
1493 try:
1494 return self._ncdf_handle.variables['dist_fact_ns{:d}'.format(ns_index)][:]
1495 except KeyError:
1496 raise ExodusError('No distribution factors defined for node set {:d}'.format(id))
1498 def get_node_set_nodes(self, id):
1499 '''Get the nodes in the specified node set
1501 Parameters
1502 ----------
1503 id : int
1504 Node set ID (not index)
1506 Returns
1507 -------
1508 nodes : np.array
1509 The node indices of nodes in the node set. Note that while the
1510 Exodus file format uses 1-based indexing, the returned nodes array
1511 is converted so it is 0-based.
1512 '''
1513 ns_ids = self.get_node_set_ids()
1514 try:
1515 ns_index = ns_ids.index(id) + 1
1516# print(ns_index)
1517 except ValueError:
1518 raise ExodusError('Invalid node set ID')
1519 try:
1520 return self._ncdf_handle.variables['node_ns{:d}'.format(ns_index)][:] - 1
1521 except KeyError:
1522 raise ExodusError('Node set {:d} not initialized correctly'.format(id))
1524 def put_node_set_info(self, id, nodes, dist_fact=None):
1525 '''Puts the node set information for a node set ID into the exodus file
1527 Parameters
1528 ----------
1529 id : int
1530 The node set ID (not index)
1531 nodes : array-like
1532 A 1d array containing the node indices. Note the array should be
1533 zero-based, it will be converted to one based when written to the
1534 exodus file.
1535 dist_fact : array-like
1536 A 1d array containing the node set distribution factors. If not
1537 specified, the distribution factors will be assumed to be 1.
1538 '''
1539 ns_ids = self.get_node_set_ids()
1540 try:
1541 ns_index = ns_ids.index(id) + 1
1542 except ValueError:
1543 raise ExodusError('Invalid node set ID')
1544 num_node_in_ns_name = 'num_nod_ns{:d}'.format(ns_index)
1545 if num_node_in_ns_name in self._ncdf_handle.dimensions:
1546 raise ExodusError(
1547 'Information for node set {:d} has already been put to the exodus file.'.format(id))
1548 nodes = np.array(nodes).flatten()
1549 if dist_fact is not None:
1550 dist_fact = np.array(dist_fact).flatten()
1551 if nodes.shape != dist_fact.shape:
1552 raise ExodusError('dist_fact must be the same size as nodes')
1553 # Create dimension
1554 self._ncdf_handle.createDimension(num_node_in_ns_name, len(nodes))
1555 # Create variables
1556 node_name = 'node_ns{:d}'.format(ns_index)
1557 self._ncdf_handle.createVariable(node_name, 'int32', (num_node_in_ns_name,))
1558 self._ncdf_handle.variables[node_name][:] = nodes + 1
1559 if dist_fact is not None:
1560 distfact_name = 'dist_fact_ns{:d}'.format(ns_index)
1561 self._ncdf_handle.createVariable(
1562 distfact_name, self.floating_point_type, (num_node_in_ns_name,))
1563 self._ncdf_handle.variables[distfact_name][:] = dist_fact
1564 # Set the ns_status to 1 when it is defined (I'm not sure this is
1565 # actually what it should be, just every exodus file I've looked at has
1566 # the eb_status as all ones...)
1567 self._ncdf_handle.variables['ns_status'][ns_index - 1] = 1
1568 r'''
1569 _____ _ _ _
1570 / ____(_) | | | |
1571 | (___ _ __| | ___ ___ ___| |_ ___
1572 \___ \| |/ _` |/ _ \/ __|/ _ \ __/ __|
1573 ____) | | (_| | __/\__ \ __/ |_\__ \
1574 |_____/|_|\__,_|\___||___/\___|\__|___/
1575 '''
1577 @property
1578 def num_side_sets(self):
1579 return self._ncdf_handle.dimensions['num_side_sets'].size
1581 def get_side_set_names(self):
1582 '''Retrieve the side set names in the exouds file.
1584 Returns
1585 -------
1586 ss_names : tuple
1587 Returns a tuple of strings containing the side set names in the
1588 exodus file
1589 '''
1590 try:
1591 raw_records = self._ncdf_handle.variables['ss_names']
1592 except KeyError:
1593 raise ExodusError('Side set names are not yet defined!')
1594 ss_names = tuple(''.join(value.decode() for value in line if not isinstance(
1595 value, np.ma.core.MaskedConstant)) for line in raw_records)
1596 return ss_names
1598 def put_side_set_names(self, ss_names):
1599 '''Puts the side set names into the exodus file
1601 Parameters
1602 ----------
1603 ss_names : sequence of strings
1604 A sequence (list/tuple/etc.) containing the side set names.
1605 '''
1606 if 'xs_names' in self._ncdf_handle.variables:
1607 raise ExodusError('Side set names have already been put to the exodus file')
1608 # Make sure that it is the correct shape
1609 num_ss = self.num_side_sets
1610 if num_ss != len(ss_names):
1611 raise ExodusError('Length of ss_names should be same as self.num_side_sets')
1612 len_name = self._ncdf_handle.dimensions['len_name'].size
1613 array = np.zeros((num_ss, len_name), dtype='|S1')
1614 for i, line in enumerate(ss_names):
1615 if len(line) > len_name:
1616 print('Warning:When setting info records, entry "{}" (record {}) will be truncated'.format(line, i))
1617 line = line[:len_name]
1618 listed_string = [val for i, val in enumerate(line)]
1619 # Create the coordinate variables
1620 array[i, :len(listed_string)] = listed_string
1621 # Create the variable
1622 self._ncdf_handle.createVariable('ss_names', 'S1', ('num_side_sets', 'len_name'))
1623 # Assign to it
1624 self._ncdf_handle.variables['ss_names'][:] = array[:]
1626 def get_side_set_ids(self):
1627 ''' Gets a list of the side set ID numbers in the exodus file
1629 Returns
1630 -------
1631 ss_ids : tuple
1632 A tuple containing the side set ID numbers
1633 '''
1634 try:
1635 return tuple(self._ncdf_handle.variables['ss_prop1'][:])
1636 except KeyError:
1637 raise ExodusError('Side set IDs are not defined!')
1639 def put_side_set_ids(self, ss_ids):
1640 '''Puts a list of the side set ID numbers in the exodus file
1642 Parameters
1643 ----------
1644 ss_ids : array-like
1645 A sequency of integers specifying the side set id numbers
1646 '''
1647 if 'ss_prop1' in self._ncdf_handle.variables:
1648 raise ExodusError('Side set IDs have already been put to the exodus file.')
1649 ss_ids = np.array(ss_ids).flatten()
1650 if len(ss_ids) != self._ncdf_handle.dimensions['num_side_sets'].size:
1651 raise ExodusError('Invalid number of side set IDs specified')
1652 self._ncdf_handle.createVariable('ss_prop1', 'int32', ('num_side_sets',))
1653 self._ncdf_handle.variables['ss_prop1'][:] = ss_ids
1654 self._ncdf_handle.variables['ss_prop1'].setncattr('name', 'ID')
1655 # Create ns_status as well
1656 self._ncdf_handle.createVariable('ss_status', 'int32', ('num_side_sets',))
1657 self._ncdf_handle.variables['ss_status'][:] = 0
1659 def get_side_set_num_faces(self, id):
1660 '''Get the number of faces in the specified side set
1662 Parameters
1663 ----------
1664 id : int
1665 Side set ID (not index)
1667 Returns
1668 -------
1669 num_faces : int
1670 The number of faces in the side set.
1671 '''
1672 ss_ids = self.get_side_set_ids()
1673 try:
1674 ss_index = ss_ids.index(id) + 1
1675 except ValueError:
1676 raise ExodusError('Invalid side set ID')
1677 try:
1678 return self._ncdf_handle.dimensions['num_side_ss{:d}'.format(ss_index)].size
1679 except KeyError:
1680 raise ExodusError('Side set {:d} not initialized correctly'.format(id))
1682 def get_side_set_dist_factors(self, id):
1683 '''Get the distribution factors of the specified side set
1685 Parameters
1686 ----------
1687 id : int
1688 Side set ID (not index)
1690 Returns
1691 -------
1692 dist_factor : np.array
1693 The distribution factors in the node set
1694 '''
1695 ss_ids = self.get_side_set_ids()
1696 try:
1697 ss_index = ss_ids.index(id) + 1
1698 except ValueError:
1699 raise ExodusError('Invalid side set ID')
1700 try:
1701 return self._ncdf_handle.variables['dist_fact_ss{:d}'.format(ss_index)][:]
1702 except KeyError:
1703 raise ExodusError('No distribution factors defined for side set {:d}'.format(id))
1705 def get_side_set_faces(self, id):
1706 '''Get the faces in the specified side set
1708 Parameters
1709 ----------
1710 id : int
1711 Side set ID (not index)
1713 Returns
1714 -------
1715 elements : np.array
1716 The element of each face in the sideset (converted from 1-based
1717 exodus indexing to 0-based python indexing)
1718 sides : np.array
1719 The element side of each face in the sideset (converted from
1720 1-based exodus indexing to 0-based python indexing)
1721 '''
1722 ss_ids = self.get_side_set_ids()
1723 try:
1724 ss_index = ss_ids.index(id) + 1
1725# print(ss_index)
1726 except ValueError:
1727 raise ExodusError('Invalid side set ID')
1728 try:
1729 return (self._ncdf_handle.variables['elem_ss{:d}'.format(ss_index)][:] - 1,
1730 self._ncdf_handle.variables['side_ss{:d}'.format(ss_index)][:] - 1)
1731 except KeyError:
1732 raise ExodusError('Side set {:d} not initialized correctly'.format(id))
1734 def put_side_set_info(self, id, elements, sides, dist_fact=None):
1735 '''Puts the side set information for a side set ID into the exodus file
1737 Parameters
1738 ----------
1739 id : int
1740 The side set ID (not index)
1741 elements : np.array
1742 The element of each face in the sideset (converted from 1-based
1743 exodus indexing to 0-based python indexing)
1744 sides : np.array
1745 The element side of each face in the sideset (converted from
1746 1-based exodus indexing to 0-based python indexing)
1747 dist_fact : array-like
1748 A 1d array containing the node set distribution factors. If not
1749 specified, the distribution factors will be assumed to be 1.
1750 '''
1751 ss_ids = self.get_side_set_ids()
1752 try:
1753 ss_index = ss_ids.index(id) + 1
1754 except ValueError:
1755 raise ExodusError('Invalid side set ID')
1756 num_side_in_ss_name = 'num_side_ss{:d}'.format(ss_index)
1757 if num_side_in_ss_name in self._ncdf_handle.dimensions:
1758 raise ExodusError(
1759 'Information for side set {:d} has already been put to the exodus file.'.format(id))
1760 elements = np.array(elements).flatten()
1761 sides = np.array(sides).flatten()
1762 if dist_fact is not None:
1763 print('Warning:Distribution factors for sidesets have not been thoroughly checked! Use with Caution!')
1764 dist_fact = np.array(dist_fact).flatten()
1765 # TODO: Check if the size of dist_fact is correct
1766# if nodes.shape != dist_fact.shape:
1767# raise ExodusError('dist_fact must be the same size as nodes')
1768 # Create dimensions
1769 self._ncdf_handle.createDimension(num_side_in_ss_name, len(elements))
1770 # Create variables
1771 elem_name = 'elem_ss{:d}'.format(ss_index)
1772 side_name = 'side_ss{:d}'.format(ss_index)
1773 self._ncdf_handle.createVariable(elem_name, 'int32', (num_side_in_ss_name,))
1774 self._ncdf_handle.createVariable(side_name, 'int32', (num_side_in_ss_name,))
1775 self._ncdf_handle.variables[elem_name][:] = elements + 1
1776 self._ncdf_handle.variables[side_name][:] = sides + 1
1777 if dist_fact is not None:
1778 distfact_name = 'dist_fact_ss{:d}'.format(ss_index)
1779 distfact_dim_name = 'num_df_ss{:d}'.format(ss_index)
1780 self._ncdf_handle.createDimension(distfact_dim_name, len(dist_fact))
1781 self._ncdf_handle.createVariable(
1782 distfact_name, self.floating_point_type, (distfact_dim_name,))
1783 self._ncdf_handle.variables[distfact_name][:] = dist_fact
1784 # Set the ns_status to 1 when it is defined (I'm not sure this is
1785 # actually what it should be, just every exodus file I've looked at has
1786 # the eb_status as all ones...)
1787 self._ncdf_handle.variables['ss_status'][ss_index - 1] = 1
1789 def get_side_set_node_list(self, id):
1790 raise NotImplementedError
1792 r'''
1793 _____ _ _ _ __ __ _ _ _
1794 / ____| | | | | | \ \ / / (_) | | | |
1795 | | __| | ___ | |__ __ _| | \ \ / /_ _ _ __ _ __ _| |__ | | ___ ___
1796 | | |_ | |/ _ \| '_ \ / _` | | \ \/ / _` | '__| |/ _` | '_ \| |/ _ \/ __|
1797 | |__| | | (_) | |_) | (_| | | \ / (_| | | | | (_| | |_) | | __/\__ \
1798 \_____|_|\___/|_.__/ \__,_|_| \/ \__,_|_| |_|\__,_|_.__/|_|\___||___/
1799 '''
1801 @property
1802 def num_global_variables(self):
1803 try:
1804 return self._ncdf_handle.dimensions['num_glo_var'].size
1805 except KeyError:
1806 raise ExodusError('Number of global variables is not defined')
1808 def put_global_variable_names(self, global_var_names):
1809 '''Puts the specified global variable names in the exodus file
1811 Parameters
1812 ----------
1813 global_var_names : tuple of strings
1814 A tuple containing the names of the global variables in the model
1815 '''
1816 if 'name_glo_var' in self._ncdf_handle.variables:
1817 raise ExodusError('Global Variable Names have already been put to the exodus file')
1818 # Make sure that it is the correct shape
1819 num_names = len(global_var_names)
1820 string_length = self._ncdf_handle.dimensions['len_name'].size
1821 array = np.zeros((num_names, string_length), dtype='|S1')
1822 for i, line in enumerate(global_var_names):
1823 if len(line) > string_length:
1824 print('Warning:When setting global variable names, entry "{}" (record {}) will be truncated'.format(
1825 line, i))
1826 line = line[:string_length]
1827 listed_string = [val for i, val in enumerate(line)]
1828 # Create the coordinate variables
1829 array[i, :len(listed_string)] = listed_string
1830 # Create the dimension
1831 self._ncdf_handle.createDimension('num_glo_var', num_names)
1832 # Create the variable
1833 self._ncdf_handle.createVariable('name_glo_var', 'S1', ('num_glo_var', 'len_name'))
1834 self._ncdf_handle.createVariable(
1835 'vals_glo_var', self.floating_point_type, ('time_step', 'num_glo_var'))
1836 # Assign to it
1837 self._ncdf_handle.variables['name_glo_var'][:] = array[:]
1839 def get_global_variable_names(self):
1840 '''Gets a tuple of global variable names from the exodus file
1842 Returns
1843 -------
1844 global_var_names : tuple of strings
1845 Returns a tuple containing the names of the global variables in the
1846 exodus file.
1847 '''
1848 try:
1849 raw_records = self._ncdf_handle.variables['name_glo_var']
1850 except KeyError:
1851 raise ExodusError('Global Variable Names are not defined!')
1852 global_var_names = tuple(''.join(value.decode() for value in line if not isinstance(
1853 value, np.ma.core.MaskedConstant)) for line in raw_records)
1854 return global_var_names
1856 def get_global_variable_values(self, name_or_index, step=None):
1857 '''Gets the global variable value for the specified timesteps
1859 Parameters
1860 ----------
1861 name_or_index : str or int
1862 Name or Index of the global variable that is desired. If
1863 type(name_or_index) == str, then it is assumed to be the name. If
1864 type(name_or_index) == int, then it is assumed to be the index.
1865 step : int
1866 Time step at which to recover the nodal variable
1868 Returns
1869 -------
1870 global_variable_values : maskedarray or float
1871 A 1d array or floating point number consisting of the variable
1872 values at the specified time steps.
1873 '''
1874 if isinstance(name_or_index, (int, np.integer)):
1875 index = name_or_index
1876 elif isinstance(name_or_index, (str, np.character)):
1877 try:
1878 index = self.get_global_variable_names().index(name_or_index)
1879 except ValueError:
1880 raise ExodusError('Name {} not found in self.get_global_variable_names(). Options are {}'.format(
1881 name_or_index, self.get_node_variable_names()))
1882 vals_glo_var_name = 'vals_glo_var'
1883 if step is not None:
1884 if step >= self.num_times:
1885 raise ExodusError('Invalid Time Step')
1886 return self._ncdf_handle.variables[vals_glo_var_name][step, index]
1887 else:
1888 return self._ncdf_handle.variables[vals_glo_var_name][:, index]
1890 def set_global_variable_values(self, name_or_index, step, value):
1891 '''Sets the global variable values for the specified timestep
1893 Parameters
1894 ----------
1895 name_or_index : str or int
1896 Name or Index of the nodal variable that is desired. If
1897 type(name_or_index) == str, then it is assumed to be the name. If
1898 type(name_or_index) == int, then it is assumed to be the index.
1899 step : int
1900 Time step at which to recover the nodal variable
1901 value : array-like
1902 A 1d array consisting of the variable values for each node at the
1903 specified time step.
1905 Notes
1906 -----
1907 If step is not a valid index for the time vector, the time vector will
1908 be expanded so that it is.
1909 '''
1910 if isinstance(name_or_index, (int, np.integer)):
1911 index = name_or_index
1912 elif isinstance(name_or_index, (str, np.character)):
1913 try:
1914 index = self.get_global_variable_names().index(name_or_index)
1915 except ValueError:
1916 raise ExodusError('Name {} not found in self.get_global_variable_names(). Options are {}'.format(
1917 name_or_index, self.get_node_variable_names()))
1918 vals_glo_var_name = 'vals_glo_var'
1919 self._ncdf_handle.variables[vals_glo_var_name][step, index] = value
1921 # TODO: Still need to do Nodeset and Sideset Variables;
1922 # Coordinate Frames (see ExodusIO.m), Element Block Attribute Names
1924 def close(self):
1925 self._ncdf_handle.close()
1927 def load_into_memory(self, close=True, variables=None, timesteps=None, blocks=None):
1928 '''Loads the exodus file into an ExodusInMemory object
1930 This function loads the exodus file into memory in an ExodusInMemory
1931 format. Not for use with large files.
1933 Parameters
1934 ----------
1935 close : bool
1936 Close the netcdf file upon loading into memory. Optional argument,
1937 default is true.
1938 variables : iterable
1939 A list of variable names that are loaded into memory. Default is
1940 to load all variables
1941 timesteps : iterable
1942 A list of timestep indices that are loaded into memory. Default is
1943 to load all timesteps.
1944 blocks : iterable
1945 A list of block ids that are loaded into memory. Default is to
1946 load all blocks.
1948 Returns
1949 -------
1950 fexo : ExodusInMemory
1951 The exodus file in an ExodusInMemory format
1953 '''
1954 fexo = ExodusInMemory(self, variables=variables, timesteps=timesteps, blocks=blocks)
1955 if close:
1956 self.close()
1957 return fexo
1959 def get_block_surface(self, block_id, keep_midside_nodes=False, warn=True):
1960 '''Gets the node indices and element connectivity of surface elements
1962 This function "skins" the element block, returning a list of node
1963 indices and a surface connectivity matrix.
1965 Parameters
1966 ----------
1967 block_id : int
1968 The ID number of the block that will be skinned.
1969 keep_midside_nodes : bool
1970 Specifies whether or not to keep midside nodes in the surface mesh.
1971 Default is False.
1972 warn : bool
1973 Specifies whether or not to warn the user if the block ID doesn't
1974 have a skinning method defined for its element type. Default is
1975 True.
1977 Returns
1978 -------
1979 element_block_information : list
1980 A list of tuples of element information. These data are
1981 element_type, node_indices, block_surface_connectivity, and
1982 block_surface_original_elements. The element_type is a string
1983 representing the new block element type ('quad4','tri3',etc.). The
1984 node_indices can be used as an index into the coordinate or nodal
1985 variable arrays to select nodes corresponding to this block. The
1986 block_surface_connectivity represents the connectivity array of the
1987 surface faces of the block. Values in this array correspond to
1988 indices into the node_indices array. To recover the connectivity
1989 array in the original node indices of the exodus file, it can be
1990 passed through the node_indices as
1991 node_indices[block_surface_connectivity]. The
1992 block_surface_original_elements array shows the original element
1993 indices of the block that each surface came from. This can be
1994 used to map element variables to the new surface mesh. This list
1995 will normally be length 1 unless an element type is processed that
1996 has two different surface elements in it (e.g. wedges have
1997 tris and quads)
1998 '''
1999 elem_type = self.get_elem_type(block_id)
2000 connectivity = self.get_elem_connectivity(block_id)
2001 face_connectivity_array = face_connectivity(elem_type, keep_midside_nodes)
2002 element_block_information = []
2003 for reduced_elem_type, reduced_connectivity in face_connectivity_array:
2004 if reduced_elem_type is None:
2005 if warn:
2006 print('Warning: Element Type {:} has no defined face reduction. Passing connectivity unchanged.'.format(
2007 elem_type))
2008 block_face_connectivity = connectivity
2009 block_elem_type = elem_type
2010 block_face_original_elements = np.arange(connectivity.shape[0])
2011 else:
2012 block_face_connectivity = connectivity[:,
2013 reduced_connectivity].reshape(-1, reduced_connectivity.shape[-1])
2014 block_face_original_elements = np.repeat(np.arange(connectivity.shape[0]),
2015 reduced_connectivity.shape[0])
2016 block_elem_type = reduced_elem_type
2017 # Remove faces that are duplicates
2018 (unique_rows, unique_row_indices, unique_inverse, unique_counts) = np.unique(np.sort(
2019 block_face_connectivity, axis=1), return_index=True, return_inverse=True, return_counts=True, axis=0)
2020 original_unique_counts = unique_counts[unique_inverse]
2021 nondupe_faces = original_unique_counts == 1
2022 block_face_connectivity = block_face_connectivity[nondupe_faces, :]
2023 block_face_original_elements = block_face_original_elements[nondupe_faces]
2025 node_indices = np.unique(block_face_connectivity)
2026 node_map = np.zeros(self.num_nodes, dtype=int)
2027 node_map[node_indices] = np.arange(len(node_indices))
2028 block_face_connectivity = node_map[block_face_connectivity]
2029 element_block_information.append((block_elem_type,
2030 node_indices,
2031 block_face_connectivity,
2032 block_face_original_elements))
2033 return element_block_information
2035 def triangulate_surface_mesh(self):
2036 '''Triangulate a surface mesh for plotting patches
2038 This function generates a triangle mesh for each block in the model if
2039 it can. If there are more than 3 nodes per element in a block, and the
2040 triangulation scheme hasn't been defined in
2041 pyexodus.mesh_triangulation_array, it will be skipped.
2043 Parameters
2044 ----------
2045 None
2047 Returns
2048 -------
2049 triangulated_mesh_info : list
2050 A list of tuples containing block id, node_indices, triangulated
2051 connectivity, and original block elements
2052 '''
2053 triangulated_mesh_info = []
2054 for block_id in self.get_elem_blk_ids():
2055 surface_mesh_info = self.get_block_surface(block_id, warn=False)
2056 for elem_type, node_indices, connectivity, original_elements in surface_mesh_info:
2057 triangulation_scheme = mesh_triangulation_array(elem_type)
2058 if triangulation_scheme is None:
2059 if connectivity.shape[-1] > 3:
2060 print('Warning: More than 3 ({:}) nodes per element in block {:}, this block will not be triangulated'.format(
2061 connectivity.shape[-1], block_id))
2062 else:
2063 triangulated_mesh_info.append(
2064 (block_id, node_indices, connectivity, original_elements))
2065 else:
2066 triangulated_mesh_info.append(
2067 (block_id,
2068 node_indices,
2069 connectivity[:,
2070 triangulation_scheme].reshape(-1, triangulation_scheme.shape[-1]),
2071 np.repeat(original_elements, triangulation_scheme.shape[0])))
2072 return triangulated_mesh_info
2074 def reduce_to_surfaces(self, *args, **kwargs):
2075 return reduce_exodus_to_surfaces(self, *args, **kwargs)
2077 def extract_sharp_edges(self, *args, **kwargs):
2078 return extract_sharp_edges(self, *args, **kwargs)
2080 def copy_structure(self, filename, variables = None, clobber = False, close = False):
2081 """
2082 Copies the structure of an exodus file into a new file.
2084 Copies the number of nodes/elements, element blocks, variable names, etc.
2085 However it does not copy the timesteps or any variable values.
2087 Parameters
2088 ----------
2089 filename : str
2090 The name of the new file into which the data will be copied.
2091 variables : list of str
2092 A list of variable names to copy over. If not specified, all
2093 variables will be copied.
2094 clobber : bool
2095 If True, if the file already exists, it will be overwritten.
2096 If False (default), an error will occur if the file already
2097 exists.
2098 close : bool
2099 If True, the new file will be closed after it is created.
2100 If False (default), the new file will remain open for further
2101 modifications.
2102 """
2103 # Create the empty exodus file
2104 exo = Exodus(filename, mode='w', clobber=clobber,
2105 title = self.title,
2106 num_dims = self.num_dimensions,
2107 num_nodes = self.num_nodes,
2108 num_elem = self.num_elems,
2109 num_blocks = self.num_blks,
2110 num_node_sets = self.num_node_sets,
2111 num_side_sets = self.num_side_sets)
2113 # Add QA Record for this operation
2114 qa_records = self.get_qa_records()
2115 qa_records = qa_records + (('sdynpy Exodus.repack', version('sdynpy'),str(datetime.datetime.now().date()),str(datetime.datetime.now().time())),)
2116 exo.put_qa_records(qa_records)
2118 # Copy over the info records
2119 exo.put_info_records(self.get_info_records())
2121 # Set up the nodes
2122 exo.put_coord_names(self.get_coord_names())
2123 exo.put_coords(self.get_coords())
2124 if 'node_num_map' in self._ncdf_handle.variables:
2125 exo.put_node_num_map(self.get_node_num_map())
2127 # Set up nodal variables
2128 if variables is None:
2129 names = self.get_node_variable_names()
2130 else:
2131 names = [name for name in self.get_node_variable_names() if name in variables]
2132 if len(names) > 0:
2133 exo.put_node_variable_names(names)
2135 # Set up elements
2136 if 'elem_num_map' in self._ncdf_handle.variables:
2137 exo.put_elem_num_map(self.get_elem_num_map())
2139 exo.put_elem_blk_ids(self.get_elem_blk_ids())
2141 for block in self.get_elem_blk_ids():
2142 elem_type, elems_in_block, nodes_per_element, attributes_per_element = self.get_elem_blk_info(block)
2143 exo.put_elem_blk_info(block, elem_type, elems_in_block, nodes_per_element, attributes_per_element)
2144 exo.set_elem_connectivity(block,self.get_elem_connectivity(block))
2145 if attributes_per_element > 0:
2146 exo.set_elem_attr(block, self.get_elem_attr(block))
2148 if variables is None:
2149 elem_variable_names = self.get_elem_variable_names()
2150 elem_variable_table = self.get_elem_variable_table()
2151 else:
2152 elem_variable_names = self.get_elem_variable_names()
2153 elem_variable_table = self.get_elem_variable_table()
2154 is_kept = [name in variables for name in elem_variable_names]
2155 elem_variable_table = elem_variable_table[:,is_kept]
2156 elem_variable_names = [name for name in elem_variable_names if name in variables]
2157 if len(elem_variable_names) > 0:
2158 exo.put_elem_variable_names(elem_variable_names,elem_variable_table)
2160 # Set up nodesets
2161 if self.num_node_sets > 0:
2162 exo.put_node_set_names(self.get_node_set_names())
2163 exo.put_node_set_ids(self.get_node_set_ids())
2165 for node_set in self.get_node_set_ids():
2166 nodes = self.get_node_set_nodes(node_set)
2167 try:
2168 dist_factor = self.get_node_set_dist_factors(node_set)
2169 except sdpy.fem.sdynpy_exodus.ExodusError:
2170 dist_factor = None
2171 exo.put_node_set_info(node_set, nodes, dist_factor)
2173 # Set up sidesets
2174 if self.num_side_sets > 0:
2175 exo.put_side_set_names(self.get_side_set_names())
2176 exo.put_side_set_ids(self.get_side_set_ids())
2178 for side_set in self.get_side_set_ids():
2179 elements,sides = self.get_side_set_faces(side_set)
2180 try:
2181 dist_factor = self.get_side_set_dist_factors(side_set)
2182 except sdpy.fem.sdynpy_exodus.ExodusError:
2183 dist_factor = None
2184 exo.put_side_set_info(side_set, elements, sides, dist_factor)
2186 # Set up global variables
2187 if variables is None:
2188 names = self.get_global_variable_names()
2189 else:
2190 names = [name for name in self.get_global_variable_names() if name in variables]
2191 if len(names) > 0:
2192 exo.put_global_variable_names(names)
2194 if close:
2195 exo.close()
2197 return exo
2199 def repack(self, filename, repack_data, variables_to_repack = None, indices_in_repack = None, clobber = False, verbose = False, close = False):
2200 if verbose:
2201 print(f'Copying {self.filename} structure to {filename}')
2202 exo = self.copy_structure(filename, variables_to_repack, clobber)
2204 if isinstance(repack_data, TimeHistoryArray):
2205 if not repack_data.validate_common_abscissa():
2206 print('Warning: repack_data does not have common abscissa. The first abscissa entry will be used.')
2207 times = repack_data.ravel()[0].abscissa
2208 repack_data = repack_data.ordinate
2209 elif isinstance(repack_data, ShapeArray):
2210 times = repack_data.frequency
2211 repack_data = repack_data.shape_matrix.T
2212 else:
2213 repack_data = np.asarray(repack_data)
2214 times = np.arange(repack_data.shape[-1])
2216 # Set the times
2217 exo.set_times(times)
2218 num_times = len(times)
2220 # Go through node variables
2221 if variables_to_repack is None:
2222 variable_names = self.get_node_variable_names()
2223 else:
2224 variable_names = [name for name in self.get_node_variable_names() if name in variables_to_repack]
2225 for variable_name in variable_names:
2226 if verbose:
2227 print(f'Repacking Nodal Variable {variable_name}')
2228 data_matrix = self.get_node_variable_values(variable_name).T
2229 if indices_in_repack is not None:
2230 data_matrix = data_matrix[...,indices_in_repack]
2231 for step in range(num_times):
2232 q = repack_data[:,step]
2233 exo.set_node_variable_values(variable_name, step, data_matrix@q)
2235 # Go through element variables
2236 if variables_to_repack is None:
2237 variable_names = self.get_elem_variable_names()
2238 else:
2239 variable_names = [name for name in self.get_elem_variable_names() if name in variables_to_repack]
2240 for variable_name in variable_names:
2241 for block_id in self.get_elem_blk_ids():
2242 if verbose:
2243 print(f'Repacking Element Variable {variable_name} for block {block_id}')
2244 data_matrix = self.get_elem_variable_values(block_id, variable_name).T
2245 if indices_in_repack is not None:
2246 data_matrix = data_matrix[...,indices_in_repack]
2247 for step in range(num_times):
2248 q = repack_data[:,step]
2249 exo.set_elem_variable_values(block_id, variable_name, step, data_matrix@q)
2251 # Go through global variables
2252 if variables_to_repack is None:
2253 variable_names = self.get_global_variable_names()
2254 else:
2255 variable_names = [name for name in self.get_global_variable_names() if name in variables_to_repack]
2256 for variable_name in variable_names:
2257 if verbose:
2258 print(f'Repacking Global Variable {variable_name}')
2259 data_matrix = self.get_global_variable_values(variable_name)
2260 if indices_in_repack is not None:
2261 data_matrix = data_matrix[...,indices_in_repack]
2262 for step in range(num_times):
2263 q = repack_data[:,step]
2264 exo.set_global_variable_values(variable_name, step, data_matrix@q)
2266 if close:
2267 exo.close()
2269 return exo
2271 def __repr__(self):
2272 return_string = 'Exodus File at {:}'.format(self.filename)
2273 try:
2274 return_string += '\n {:} Timesteps'.format(self.num_times)
2275 except Exception:
2276 pass
2277 try:
2278 return_string += '\n {:} Nodes'.format(self.num_nodes)
2279 except Exception:
2280 pass
2281 try:
2282 return_string += '\n {:} Elements'.format(self.num_elems)
2283 except Exception:
2284 pass
2285 try:
2286 return_string += '\n Blocks: {:}'.format(', '.join(str(v)
2287 for v in self.get_elem_blk_ids()))
2288 except Exception:
2289 pass
2290 try:
2291 return_string += '\n Node Variables: {:}'.format(
2292 ', '.join(self.get_node_variable_names()))
2293 except Exception:
2294 pass
2295 try:
2296 return_string += '\n Element Variables: {:}'.format(
2297 ', '.join(self.get_elem_variable_names()))
2298 except Exception:
2299 pass
2300 try:
2301 return_string += '\n Global Variables: {:}'.format(
2302 ', '.join(self.get_global_variable_names()))
2303 except Exception:
2304 pass
2305 return return_string
2307 # def plot_mesh(self,surface_kwargs={'color':(0,0,1)},
2308 # bar_kwargs = {'color':(0,1,0),'tube_radius':None},
2309 # point_kwargs = {'color':(1,0,0),'scale_factor':0.1},
2310 # plot_surfaces = True, plot_bars = True, plot_points = True,
2311 # plot_edges = False):
2312 # '''Skins, triangulates, and plots a 3D representation of the mesh.
2314 # Parameters
2315 # ----------
2316 # mesh_kwargs : dict
2317 # A dictionary of keyword arguments for the rendering of surface
2318 # patches.
2319 # bar_kwargs : dict
2320 # A dictionary of keyword arguments for the rendering of lines (bar,
2321 # beam elements)
2322 # point_kwargs : dict
2323 # A dictionary of keyword arguments for the rendering of points
2324 # (sphere elements)
2325 # show : bool
2326 # A flag to specify whether or not to show the window.
2328 # Returns
2329 # -------
2330 # window : GLViewWidget
2331 # A reference to the view that the 3D content is plotted in.
2332 # meshes : list
2333 # A list of the mesh objects that have been added to the view.
2335 # '''
2336 # if mlab is None:
2337 # raise ModuleNotFoundError('Mayavi not installed!')
2338 # triangulated_mesh_info = self.triangulate_surface_mesh()
2339 # meshes = []
2340 # for block,node_indices,faces,original_elements in triangulated_mesh_info:
2341 # if faces.shape[1] == 3:
2342 # if plot_surfaces:
2343 # meshes.append(mlab.triangular_mesh(
2344 # *(self.get_coords()[:,node_indices]),faces,
2345 # **surface_kwargs))
2346 # if plot_edges:
2347 # meshes.append(mlab.triangular_mesh(
2348 # *(self.get_coords()[:,node_indices]),faces,
2349 # representation='wireframe',color=(0,0,0)))
2350 # elif faces.shape[1] == 2:
2351 # if plot_bars:
2352 # line_items = []
2353 # xs,ys,zs = self.get_coords()[:,node_indices][:,faces]
2354 # for x,y,z in zip(xs,ys,zs):
2355 # line_items.append(mlab.plot3d(x,y,z,**bar_kwargs))
2356 # meshes.append(line_items)
2357 # elif faces.shape[1] == 1:
2358 # if plot_points:
2359 # meshes.append(mlab.points3d(*(self.get_coords()[:,node_indices]),
2360 # **point_kwargs))
2361 # return meshes
2364class subfield(SimpleNamespace):
2365 def __init__(self, *args, **kwargs):
2366 super().__init__(*args, **kwargs)
2368 def __repr__(self):
2369 return_string = 'Field with subfields: {:}'.format(', '.join(['{:}'.format(key) if type(
2370 val) is np.ndarray else '{:}={:}'.format(key, val) for key, val in self.__dict__.items()]))
2371 return return_string
2373 def __str__(self):
2374 return repr(self)
2377class ExodusInMemory:
2378 '''Read or write exodus files loaded into memory
2380 This is a convenience class wrapped around the exodus class to enable
2381 easier manipulation of exodus files that fit entirely into memory
2383 Parameters
2384 ----------
2385 exo : exodus
2386 An exodus object, if you want to load an exodus file into memory, or
2387 None if you want to create an empty exodus file.
2388 variables : iterable
2389 A list of variable names that are loaded into memory. Default is
2390 to load all variables
2391 timesteps : iterable
2392 A list of timestep indices that are loaded into memory. Default is
2393 to load all timesteps.
2394 blocks : iterable
2395 A list of block ids that are loaded into memory. Default is to load
2396 all blocks.
2397 '''
2399 def __init__(self, exo=None, variables=None, timesteps=None, blocks=None):
2400 if exo is not None:
2401 self.load_from_exodus(exo, variables=variables, timesteps=timesteps, blocks=blocks)
2402 else:
2403 # Title
2404 self.title = ''
2405 # Attributes
2406 attribute_dict = {}
2407 attribute_dict['api_version'] = 5.22
2408 attribute_dict['version'] = 5.22
2409 attribute_dict['floating_point_word_size'] = 8
2410 attribute_dict['file_size'] = 1
2411 attribute_dict['int64_status'] = 0
2412 self.attributes = subfield(**attribute_dict)
2413 # QA Records
2414 self.qa_records = []
2415 # Info records
2416 self.info_records = [""]
2417 # Nodes
2418 node_dict = {}
2419 node_dict['coordinates'] = np.empty(0, dtype=float)
2420 node_dict['node_num_map'] = None
2421 node_dict['names'] = ('x', 'y', 'z')
2422 self.nodes = subfield(**node_dict)
2423 # TODO Coordinate Frames
2424 # Blocks
2425 self.blocks = []
2426 # Element Map
2427 self.element_map = None
2428 # Nodesets
2429 self.nodesets = None
2430 # Sidesets
2431 self.sidesets = None
2432 # Global Variables
2433 self.global_vars = None
2434 # Element Variables
2435 self.elem_vars = None
2436 # Nodal Variables
2437 self.nodal_vars = None
2438 # Time Steps
2439 self.time = []
2441 def load_from_exodus(self, exo, variables=None, timesteps=None, blocks=None):
2442 # Title
2443 self.title = exo.title
2444 # Attributes
2445 attribute_dict = {}
2446 attribute_dict['api_version'] = exo._ncdf_handle.api_version
2447 attribute_dict['version'] = exo._ncdf_handle.version
2448 attribute_dict['floating_point_word_size'] = exo._ncdf_handle.floating_point_word_size
2449 attribute_dict['file_size'] = exo._ncdf_handle.file_size
2450 try:
2451 attribute_dict['int64_status'] = exo._ncdf_handle.int64_status
2452 except AttributeError:
2453 attribute_dict['int64_status'] = 0
2454 self.attributes = subfield(**attribute_dict)
2456 # QA Records
2457 try:
2458 self.qa_records = exo.get_qa_records()
2459 except ExodusError:
2460 self.qa_records = []
2461 # Info Records
2462 try:
2463 self.info_records = exo.get_info_records()
2464 except ExodusError:
2465 self.info_records = [""]
2466 # Get Nodes
2467 node_dict = {}
2468 node_dict['coordinates'] = exo.get_coords().T
2469 node_dict['node_num_map'] = exo.get_node_num_map().data
2470 node_dict['names'] = exo.get_coord_names()
2471 self.nodes = subfield(**node_dict)
2472 # TODO Read Coordinate Frames
2473 # Blocks
2474 if blocks is None:
2475 block_ids = exo.get_elem_blk_ids()
2476 self.element_map = None
2477 else:
2478 block_ids = [block_id for block_id in exo.get_elem_blk_ids() if block_id in blocks]
2479 if timesteps is None:
2480 time_indices = np.arange(exo.get_times().size)
2481 else:
2482 time_indices = list(timesteps)
2483 self.blocks = []
2484 for index, blk_id in enumerate(block_ids):
2485 blk_dict = {}
2486 blk_dict['id'] = blk_id
2487 if not exo.num_elems == 0:
2488 blk_dict['connectivity'] = exo.get_elem_connectivity(blk_id).data
2489 blk_dict['elem_type'] = exo.get_elem_type(blk_id)
2490 blk_dict['name'] = ''.join(s.decode() for s in exo._ncdf_handle.variables['eb_names']
2491 [index, :] if not isinstance(s, np.ma.core.MaskedConstant))
2492 blk_dict['status'] = exo._ncdf_handle.variables['eb_status'][:].data[index]
2493 try:
2494 blk_dict['attributes'] = exo.get_elem_attr(blk_id)
2495 except ExodusError:
2496 blk_dict['attributes'] = None
2497 # TODO: Add attribute names once implemented
2498 try:
2499 blk_dict['attributes_name'] = None # exo.get_elem_attr_names(blk_id)
2500 except ExodusError:
2501 blk_dict['attributes_name'] = None
2502 self.blocks.append(subfield(**blk_dict))
2503 # Element Map
2504 self.element_map = exo.get_elem_num_map().data
2506 # Get Nodesets
2507 try:
2508 nodeset_ids = exo.get_node_set_ids()
2509 self.nodesets = []
2510 for index, ns_id in enumerate(nodeset_ids):
2511 ns_dict = {}
2512 ns_dict['id'] = ns_id
2513 ns_dict['nodes'] = exo.get_node_set_nodes(ns_id)
2514 try:
2515 ns_dict['dist_factors'] = exo.get_node_set_dist_factors(ns_id)
2516 except ExodusError:
2517 ns_dict['dist_factors'] = None
2518 try:
2519 ns_dict['name'] = exo.get_node_set_names()[index]
2520 except ExodusError:
2521 ns_dict['name'] = None
2522 ns_dict['status'] = exo._ncdf_handle.variables['ns_status'][:].data[index]
2523 self.nodesets.append(subfield(**ns_dict))
2524 except ExodusError:
2525 self.nodesets = None
2527 # Get Sidesets
2528 try:
2529 sideset_ids = exo.get_side_set_ids()
2530 self.sidesets = []
2531 for index, ss_id in enumerate(sideset_ids):
2532 ss_dict = {}
2533 ss_dict['id'] = ss_id
2534 ss_dict['elements'], ss_dict['sides'] = exo.get_side_set_faces(ss_id)
2535 try:
2536 ss_dict['dist_factors'] = exo.get_side_set_dist_factors(ss_id)
2537 except ExodusError:
2538 ss_dict['dist_factors'] = None
2539 ss_dict['status'] = exo._ncdf_handle.variables['ss_status'][:].data[index]
2540 try:
2541 ss_dict['name'] = exo.get_side_set_names()[index]
2542 except ExodusError:
2543 ss_dict['name'] = None
2544 self.sidesets.append(subfield(**ss_dict))
2545 except ExodusError:
2546 self.sidesets = None
2548 # Get Global Variables
2549 try:
2550 num_glob_vars = exo.num_global_variables
2551 if num_glob_vars == 0:
2552 self.global_vars = None
2553 else:
2554 self.global_vars = []
2555 global_var_names = exo.get_global_variable_names()
2556 if variables is not None:
2557 global_var_names = [name for name in global_var_names if name in variables]
2558 if len(global_var_names) == 0:
2559 raise ExodusError('No Global Variables to Load')
2560 for index, name in enumerate(global_var_names):
2561 global_var_dict = {}
2562 global_var_dict['name'] = name
2563 global_var_dict['data'] = exo.get_global_variable_values(name)[time_indices]
2564 self.global_vars.append(subfield(**global_var_dict))
2565 except ExodusError:
2566 self.global_vars = None
2568 # Get Nodal Variables
2569 try:
2570 num_node_vars = exo.num_node_variables
2571 if num_node_vars == 0:
2572 self.nodal_vars = None
2573 else:
2574 self.nodal_vars = []
2575 node_var_names = exo.get_node_variable_names()
2576 if variables is not None:
2577 node_var_names = [name for name in node_var_names if name in variables]
2578 if len(node_var_names) == 0:
2579 raise ExodusError('No Nodal Variables to Load')
2580 for index, name in enumerate(node_var_names):
2581 node_var_dict = {}
2582 node_var_dict['name'] = name
2583 node_var_dict['data'] = exo.get_node_variable_values(name).data[time_indices, :]
2584 self.nodal_vars.append(subfield(**node_var_dict))
2585 except ExodusError:
2586 self.nodal_vars = None
2588 # Get Element Variables
2589 try:
2590 num_elem_vars = exo.num_elem_variables
2591 if num_elem_vars == 0:
2592 self.elem_vars = None
2593 else:
2594 self.elem_vars = []
2595 elem_var_names = exo.get_elem_variable_names()
2596 variable_truth_table = exo.get_elem_variable_table().data
2597 for blk_index, blk_id in enumerate(exo.get_elem_blk_ids()):
2598 if blk_id not in block_ids:
2599 continue
2600 block_dict = {}
2601 block_dict['id'] = blk_id
2602 block_dict['elem_var_data'] = []
2603 for elem_index, elem_name in enumerate(elem_var_names):
2604 if (variables is not None) and (elem_name not in variables):
2605 continue
2606 if variable_truth_table[blk_index, elem_index] == 0:
2607 block_dict['elem_var_data'].append(None)
2608 else:
2609 elem_var_dict = {}
2610 elem_var_dict['name'] = elem_name
2611 elem_var_dict['data'] = exo.get_elem_variable_values(
2612 blk_id, elem_index).data[time_indices, :]
2613 block_dict['elem_var_data'].append(subfield(**elem_var_dict))
2614 if len(block_dict['elem_var_data']) == 0:
2615 raise ExodusError('No Element Variables to Load')
2616 self.elem_vars.append(subfield(**block_dict))
2617 except ExodusError:
2618 self.elem_vars = None
2620 # Get Time Steps
2621 self.time = exo.get_times(time_indices).data
2623 @staticmethod
2624 def from_sdynpy(geometry, displacement_data=None):
2625 fexo = ExodusInMemory()
2627 # Put in the node positions
2628 node_data = np.sort(geometry.node)
2629 definition_coordinate_systems = geometry.coordinate_system(node_data.def_cs)
2630 global_positions = global_coord(definition_coordinate_systems, node_data.coordinate)
2631 fexo.nodes.coordinates = global_positions
2632 fexo.nodes.node_num_map = node_data.id
2633 node_map_dict = {node_id: node_index for node_index, node_id in enumerate(node_data.id)}
2634 node_map_dict[0] = -1 # For tracelines, map 0 (pick up pen) to -1 (something easy to throw out)
2635 node_index_map = np.vectorize(
2636 node_map_dict.__getitem__)
2637 # Now go through and add in the elements
2638 element_data = geometry.element
2639 element_types = np.unique(element_data.type)
2640 elem_colors = []
2641 for elem_type in element_types:
2642 blk_dict = {}
2643 blk_dict['id'] = elem_type
2644 blk_dict['name'] = 'block_{:}'.format(elem_type)
2645 elements_this_type = element_data[element_data.type == elem_type]
2646 connectivity = np.stack(elements_this_type.connectivity)
2647 blk_dict['connectivity'] = node_index_map(connectivity)
2648 blk_dict['elem_type'] = _inverse_exodus_elem_type_map[elem_type]
2649 blk_dict['status'] = 1
2650 blk_dict['attributes'] = None
2651 blk_dict['attributes_name'] = None
2652 elem_colors.append(elements_this_type.color)
2653 fexo.blocks.append(subfield(**blk_dict))
2655 # Now do the same for the tracelines
2656 traceline_data = geometry.traceline
2657 for i, (key, traceline) in enumerate(traceline_data.ndenumerate()):
2658 connectivity_1d = node_index_map(np.stack(traceline.connectivity))
2659 connectivity = np.concatenate((connectivity_1d[:-1, np.newaxis],
2660 connectivity_1d[1:, np.newaxis]), axis=-1)
2661 connectivity = connectivity[~np.any(connectivity == -1, axis=1)] # We mapped 0s to -1, so we want to throw those out
2662 blk_dict = {}
2663 blk_dict['id'] = i + 500
2664 blk_dict['name'] = traceline.description
2665 blk_dict['connectivity'] = connectivity
2666 blk_dict['elem_type'] = 'bar'
2667 blk_dict['status'] = 1
2668 blk_dict['attributes'] = None
2669 blk_dict['attributes_name'] = None
2670 elem_colors.append(np.ones(connectivity.shape[0]) * traceline.color)
2671 fexo.blocks.append(subfield(**blk_dict))
2673 from ..core.sdynpy_shape import ShapeArray
2674 from ..core.sdynpy_data import NDDataArray, TimeHistoryArray
2675 if isinstance(displacement_data, ShapeArray):
2676 shapes = displacement_data.flatten()
2677 fexo.time = shapes.frequency
2678 # Add global variables frequency and damping
2679 for name in ['frequency', 'damping']:
2680 global_var_dict = {}
2681 global_var_dict['name'] = name
2682 global_var_dict['data'] = shapes[name]
2683 # Now add displacement data for the shapes
2684 coordinates = shapes[0].coordinate.copy()
2685 indices = np.isin(coordinates.node, node_data.id)
2686 coordinates = coordinates[indices]
2687 # Keep only displacements
2688 coordinates = coordinates[(abs(coordinates.direction) <= 3) &
2689 (abs(coordinates.direction) >= 1)]
2690 local_deformation_directions = coordinates.local_direction()
2691 coordinate_nodes = node_data(coordinates.node)
2692 displacement_coordinate_systems = geometry.coordinate_system(coordinate_nodes.disp_cs)
2693 shape_displacements = shapes[coordinates]
2694 coordinates.node = node_index_map(coordinates.node)
2695 points = global_positions[coordinates.node]
2696 global_deformation_directions = global_deflection(displacement_coordinate_systems,
2697 local_deformation_directions,
2698 points)
2699 node_displacements = np.zeros((shape_displacements.shape[0],)
2700 + node_data.shape
2701 + (3,), dtype=shape_displacements.dtype)
2703 for coordinate, direction, shape_coefficients in zip(coordinates, global_deformation_directions, shape_displacements.T):
2704 node_index = coordinate.node
2705 displacements = shape_coefficients[:, np.newaxis] * direction
2706 node_displacements[:, node_index, :] += displacements
2708 # Now add the variables
2709 fexo.nodal_vars = []
2710 for i, name in enumerate(('DispX', 'DispY', 'DispZ')):
2711 node_var_dict = {}
2712 node_var_dict['name'] = name
2713 node_var_dict['data'] = node_displacements[..., i]
2714 fexo.nodal_vars.append(subfield(**node_var_dict))
2715 elif isinstance(displacement_data, NDDataArray):
2716 if isinstance(displacement_data, TimeHistoryArray):
2717 data = displacement_data.flatten()
2718 fexo.time = data[0].abscissa
2719 # Now add displacement data for the shapes
2720 coordinates = data.response_coordinate
2721 indices = np.isin(coordinates.node, node_data.id)
2722 coordinates = coordinates[indices]
2723 # Keep only displacements
2724 coordinates = coordinates[(abs(coordinates.direction) <= 3)
2725 & (abs(coordinates.direction) >= 1)]
2726 local_deformation_directions = coordinates.local_direction()
2727 coordinate_nodes = node_data(coordinates.node)
2728 displacement_coordinate_systems = geometry.coordinate_system(
2729 coordinate_nodes.disp_cs)
2730 data_displacements = data[coordinates[:, np.newaxis]].ordinate.T
2731 coordinates.node = node_index_map(coordinates.node)
2732 points = global_positions[coordinates.node]
2733 global_deformation_directions = global_deflection(displacement_coordinate_systems,
2734 local_deformation_directions,
2735 points)
2736 node_displacements = np.zeros((data_displacements.shape[0],)
2737 + node_data.shape
2738 + (3,), dtype=data_displacements.dtype)
2739 for coordinate, direction, shape_coefficients in zip(coordinates, global_deformation_directions, data_displacements.T):
2740 node_index = coordinate.node
2741 displacements = shape_coefficients[:, np.newaxis] * direction
2742 node_displacements[:, node_index, :] += displacements
2744 # Now add the variables
2745 fexo.nodal_vars = []
2746 for i, name in enumerate(('DispX', 'DispY', 'DispZ')):
2747 node_var_dict = {}
2748 node_var_dict['name'] = name
2749 node_var_dict['data'] = node_displacements[..., i]
2750 fexo.nodal_vars.append(subfield(**node_var_dict))
2751 else:
2752 raise NotImplementedError(
2753 'Saving Functions with type {:} to Exodus Data is Not Implemented Yet'.format(type(displacement_data)))
2754 else:
2755 fexo.time = np.array([0])
2756 fexo.nodal_vars = []
2758 # Now add color information
2759 node_var_dict = {}
2760 node_var_dict['name'] = 'NodeColor'
2761 node_var_dict['data'] = np.tile(node_data.color[np.newaxis, :], (len(fexo.time), 1))
2762 fexo.nodal_vars.append(subfield(**node_var_dict))
2764 fexo.elem_vars = []
2765 for block, color in zip(fexo.blocks, elem_colors):
2766 block_dict = {}
2767 block_dict['id'] = block.id
2768 elem_var_dict = {}
2769 elem_var_dict['name'] = 'ElemColor'
2770 elem_var_dict['data'] = color * np.ones((len(fexo.time), 1))
2772 block_dict['elem_var_data'] = [subfield(**elem_var_dict)]
2773 fexo.elem_vars.append(subfield(**block_dict))
2775 return fexo
2777 def write_to_file(self, filename, clobber=False):
2778 num_nodes, num_dims = self.nodes.coordinates.shape
2779 num_elem = sum([block.connectivity.shape[0] for block in self.blocks])
2780 num_blocks = len(self.blocks)
2781 if self.nodesets is None:
2782 num_node_sets = 0
2783 else:
2784 num_node_sets = len(self.nodesets)
2785 if self.sidesets is None:
2786 num_side_sets = 0
2787 else:
2788 num_side_sets = len(self.sidesets)
2789 exo_out = Exodus(filename, 'w', self.title,
2790 num_dims, num_nodes, num_elem, num_blocks, num_node_sets,
2791 num_side_sets, clobber=clobber)
2793 # Write time steps
2794 exo_out.set_times(self.time)
2796 # Initialize Information
2797 exo_out.put_info_records(self.info_records)
2798 qa_records = list(self.qa_records)
2799 qa_records.append(('pyexodus', __version__, str(
2800 datetime.datetime.now().date()), str(datetime.datetime.now().time())))
2801 exo_out.put_qa_records(qa_records)
2803 # Write the nodes
2804 exo_out.put_coord_names(self.nodes.names)
2805 exo_out.put_coords(self.nodes.coordinates.T)
2806 if self.nodes.node_num_map is not None:
2807 exo_out.put_node_num_map(self.nodes.node_num_map)
2809 # Initialize blocks
2810 exo_out.put_elem_blk_ids([block.id for block in self.blocks])
2811 for block in self.blocks:
2812 if block.attributes is None:
2813 nattr_per_elem = 0
2814 else:
2815 nattr_per_elem = block.attributes.shape[1]
2816 exo_out.put_elem_blk_info(block.id,
2817 block.elem_type,
2818 block.connectivity.shape[0],
2819 block.connectivity.shape[1],
2820 nattr_per_elem)
2821 exo_out.set_elem_connectivity(block.id, block.connectivity)
2822 if nattr_per_elem > 0:
2823 exo_out.set_elem_attr(block.id, block.attributes)
2825 # Element Map
2826 if self.element_map is not None:
2827 exo_out.put_elem_num_map(self.element_map)
2829 # Nodesets
2830 if self.nodesets is not None:
2831 exo_out.put_node_set_ids([ns.id for ns in self.nodesets])
2832 exo_out.put_node_set_names(
2833 ['' if nodeset.name is None else nodeset.name for nodeset in self.nodesets])
2834 for nodeset in self.nodesets:
2835 exo_out.put_node_set_info(nodeset.id, nodeset.nodes, nodeset.dist_factors)
2837 # Sidesets
2838 if self.sidesets is not None:
2839 exo_out.put_side_set_ids([ns.id for ns in self.sidesets])
2840 exo_out.put_side_set_names(
2841 ['' if sideset.name is None else sideset.name for sideset in self.sidesets])
2842 for sideset in self.sidesets:
2843 exo_out.put_side_set_info(sideset.id, sideset.elements,
2844 sideset.sides, sideset.dist_factors)
2846 # Global Variables #TODO: Make it so I can write all timesteps at once.
2847 if self.global_vars is not None:
2848 exo_out.put_global_variable_names([var.name for var in self.global_vars])
2849 for i, var in enumerate(self.global_vars):
2850 for j, time in enumerate(self.time):
2851 exo_out.set_global_variable_values(i, j, self.global_vars[i].data[j])
2853 # Nodal Variables #TODO: Make it so I can write all timesteps at once.
2854 if self.nodal_vars is not None:
2855 exo_out.put_node_variable_names([var.name for var in self.nodal_vars])
2856 for i, var in enumerate(self.nodal_vars):
2857 for j, time in enumerate(self.time):
2858 exo_out.set_node_variable_values(i, j, self.nodal_vars[i].data[j, :])
2860 # Element Variables # TODO: Make it so I can write all timesteps at once
2861 if self.elem_vars is not None:
2862 nblocks = len(self.elem_vars)
2863 nelemvars = len(self.elem_vars[0].elem_var_data)
2864 elem_table = np.zeros((nblocks, nelemvars), dtype=int)
2865 for i in range(nblocks):
2866 for j in range(nelemvars):
2867 elem_table[i, j] = not self.elem_vars[i].elem_var_data[j] is None
2868 # Detect if any variables are not used ever
2869 variable_names = []
2870 keep_booleans = []
2871 for j in range(nelemvars):
2872 indices = np.nonzero(elem_table[:, j])[0]
2873 if len(indices) == 0:
2874 variable_names.append(None)
2875 keep_booleans.append(False)
2876 else:
2877 variable_names.append(self.elem_vars[indices[0]].elem_var_data[j].name)
2878 keep_booleans.append(True)
2879 elem_table = elem_table[:, keep_booleans]
2880 variable_names = [name for name, boolean in zip(
2881 variable_names, keep_booleans) if boolean]
2882 keep_indices = np.nonzero(keep_booleans)[0]
2883 exo_out.put_elem_variable_names(variable_names, elem_table)
2884 for i, block in enumerate(self.blocks):
2885 for j, name in enumerate(variable_names):
2886 if elem_table[i, j]:
2887 for k, time in enumerate(self.time):
2888 exo_out.set_elem_variable_values(block.id, j, k,
2889 self.elem_vars[i].elem_var_data[keep_indices[j]].data[k, :])
2891 exo_out.close()
2893 def repack(self, q, modes=None):
2894 '''Repackages an exodus file as a linear combination of itself'''
2895 if modes is None:
2896 mode_indices = np.arange(len(self.time))
2897 else:
2898 mode_indices = modes
2900 if q.shape[0] != len(mode_indices):
2901 raise ValueError(
2902 'Number of rows in q must be equivalent to number of time steps in exodus file')
2904 out_fexo = copy.deepcopy(self)
2906 out_fexo.time = np.zeros((q.shape[1]), dtype='float')
2908 if self.global_vars is not None:
2909 for i, global_var in enumerate(self.global_vars):
2910 out_fexo.global_vars[i].data = (
2911 global_var.data[np.newaxis, mode_indices] @ q).squeeze()
2913 if self.nodal_vars is not None:
2914 for i, nodal_var in enumerate(self.nodal_vars):
2915 out_fexo.nodal_vars[i].data = (nodal_var.data.T[:, mode_indices] @ q).T
2917 if self.elem_vars is not None:
2918 for i, block in enumerate(self.elem_vars):
2919 for j, elem_vars in enumerate(block.elem_var_data):
2920 if elem_vars is not None:
2921 out_fexo.elem_vars[i].elem_var_data[j].data = (
2922 elem_vars.data.T[:, mode_indices] @ q).T
2924 return out_fexo
2926 def get_block_surface(self, block_id, keep_midside_nodes=False, warn=True):
2927 '''Gets the node indices and element connectivity of surface elements
2929 This function "skins" the element block, returning a list of node
2930 indices and a surface connectivity matrix.
2932 Parameters
2933 ----------
2934 block_id : int
2935 The ID number of the block that will be skinned.
2936 keep_midside_nodes : bool
2937 Specifies whether or not to keep midside nodes in the surface mesh.
2938 Default is False.
2939 warn : bool
2940 Specifies whether or not to warn the user if the block ID doesn't
2941 have a skinning method defined for its element type. Default is
2942 True.
2944 Returns
2945 -------
2946 element_block_information : list
2947 A list of tuples of element information. These data are
2948 element_type, node_indices, block_surface_connectivity, and
2949 block_surface_original_elements. The element_type is a string
2950 representing the new block element type ('quad4','tri3',etc.). The
2951 node_indices can be used as an index into the coordinate or nodal
2952 variable arrays to select nodes corresponding to this block. The
2953 block_surface_connectivity represents the connectivity array of the
2954 surface faces of the block. Values in this array correspond to
2955 indices into the node_indices array. To recover the connectivity
2956 array in the original node indices of the exodus file, it can be
2957 passed through the node_indices as
2958 node_indices[block_surface_connectivity]. The
2959 block_surface_original_elements array shows the original element
2960 indices of the block that each surface came from. This can be
2961 used to map element variables to the new surface mesh. This list
2962 will normally be length 1 unless an element type is processed that
2963 has two different surface elements in it (e.g. wedges have
2964 tris and quads)
2965 '''
2966 block_index = [block.id for block in self.blocks].index(block_id)
2967 elem_type = self.blocks[block_index].elem_type
2968 connectivity = self.blocks[block_index].connectivity
2969 face_connectivity_array = face_connectivity(elem_type, keep_midside_nodes)
2970 element_block_information = []
2971 for reduced_elem_type, reduced_connectivity in face_connectivity_array:
2972 if reduced_elem_type is None:
2973 if warn:
2974 print('Warning: Element Type {:} has no defined face reduction. Passing connectivity unchanged.'.format(
2975 elem_type))
2976 block_face_connectivity = connectivity
2977 block_elem_type = elem_type
2978 block_face_original_elements = np.arange(connectivity.shape[0])
2979 else:
2980 block_face_connectivity = connectivity[:,
2981 reduced_connectivity].reshape(-1, reduced_connectivity.shape[-1])
2982 block_face_original_elements = np.repeat(np.arange(connectivity.shape[0]),
2983 reduced_connectivity.shape[0])
2984 block_elem_type = reduced_elem_type
2985 # Remove faces that are duplicates
2986 (unique_rows, unique_row_indices, unique_inverse, unique_counts) = np.unique(np.sort(
2987 block_face_connectivity, axis=1), return_index=True, return_inverse=True, return_counts=True, axis=0)
2988 original_unique_counts = unique_counts[unique_inverse]
2989 nondupe_faces = original_unique_counts == 1
2990 block_face_connectivity = block_face_connectivity[nondupe_faces, :]
2991 block_face_original_elements = block_face_original_elements[nondupe_faces]
2993 node_indices = np.unique(block_face_connectivity)
2994 node_map = np.zeros(self.nodes.coordinates.shape[0], dtype=int)
2995 node_map[node_indices] = np.arange(len(node_indices))
2996 block_face_connectivity = node_map[block_face_connectivity]
2997 element_block_information.append((block_elem_type,
2998 node_indices,
2999 block_face_connectivity,
3000 block_face_original_elements))
3001 return element_block_information
3003 def triangulate_surface_mesh(self):
3004 '''Triangulate a surface mesh for plotting patches
3006 This function generates a triangle mesh for each block in the model if
3007 it can. If there are more than 3 nodes per element in a block, and the
3008 triangulation scheme hasn't been defined in
3009 pyexodus.mesh_triangulation_array, it will be skipped.
3011 Parameters
3012 ----------
3013 None
3015 Returns
3016 -------
3017 triangulated_mesh_info : list
3018 A list of tuples containing block id, node_indices, triangulated
3019 connectivity, and original block elements
3020 '''
3021 triangulated_mesh_info = []
3022 for block in self.blocks:
3023 surface_mesh_info = self.get_block_surface(block.id, warn=False)
3024 for elem_type, node_indices, connectivity, original_elements in surface_mesh_info:
3025 triangulation_scheme = mesh_triangulation_array(elem_type)
3026 if triangulation_scheme is None:
3027 if connectivity.shape[-1] > 3:
3028 print('Warning: More than 3 ({:}) nodes per element in block {:}, this block will not be triangulated'.format(
3029 connectivity.shape[-1], block.id))
3030 else:
3031 triangulated_mesh_info.append(
3032 (block.id, node_indices, connectivity, original_elements))
3033 else:
3034 triangulated_mesh_info.append(
3035 (block.id,
3036 node_indices,
3037 connectivity[:,
3038 triangulation_scheme].reshape(-1, triangulation_scheme.shape[-1]),
3039 np.repeat(original_elements, triangulation_scheme.shape[0])))
3040 return triangulated_mesh_info
3042 # def plot_mesh(self,surface_kwargs={'color':(0,0,1)},
3043 # bar_kwargs = {'color':(0,1,0),'tube_radius':None},
3044 # point_kwargs = {'color':(1,0,0),'scale_factor':0.1},
3045 # plot_surfaces = True, plot_bars = True, plot_points = True,
3046 # plot_edges = False):
3047 # '''Skins, triangulates, and plots a 3D representation of the mesh.
3049 # Parameters
3050 # ----------
3051 # mesh_kwargs : dict
3052 # A dictionary of keyword arguments for the rendering of surface
3053 # patches.
3054 # bar_kwargs : dict
3055 # A dictionary of keyword arguments for the rendering of lines (bar,
3056 # beam elements)
3057 # point_kwargs : dict
3058 # A dictionary of keyword arguments for the rendering of points
3059 # (sphere elements)
3060 # show : bool
3061 # A flag to specify whether or not to show the window.
3063 # Returns
3064 # -------
3065 # window : GLViewWidget
3066 # A reference to the view that the 3D content is plotted in.
3067 # meshes : list
3068 # A list of the mesh objects that have been added to the view.
3070 # '''
3071 # if mlab is None:
3072 # raise ModuleNotFoundError('Mayavi not installed!')
3073 # triangulated_mesh_info = self.triangulate_surface_mesh()
3074 # meshes = []
3075 # for block,node_indices,faces,original_elements in triangulated_mesh_info:
3076 # if faces.shape[1] == 3:
3077 # if plot_surfaces:
3078 # meshes.append(mlab.triangular_mesh(
3079 # *(self.nodes.coordinates.T[:,node_indices]),faces,
3080 # **surface_kwargs))
3081 # if plot_edges:
3082 # meshes.append(mlab.triangular_mesh(
3083 # *(self.nodes.coordinates.T[:,node_indices]),faces,
3084 # representation='wireframe',color=(0,0,0)))
3085 # elif faces.shape[1] == 2:
3086 # if plot_bars:
3087 # line_items = []
3088 # xs,ys,zs = self.nodes.coordinates.T[:,node_indices][:,faces]
3089 # for x,y,z in zip(xs,ys,zs):
3090 # line_items.append(mlab.plot3d(x,y,z,**bar_kwargs))
3091 # meshes.append(line_items)
3092 # elif faces.shape[1] == 1:
3093 # if plot_points:
3094 # meshes.append(mlab.points3d(*(self.nodes.coordinates.T[:,node_indices]),
3095 # **point_kwargs))
3096 # return meshes
3098 def reduce_to_surfaces(self,*args,**kwargs):
3099 return reduce_exodus_to_surfaces(self,*args,**kwargs)
3101 def extract_sharp_edges(self,*args,**kwargs):
3102 return extract_sharp_edges(self, *args,**kwargs)
3105def reduce_exodus_to_surfaces(fexo, blocks_to_transform=None, variables_to_transform=None, keep_midside_nodes=False, verbose=False):
3106 '''Convert exodus finite element models to surface elements
3108 This function converts specified volume meshes in an exodus file to surface
3109 meshes to ease computational complexity.
3111 Parameters
3112 ----------
3113 fexo : ExodusInMemory object or Exodus object
3114 fexo should be an in-memory representation of the finite element mesh
3115 blocks_to_transform : iterable
3116 blocks_to_transform includes all of the blocks that will be included
3117 in the output model
3118 variables_to_transform : iterable
3119 variables_to_transform is a case-insensitive list of all variable names
3120 that will be kept in the final model
3121 keep_midside_nodes : bool
3122 keep_midside_nodes specifies whether or not to transform quadratic
3123 elements into linear elements, discarding any non-corner nodes.
3125 Returns
3126 -------
3127 fexo_out : ExodusInMemory object
3128 An equivalent fexo reduced to the surface geometry.
3130 Notes
3131 -----
3132 TO ADD MORE ELEMENT TYPES, we need to define a volume element name that
3133 we will find in the exodus file, for example 'tetra4' or 'hex20', the
3134 number of nodes per face on the element, and a connectivity matrix that
3135 specifies how the faces are made from the nodes of the element. For
3136 example a hex8 element has the following structure::
3138 8--------------7
3139 /| /|
3140 / | / |
3141 / | / |
3142 5--------------6 |
3143 | 4----------|---3
3144 | / | /
3145 | / | /
3146 |/ |/
3147 1--------------2
3149 So the 6 faces are as follows::
3151 1,2,6,5
3152 2,3,7,6
3153 3,4,8,7
3154 4,1,5,8
3155 4,3,2,1
3156 5,6,7,8
3158 We create the element face connectivity by simply mashing these together
3159 one after another like so: 1,2,6,5,2,3,7,6,3,4,8,7,4,1,5,8,4,3,2,1,5,6,7,8
3160 We lastly need to specify what we are turning the element into, in this
3161 case, quad4.
3162 So we add an entry to the reduce_element_types,
3163 reduce_element_nodes_per_face, reduce_element_face, and
3164 reduce_element_substitute_type variables containing this information
3165 hex8::
3167 reduce_element_types{4} = 'hex8'
3168 reduce_element_nodes_per_face{4} = 4
3169 reduce_element_face{4} = [1,2,6,5,2,3,7,6,3,4,8,7,4,1,5,8,4,3,2,1,5,6,7,8]
3170 reduce_element_substitute_type{4} = 'quad4'
3172 To see what other elements look like, see page 18 of the exodusII manual:
3173 http://prod.sandia.gov/techlib/access-control.cgi/1992/922137.pdf
3174 '''
3176 if isinstance(fexo, Exodus):
3177 exo_in_mem = False
3178 else:
3179 exo_in_mem = True
3181 # If blocks_to_transform is not specified, keep all blocks
3182 if blocks_to_transform is None:
3183 if exo_in_mem:
3184 blocks_to_transform = [block.id for block in fexo.blocks]
3185 else:
3186 blocks_to_transform = fexo.get_elem_blk_ids()
3188 # if variables_to_keep is not specified, keep all variables
3189 if variables_to_transform is None:
3190 if exo_in_mem:
3191 try:
3192 node_variables_to_transform = np.arange(len(fexo.nodal_vars))
3193 except TypeError:
3194 node_variables_to_transform = []
3195 try:
3196 elem_variables_to_transform = np.arange(len(fexo.elem_vars[0].elem_var_data))
3197 except TypeError:
3198 elem_variables_to_transform = []
3199 else:
3200 try:
3201 node_variables_to_transform = np.arange(fexo.num_node_variables)
3202 except ExodusError:
3203 node_variables_to_transform = []
3204 try:
3205 elem_variables_to_transform = np.arange(fexo.num_elem_variables)
3206 except ExodusError:
3207 elem_variables_to_transform = []
3208 else:
3209 if exo_in_mem:
3210 try:
3211 node_variables_to_transform = np.array([i for i, nodal_variable in enumerate(
3212 fexo.nodal_vars) if nodal_variable.lower() in [var.lower() for var in variables_to_transform]])
3213 except TypeError:
3214 node_variables_to_transform = []
3215 # Need to do extra processing for element variable names because
3216 # not all variable names are in all blocks (might be none)
3217 try:
3218 elem_var_names = []
3219 for i in range(len(fexo.elem_vars[0].elem_var_data)):
3220 var_name = None
3221 for elem_var in fexo.elem_vars:
3222 try:
3223 var_name = elem_var.elem_var_data[i].name
3224 break
3225 except AttributeError:
3226 continue
3227 elem_var_names.append(var_name)
3228 elem_variables_to_transform = np.array([i for i, name in enumerate(
3229 elem_var_names) if name.lower() in [var.lower() for var in variables_to_transform]])
3230 except TypeError:
3231 elem_variables_to_transform = []
3232 else:
3233 try:
3234 node_variables_to_transform = np.array([i for i, name in enumerate(
3235 fexo.get_node_variable_names()) if name.lower() in [var.lower() for var in variables_to_transform]])
3236 except ExodusError:
3237 node_variables_to_transform = []
3238 try:
3239 elem_variables_to_transform = np.array([i for i, name in enumerate(
3240 fexo.get_elem_variable_names()) if name.lower() in [var.lower() for var in variables_to_transform]])
3241 except ExodusError:
3242 elem_variables_to_transform = []
3243 if verbose:
3244 print('Keeping {:} nodal variables and {:} element variables'.format(
3245 len(node_variables_to_transform), len(elem_variables_to_transform)))
3247 # Define element types and what we will do to reduce them
3248 reduced_element_types = []
3249 reduced_element_nodes_per_face = []
3250 reduced_element_face_connectivity = []
3251 reduced_element_substitute_type = []
3253 # 10-node tetrahedral
3254 reduced_element_types.append('tetra10')
3255 if keep_midside_nodes:
3256 reduced_element_nodes_per_face.append(6)
3257 reduced_element_face_connectivity.append([0, 1, 3, 4, 8, 7,
3258 1, 2, 3, 5, 9, 8,
3259 2, 0, 3, 6, 7, 9,
3260 0, 2, 1, 6, 5, 4])
3261 reduced_element_substitute_type.append('tri6')
3262 else:
3263 reduced_element_nodes_per_face.append(3)
3264 reduced_element_face_connectivity.append([0, 1, 3,
3265 1, 2, 3,
3266 2, 0, 3,
3267 0, 2, 1])
3268 reduced_element_substitute_type.append('tri3')
3269 # 4-node tetrahedral
3270 reduced_element_types.append('tetra4')
3271 reduced_element_nodes_per_face.append(3)
3272 reduced_element_face_connectivity.append([0, 1, 3,
3273 1, 2, 3,
3274 2, 0, 3,
3275 0, 2, 1])
3276 reduced_element_substitute_type.append('tri3')
3277 # 20-node hexahedral
3278 reduced_element_types.append('hex20')
3279 if keep_midside_nodes:
3280 reduced_element_nodes_per_face.append(8)
3281 reduced_element_face_connectivity.append([0, 1, 5, 4, 8, 13, 16, 12,
3282 1, 2, 6, 5, 9, 14, 17, 13,
3283 2, 3, 7, 6, 10, 15, 18, 14,
3284 3, 0, 4, 7, 11, 12, 19, 15,
3285 3, 2, 1, 0, 11, 8, 9, 10,
3286 4, 5, 6, 7, 16, 17, 18, 19])
3287 reduced_element_substitute_type.append('quad8')
3288 else:
3289 reduced_element_nodes_per_face.append(4)
3290 reduced_element_face_connectivity.append([0, 1, 5, 4,
3291 1, 2, 6, 5,
3292 2, 3, 7, 6,
3293 3, 0, 4, 7,
3294 3, 2, 1, 0,
3295 4, 5, 6, 7])
3296 reduced_element_substitute_type.append('quad4')
3297 # 8-node hexahedral
3298 reduced_element_types.append('hex8')
3299 reduced_element_nodes_per_face.append(4)
3300 reduced_element_face_connectivity.append([0, 1, 5, 4,
3301 1, 2, 6, 5,
3302 2, 3, 7, 6,
3303 3, 0, 4, 7,
3304 3, 2, 1, 0,
3305 4, 5, 6, 7])
3306 reduced_element_substitute_type.append('quad4')
3307 reduced_element_types.append('hex')
3308 reduced_element_nodes_per_face.append(4)
3309 reduced_element_face_connectivity.append([0, 1, 5, 4,
3310 1, 2, 6, 5,
3311 2, 3, 7, 6,
3312 3, 0, 4, 7,
3313 3, 2, 1, 0,
3314 4, 5, 6, 7])
3315 reduced_element_substitute_type.append('quad4')
3317 # Keep track of which nodes are used in the final model
3318 if exo_in_mem:
3319 keep_nodes = np.zeros(fexo.nodes.coordinates.shape[0], dtype=bool)
3320 else:
3321 keep_nodes = np.zeros(fexo.num_nodes, dtype=bool)
3323 # Create the output variable
3324 fexo_out = ExodusInMemory()
3325 # Set the same times
3326 if exo_in_mem:
3327 fexo_out.time = copy.deepcopy(fexo.time)
3328 else:
3329 fexo_out.time = fexo.get_times().data
3331 # Since we are just visualizing, eliminate exodus features that we don't
3332 # care about
3333 fexo_out.element_map = None
3334 fexo_out.sidesets = None
3335 fexo_out.nodesets = None
3337 if exo_in_mem:
3338 block_inds_to_transform = [(block_index, block.id) for block_index, block in enumerate(
3339 fexo.blocks) if block.id in blocks_to_transform]
3340 else:
3341 block_inds_to_transform = [(block_index, block_id) for block_index, block_id in enumerate(
3342 fexo.get_elem_blk_ids()) if block_id in blocks_to_transform]
3344 # If the exodus file has element variables, we are only keeping the ones
3345 # corresponding to the blocks we are keeping
3346 if exo_in_mem:
3347 has_elem_vars = (fexo.elem_vars is not None) and len(node_variables_to_transform) > 0
3348 else:
3349 try:
3350 has_elem_vars = fexo.num_elem_variables > 0 and len(elem_variables_to_transform) > 0
3351 except ExodusError:
3352 has_elem_vars = False
3353 if has_elem_vars:
3354 fexo_out.elem_vars = []
3356 # Now we loop through each block to reduce the elements to faces
3357 for block_index, block_id in block_inds_to_transform:
3358 # Determine the element type
3359 if exo_in_mem:
3360 elem_type = copy.deepcopy(fexo.blocks[block_index].elem_type.lower())
3361 connectivity = copy.deepcopy(fexo.blocks[block_index].connectivity)
3362 else:
3363 elem_type = copy.deepcopy(fexo.get_elem_blk_info(block_id)[0].lower())
3364 connectivity = copy.deepcopy(fexo.get_elem_connectivity(block_id).data)
3365 if verbose:
3366 print('Analyzing Block {}, type {}'.format(block_id, elem_type))
3367 if elem_type in reduced_element_types:
3368 block_type_ind = reduced_element_types.index(elem_type)
3369 # If the element type IS in our list of types to reduce, then we
3370 # will reduce to only surface nodes and elements.
3371 # Here we are computing all the faces in the block and creating a
3372 # face connectivity matrix (nfaces x nnodes per face) from an
3373 # element connectivity matrix (nelems x nnodes per element). It is
3374 # a complicated one-liner but boy is it fast. The basic steps are:
3375 # 1. Pass the element reduce_element_face vector in as the column
3376 # indices of the element connectivity matrix, which selects the face
3377 # nodes of all faces, making each row of the block connectivity
3378 # matrix contain all faces
3379 # block.connectivity[:,reduced_element_face_connectivity[block_type_ind]]
3380 # 2. Since each row now contains a number of faces, we want to
3381 # split it out into each face using the reshape command with the
3382 # number of nodes per face. After we do this, we have a connectivity
3383 # matrix where each row is a face of the element and the columns are
3384 # the nodes in that face.
3385 face_connectivity = connectivity[:, reduced_element_face_connectivity[block_type_ind]].reshape(
3386 (-1, reduced_element_nodes_per_face[block_type_ind]))
3387 # We will create a list that ties each face back to the original
3388 # element. This is used to map the element variables that were
3389 # defined for the volume elements to the newly created surface
3390 # elements. Basically this vector will be increasing by 1, with
3391 # each number repeated for the number of faces in the element. For
3392 # example if the elements are tet4s with 4 faces per element, the
3393 # vector will look like
3394 # [0,0,0,0, 1,1,1,1, 2,2,2,2, 3,3,3,3, 4,4,4,4, ... ]
3395 face_original_elements = np.repeat(np.arange(connectivity.shape[0]), len(
3396 reduced_element_face_connectivity[block_type_ind]) / reduced_element_nodes_per_face[block_type_ind])
3397 # To find EXTERIOR faces, we want to look for faces that only exist
3398 # once. If two elements have the same face, then that face is
3399 # necessarily on the interior of the block. Because the ordering
3400 # of the nodes is somewhat arbitrary, we first sort the nodes so
3401 # they are in ascending order for each face, then we can directly
3402 # compare the rows because two rows consisting of the same nodes
3403 # will look the same if they are sorted. This gives us the unique
3404 # rows (but we don't really care about them because they are
3405 # sorted, so we ditch the first output argument), the unique row
3406 # indices, which is a list of the first time that index appeared
3407 # for each unique entry, and a unique row mapping, which is the
3408 # location each row in the original face connectivity matrix ends
3409 # up in the unique set of rows (see help unique for more info)
3410 (unique_rows, unique_row_indices, unique_inverse, unique_counts) = np.unique(np.sort(
3411 face_connectivity, axis=1), return_index=True, return_inverse=True, return_counts=True, axis=0)
3412 # We can compute the values that are duplicates by using the optional
3413 # output argumenst of unique. The counts variable will enable us to
3414 # see which values were repeated, and the unique_inverse will allow us
3415 # to reconstruct the original array from the unique array.
3416 original_unique_counts = unique_counts[unique_inverse]
3417 # The non-duplicated faces are then the ones that only exist once in
3418 # the array
3419 nondupe_faces = original_unique_counts == 1
3420 # We then reduce the face connectivity matrix, keeping only the faces
3421 # that are not duplicated.
3422 face_connectivity = face_connectivity[nondupe_faces, :]
3423 # We also reduce our vector that ties the faces back to the original
3424 # values to keep it consistent with the face connectivity matrix
3425 face_original_elements = face_original_elements[nondupe_faces]
3426 # Now we need to also reduce the element variables in the same way,
3427 # mapping the original elements variables to the faces that correspond
3428 # to that element.
3429 new_elem_type = reduced_element_substitute_type[block_type_ind]
3430 else:
3431 # If we don't reduce,
3432 if verbose:
3433 print(' No Reduction')
3434 face_connectivity = connectivity
3435 face_original_elements = np.arange(connectivity.shape[0])
3436 new_elem_type = elem_type
3437 # We then update the connectivity in the block to make it match our new
3438 # face connectivity
3439 blk_dict = {}
3440 blk_dict['id'] = block_id
3441 blk_dict['connectivity'] = face_connectivity
3442 blk_dict['elem_type'] = new_elem_type
3443 blk_dict['name'] = 'block_{:}'.format(block_id)
3444 blk_dict['status'] = 1
3445 blk_dict['attributes'] = None
3446 fexo_out.blocks.append(subfield(**blk_dict))
3447 if verbose:
3448 print(' Converted to {}'.format(blk_dict['elem_type']))
3449 # The nodes in the block are now all the unique nodes in the face
3450 # connectivity matrix, so we make sure to keep those nodes in the final
3451 # set of nodes.
3452 block_nodes = np.unique(face_connectivity)
3453 keep_nodes[block_nodes] = True
3455 # Update element variables if necessary
3456 if has_elem_vars:
3457 var_names = []
3458 var_data = []
3459 if exo_in_mem:
3460 # Select the correct variable
3461 face_elem_vars = fexo.elem_vars[[var.id for var in fexo.elem_vars].index(block_id)]
3462 for elem_index, elem_var_data in enumerate(face_elem_vars.elem_var_data):
3463 # Check if we need to skip it
3464 if elem_index not in elem_variables_to_transform:
3465 print(' Skipping Element Variable {:}'.format(elem_index))
3466 continue
3467 # If the variable doesn't exist it will be None (for example a
3468 # volume quantity on a surface element). We only operate if
3469 # the variable exists (is not None)
3470 if elem_var_data is not None:
3471 # We use the face's original elements vector to map the
3472 # element data to the face connectivity matrix
3473 if verbose:
3474 print(' Processing Element Variable {:} for block {:}'.format(
3475 elem_var_data.name, block_id))
3476 var_data.append(elem_var_data.data[:, face_original_elements])
3477 var_names.append(elem_var_data.name)
3478 else:
3479 if verbose:
3480 print(' Element Variable {:} not defined for block {:}'.format(
3481 elem_index, block_id))
3482 var_data.append(None)
3483 var_names.append(None)
3484 else:
3485 # Get the element truth table
3486 elem_truth_table = fexo.get_elem_variable_table()
3487 # Get the element variable names
3488 elem_var_names = fexo.get_elem_variable_names()
3489 # Loop through and pick element variables that exist in the truth table
3490 for elem_index, elem_name in enumerate(elem_var_names):
3491 # Check if we need to skip it
3492 if elem_index not in elem_variables_to_transform:
3493 if verbose:
3494 print(' Skipping Element Variable {:}'.format(elem_name))
3495 continue
3496 if elem_truth_table[block_index, elem_index] == 0:
3497 if verbose:
3498 print(' Element Variable {:} not defined for block {:}'.format(
3499 elem_name, block_id))
3500 var_data.append(None)
3501 var_names.append(None)
3502 else:
3503 # We use the face's original elements vector to map the
3504 # element data to the face connectivity matrix
3505 if verbose:
3506 print(' Processing Element Variable {:} for block {:}'.format(
3507 elem_name, block_id))
3508 var_data.append(fexo.get_elem_variable_values(
3509 block_id, elem_index).data[:, face_original_elements])
3510 var_names.append(elem_name)
3511 # Now add the element variable
3512 if verbose:
3513 print(' Adding Element Variables to output exodus file')
3514 blk_dict = {}
3515 blk_dict['id'] = block_id
3516 blk_dict['elem_var_data'] = []
3517 for name, data in zip(var_names, var_data):
3518 if name is None or data is None:
3519 blk_dict['elem_var_data'].append(None)
3520 else:
3521 elem_var_dict = {}
3522 elem_var_dict['name'] = name
3523 elem_var_dict['data'] = data
3524 blk_dict['elem_var_data'].append(subfield(**elem_var_dict))
3525 fexo_out.elem_vars.append(subfield(**blk_dict))
3527 # We have reduced the face connectivity, so now we need to reduce the nodes
3528 # as well. We create a node map vector that maps the initial node indices
3529 # to the new node indices
3530 # i.e. the new_node_index = node_map[old_node_index]
3531 # We will need to do this because once we reduce the nodes, we will need to
3532 # re-update the element connectivity matrices with the updated node numbers
3533 # since the element connectivity matrix refers to the node index, not a
3534 # node number. By eliminating nodes, we are effectively renumbering the
3535 # nodes.
3536 if verbose:
3537 print('Analyzing Nodes')
3538 if exo_in_mem:
3539 original_nodes = np.arange(len(keep_nodes)) + \
3540 1 if fexo.nodes.node_num_map is None else fexo.nodes.node_num_map
3541 else:
3542 original_nodes = fexo.get_node_num_map()
3543 fexo_out.nodes.node_num_map = original_nodes[keep_nodes]
3544 node_map = np.zeros(keep_nodes.shape, dtype=int)
3545 node_map[keep_nodes] = np.arange(sum(keep_nodes))
3546 # Reduce the nodes, keeping only the ones where keep_nodes is true
3547 if verbose:
3548 print(' Reducing Coordinates')
3549 if exo_in_mem:
3550 fexo_out.nodes.coordinates = copy.deepcopy(fexo.nodes.coordinates[keep_nodes, :])
3551 else:
3552 fexo_out.nodes.coordinates = fexo.get_coords()[:, keep_nodes].T
3553 # We do the same thing to the nodal variables if they exist
3554 if verbose:
3555 print(' Reducing Nodal Variables')
3556 if exo_in_mem:
3557 if len(node_variables_to_transform) > 0:
3558 fexo_out.nodal_vars = []
3559 for node_var_index, node_var in enumerate(fexo.nodal_vars):
3560 if node_var_index not in node_variables_to_transform:
3561 if verbose:
3562 print(' Skipping Nodal Variable {:}'.format(node_var.name))
3563 continue
3564 node_dict = {}
3565 node_dict['name'] = node_var.name
3566 node_dict['data'] = copy.deepcopy(node_var.data[:, keep_nodes])
3567 fexo_out.nodal_vars.append(subfield(**node_dict))
3568 else:
3569 try:
3570 if len(node_variables_to_transform) > 0:
3571 fexo_out.nodal_vars = []
3572 node_vars = fexo.get_node_variable_names()
3573 for node_var_index, node_var_name in enumerate(node_vars):
3574 if node_var_index not in node_variables_to_transform:
3575 if verbose:
3576 print(' Skipping Nodal Variable {:}'.format(node_var_name))
3577 continue
3578 if verbose:
3579 print(' Variable {:}: {:}'.format(node_var_index, node_var_name))
3580 node_dict = {}
3581 node_dict['name'] = node_var_name
3582 node_dict['data'] = np.empty((len(fexo_out.time), sum(keep_nodes)))
3583 for i, time in enumerate(fexo_out.time):
3584 if verbose:
3585 print(' Step {:}'.format(i))
3586 node_dict['data'][i, :] = fexo.get_node_variable_values(
3587 node_var_index, step=i).data[keep_nodes]
3588 fexo_out.nodal_vars.append(subfield(**node_dict))
3589 except ExodusError:
3590 fexo_out.nodal_vars = None
3592 # Finally we go through each element block and map the old node indices to
3593 # the new node indices
3594 for block in fexo_out.blocks:
3595 block.connectivity = node_map[block.connectivity]
3597 # Return our new fexo
3598 return fexo_out
3601def extract_sharp_edges(exo, edge_threshold=60, **kwargs):
3602 dot_threshold = np.cos(edge_threshold * np.pi / 180)
3603 surface_fexo = reduce_exodus_to_surfaces(exo, **kwargs)
3604 edge_fexo = copy.deepcopy(surface_fexo)
3605 keep_nodes = np.zeros(surface_fexo.nodes.coordinates.shape[0], dtype=bool)
3606 for block_ind, block in enumerate(surface_fexo.blocks):
3607 connectivity = block.connectivity
3608 if connectivity.shape[-1] < 3 or 'bar' in block.elem_type.lower() or 'beam' in block.elem_type.lower():
3609 edge_fexo.blocks[block_ind].id = None
3610 continue
3611 print('Reducing Block {:}'.format(block.id))
3612 edge_connectivity = np.concatenate((connectivity[..., np.newaxis],
3613 np.roll(connectivity, -1, -1)[..., np.newaxis]), -1).reshape(-1, 2)
3614 edge_original_elements = np.repeat(np.arange(connectivity.shape[0]), connectivity.shape[-1])
3615 (unique_rows, unique_row_indices, unique_inverse, unique_counts) = np.unique(np.sort(
3616 edge_connectivity, axis=1), return_index=True, return_inverse=True, return_counts=True, axis=0)
3617 # Now we go through each edge and compute the angle between the two original faces
3618 keep_edges = np.zeros(unique_rows.shape[0], dtype=bool)
3619 for i, (edge, count) in enumerate(zip(unique_rows, unique_counts)):
3620 if count != 2:
3621 keep_edges[i] = True
3622 else:
3623 # Find faces that the edges are in
3624 edge_indices = np.nonzero(unique_inverse == i)[0]
3625 face_indices = edge_original_elements[edge_indices]
3626 face_connectivities = connectivity[face_indices]
3627 node_coords = surface_fexo.nodes.coordinates[face_connectivities[:, :3], :]
3628 vectors_a = node_coords[:, 2] - node_coords[:, 0]
3629 vectors_b = node_coords[:, 1] - node_coords[:, 0]
3630 normals = np.cross(vectors_a, vectors_b)
3631 normals /= np.linalg.norm(normals, axis=-1)[:, np.newaxis]
3632 if np.dot(*normals) < dot_threshold:
3633 keep_edges[i] = True
3634 print(' Kept {:} of {:} edges'.format(keep_edges.sum(), keep_edges.size))
3635 edge_connectivity = unique_rows[keep_edges]
3636 edge_fexo.blocks[block_ind].connectivity = edge_connectivity
3637 edge_fexo.blocks[block_ind].elem_type = 'bar'
3638 block_nodes = np.unique(edge_connectivity)
3639 keep_nodes[block_nodes] = True
3640 node_map = np.zeros(keep_nodes.shape, dtype=int)
3641 node_map[keep_nodes] = np.arange(sum(keep_nodes))
3642 edge_fexo.nodes.coordinates = edge_fexo.nodes.coordinates[keep_nodes, :]
3643 for node_var_index, node_var in enumerate(edge_fexo.nodal_vars):
3644 node_var.data = node_var.data[:, keep_nodes]
3645 edge_fexo.elem_vars = None
3646 edge_fexo.blocks = [block for block in edge_fexo.blocks if block.id is not None]
3647 for block in edge_fexo.blocks:
3648 block.connectivity = node_map[block.connectivity]
3649 edge_fexo.nodes.node_num_map = surface_fexo.nodes.node_num_map[keep_nodes]
3650 return edge_fexo
3652def read_sierra_matlab_matrix_file(file):
3653 """
3654 Reads a matrix.m file from Sierra Structural Dynamics.
3656 If the KAA or MAA outputs are used, this function can read those matrices.
3658 Parameters
3659 ----------
3660 file : str
3661 The name of the file to read.
3663 Returns
3664 -------
3665 matrix : np.ndarray
3666 A NumPy array containing the matrix.
3667 """
3668 with open(file,'r') as f:
3669 f.readline()
3670 num_rows = int(f.readline().replace(';','').split()[-1])
3671 matrix = np.zeros((num_rows,num_rows))
3672 f.readline()
3673 while True:
3674 try:
3675 row,col,data = f.readline().replace(';','').split()
3676 matrix[int(row)-1,int(col)-1] = float(data)
3677 if int(row) != int(col):
3678 matrix[int(col)-1,int(row)-1] = float(data)
3679 except ValueError:
3680 break
3681 return matrix
3683def read_sierra_matlab_map_file(file,return_inverse=False):
3684 """
3685 Reads a map file from Sierra Structural Dynamics.
3687 If the MFILE output is used, this can read the global node ID and aset map arrays that
3688 are output from Sierra SD.
3690 Parameters
3691 ----------
3692 file : str
3693 File containing the map data.
3694 return_inverse : bool, optional
3695 If True, will return both the forward map and the inverse map. The default is False.
3697 Returns
3698 -------
3699 array : np.ndarray
3700 Returns the same map as in the file.
3701 array_inv : np.ndarray
3702 Returns the inverse map of the file. Only returned if return_inverse is True.
3704 """
3705 with open(file,'r') as f:
3706 f.readline()
3707 f.readline()
3708 shape = [int(v) for v in f.readline().strip().split('(')[1][:-2].split(',')]
3709 array = -1*np.ones(shape,dtype=int)
3710 while True:
3711 try:
3712 lhs,rhs = f.readline().split('=')
3713 index = int(lhs.strip()[2:-1])-1
3714 value = int(rhs.strip()[:-1])
3715 array[index] = value
3716 except ValueError:
3717 break
3718 array = array.flatten()
3719 if not return_inverse:
3720 return array
3721 array_inv = -1*np.ones(array.max()+1,dtype=int)
3722 for index,value in enumerate(array):
3723 if value == -1:
3724 continue
3725 array_inv[value] = index
3726 return array,array_inv