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