From 0ef835c45aaf8c3d6632b40f84c0e243b622eac2 Mon Sep 17 00:00:00 2001 From: jfhuang Date: Fri, 10 Jan 2025 14:06:54 -0800 Subject: [PATCH 1/6] Sync to include recent updates --- convert.py | 1252 ++++++++++++++++++----------------- template/motion_setting.csv | 2 +- 2 files changed, 631 insertions(+), 623 deletions(-) diff --git a/convert.py b/convert.py index 16428e7..980323c 100755 --- a/convert.py +++ b/convert.py @@ -1,19 +1,19 @@ #!/usr/bin/env python3 # from genericpath import exists -import os -import shutil -# import sys -import argparse -import h5py -import math - -import scipy -try: - cumtrapz = scipy.integrate.cumulative_trapezoid -except AttributeError: - cumtrapz = scipy.integrate.cumtrapz -from scipy.interpolate import interpn +import os +import shutil +# import sys +import argparse +import h5py +import math + +import scipy +try: + cumtrapz = scipy.integrate.cumulative_trapezoid +except AttributeError: + cumtrapz = scipy.integrate.cumtrapz +from scipy.interpolate import interpn import numpy as np @@ -27,103 +27,111 @@ from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection from mpi4py import MPI -import functools -print = functools.partial(print, flush=True) # Don't buffer print -import hdf5plugin # Use this when SW4 output uses ZFP compression, can be installed with "pip install hdf5plugin" - - -def build_arg_parser(): - parser = argparse.ArgumentParser() - parser.add_argument("-c", "--csv", help="full path to the CSV setting file", default="") - parser.add_argument("-d", "--drm", help="full path to the DRM file with node coordinates", default="") - parser.add_argument("-h5", "--hdf5", help="full path to the hdf5 file with node coordinates", default="") - parser.add_argument("-e", "--ssi", "--essi", dest="ssi", - help="full path to the SW4 SSI output file", default="") - parser.add_argument("-t", "--template", - help="full path to the SSI template file with node coordinates", default="") - parser.add_argument("-p", "--plotonly", help="only generate plots of the input nodes", action="store_true") - parser.add_argument("-r", "--reference", help="reference node coordinate offset, default 0 0 0", - nargs='+', type=float) - # parser.add_argument("-s", "--steprange", help="timestep range, default 0 total_steps", nargs='+', type=int) - parser.add_argument("-s", "--timerange", - help="time range, will return all steps after the lower limit for equal upper and lower limit", - nargs='+', type=int) - parser.add_argument("-R", "--rotateanlge", help="rotate angle for node coordinate and motion: [0, 360)", - type=float) - parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") - parser.add_argument("-P", "--savepath", help="full path for saving the result files", default="") - parser.add_argument( - "-o", - "--output-format", - dest="output_format", - choices=("point", "opensees", "essi"), - help="output format override: point motions, OpenSees DRM, or ESSI template output", - default=None, - ) - parser.add_argument("-z", "--zeroMotionDir", - help="direction for zeroing out motion and enforce same motion across nodes in that direction: None(default), x, y, z", - default="") - return parser - - -def resolve_output_mode(input_kind, requested_output_mode): - defaults = { - 'drm': 'opensees', - 'h5': 'point', - 'csv': 'point', - 'template': 'essi', - } - supported = { - 'drm': {'point', 'opensees', 'essi'}, - 'h5': {'point'}, - 'csv': {'point'}, - 'template': {'point', 'opensees', 'essi'}, - } - - if requested_output_mode is None: - return defaults[input_kind] - - output_mode = requested_output_mode.lower() - if output_mode not in supported[input_kind]: - supported_modes = ', '.join(sorted(supported[input_kind])) - raise ValueError( - f'Output format "{output_mode}" is not supported for {input_kind} input. ' - f'Supported formats: {supported_modes}.' - ) - return output_mode - - -def get_output_filename(input_kind, output_mode, save_path, source_fname, explicit_output_mode=False): - if output_mode == 'opensees': - return os.path.join(save_path, 'OpenSeesDRMinput.h5drm') - - if output_mode == 'point': - if input_kind == 'h5': - output_stem = os.path.splitext(os.path.basename(source_fname))[0] - return os.path.join(save_path, f'{output_stem}_motion.h5') - if input_kind == 'csv': - return os.path.join(save_path, 'csvNodeMotion.h5') - if input_kind == 'drm': - return os.path.join(save_path, 'drmNodeMotion.h5') - if input_kind == 'template': - output_stem = os.path.splitext(os.path.basename(source_fname))[0] - return os.path.join(save_path, f'{output_stem}_motion.h5') - - if output_mode == 'essi': - if input_kind == 'template' and not explicit_output_mode: - return source_fname - return os.path.join(save_path, os.path.basename(source_fname)) - - raise ValueError(f'Unsupported output format: {output_mode}') - - -def prepare_essi_output_file(source_fname, output_fname): - if os.path.abspath(source_fname) != os.path.abspath(output_fname): - shutil.copyfile(source_fname, output_fname) - -# from scipy.signal import butter,filtfilt -# def butter_lowpass_filter(data, cutoff, nyq, order): -# normal_cutoff = cutoff / nyq +import functools +print = functools.partial(print, flush=True) # Don't buffer print +import hdf5plugin # Use this when SW4 output uses ZFP compression, can be installed with "pip install hdf5plugin" + + +def build_arg_parser(): + parser = argparse.ArgumentParser() + parser.add_argument("-c", "--csv", help="full path to the CSV setting file", default="") + parser.add_argument("-d", "--drm", help="full path to the DRM file with node coordinates", default="") + parser.add_argument("-h5", "--hdf5", help="full path to the hdf5 file with node coordinates", default="") + parser.add_argument("-e", "--ssi", "--essi", dest="ssi", + help="full path to the SW4 SSI output file", default="") + parser.add_argument("-t", "--template", + help="full path to the SSI template file with node coordinates", default="") + parser.add_argument("-p", "--plotonly", help="only generate plots of the input nodes", action="store_true") + parser.add_argument("-r", "--reference", help="reference node coordinate offset, default 0 0 0", + nargs='+', type=float) + # parser.add_argument("-s", "--steprange", help="timestep range, default 0 total_steps", nargs='+', type=int) + parser.add_argument("-s", "--timerange", + help="time range, will return all steps after the lower limit for equal upper and lower limit", + nargs='+', type=int) + parser.add_argument("-R", "--rotateanlge", help="rotate angle for node coordinate and motion: [0, 360)", + type=float) + parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") + parser.add_argument("-P", "--savepath", help="full path for saving the result files", default="") + parser.add_argument( + "-o", + "--output-format", + dest="output_format", + choices=("point", "opensees", "essi"), + help="output format override: point motions, OpenSees DRM, or ESSI template output", + default=None, + ) + parser.add_argument("-z", "--zeroMotionDir", + help="direction for zeroing out motion and enforce same motion across nodes in that direction: None(default), x, y, z", + default="") + return parser + + +def resolve_output_mode(input_kind, requested_output_mode): + defaults = { + 'drm': 'opensees', + 'h5': 'point', + 'csv': 'point', + 'template': 'essi', + } + supported = { + 'drm': {'point', 'opensees', 'essi'}, + 'h5': {'point'}, + 'csv': {'point'}, + 'template': {'point', 'opensees', 'essi'}, + } + + if requested_output_mode is None: + return defaults[input_kind] + + output_mode = requested_output_mode.lower() + if output_mode not in supported[input_kind]: + supported_modes = ', '.join(sorted(supported[input_kind])) + raise ValueError( + f'Output format "{output_mode}" is not supported for {input_kind} input. ' + f'Supported formats: {supported_modes}.' + ) + return output_mode + + +def get_output_filename(input_kind, output_mode, save_path, source_fname, explicit_output_mode=False): + if output_mode == 'opensees': + return os.path.join(save_path, 'OpenSeesDRMinput.h5drm') + + if output_mode == 'point': + if input_kind == 'h5': + output_stem = os.path.splitext(os.path.basename(source_fname))[0] + return os.path.join(save_path, f'{output_stem}_motion.h5') + if input_kind == 'csv': + return os.path.join(save_path, 'csvNodeMotion.h5') + if input_kind == 'drm': + return os.path.join(save_path, 'drmNodeMotion.h5') + if input_kind == 'template': + output_stem = os.path.splitext(os.path.basename(source_fname))[0] + return os.path.join(save_path, f'{output_stem}_motion.h5') + + if output_mode == 'essi': + if input_kind == 'template' and not explicit_output_mode: + return source_fname + return os.path.join(save_path, os.path.basename(source_fname)) + + raise ValueError(f'Unsupported output format: {output_mode}') + + +def prepare_essi_output_file(source_fname, output_fname): + if os.path.abspath(source_fname) != os.path.abspath(output_fname): + shutil.copyfile(source_fname, output_fname) + +# from scipy.signal import butter,filtfilt +# def butter_lowpass_filter(data, cutoff, nyq, order): +# normal_cutoff = cutoff / nyq +# # Get the filter coefficients +# b, a = butter(order, normal_cutoff, btype='low', analog=False) +# y = filtfilt(b, a, data, axis=0, method="gust") +# return y + +# from scipy.signal import butter,filtfilt +# def butter_lowpass_filter(data, cutoff, nyq, order): +# normal_cutoff = cutoff / nyq # # Get the filter coefficients # b, a = butter(order, normal_cutoff, btype='low', analog=False) # y = filtfilt(b, a, data, axis=0, method="gust") @@ -218,23 +226,23 @@ def plot_cube(save_path, cube_definition, x, y, z, view): ax.yaxis.set_tick_params(labelsize=lblsize) ax.xaxis.set_tick_params(labelsize=lblsize) - # ax.dist = 12 - #ax.set_aspect('equal') - if hasattr(ax, 'set_box_aspect'): - box_aspect = ( - max(x_plot_max - x_plot_min, 1e-12), - max(y_plot_max - y_plot_min, 1e-12), - max(z_plot_max - z_plot_min, 1e-12), - ) - try: - ax.set_box_aspect(box_aspect, zoom=0.5) - except TypeError: - try: - ax.set_box_aspect(box_aspect) - except TypeError: - ax.dist = 12 - else: - ax.dist = 12 + # ax.dist = 12 + #ax.set_aspect('equal') + if hasattr(ax, 'set_box_aspect'): + box_aspect = ( + max(x_plot_max - x_plot_min, 1e-12), + max(y_plot_max - y_plot_min, 1e-12), + max(z_plot_max - z_plot_min, 1e-12), + ) + try: + ax.set_box_aspect(box_aspect, zoom=0.5) + except TypeError: + try: + ax.set_box_aspect(box_aspect) + except TypeError: + ax.dist = 12 + else: + ax.dist = 12 ax.text(cube_definition[2][0], cube_definition[2][1], cube_definition[2][2]-z_len*.05, 'SW4-ESSI domain', fontsize=7) @@ -290,95 +298,95 @@ def plot_coords(essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, us plot_cube(save_path, cube_definition, user_essi_x, user_essi_y, user_essi_z, 'XZ') plot_cube(save_path, cube_definition, user_essi_x, user_essi_y, user_essi_z, 'XY') -def read_coord_drm(drm_filename, verbose): - if verbose: - print('Reading coordinates from input file [%s]' % drm_filename) - - # Get the coordinates from DRM file - with h5py.File(drm_filename, 'r') as drm_file: - coordinates = drm_file['Coordinates'] - isboundary = drm_file['Is Boundary Node'][:] - - if coordinates.ndim == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('Coordinates dataset length must be divisible by 3 for 1D input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('Coordinates dataset length must be divisible by 3 for single-column input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 3: - n_coord = coordinates.shape[0] - else: - raise ValueError('Coordinates dataset must have shape (3*n,), (3*n, 1), or (n, 3)') - - drm_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on - drm_y = np.zeros(n_coord, dtype='f8') - drm_z = np.zeros(n_coord, dtype='f8') - - if coordinates.ndim == 1: - for i in range(0, n_coord): - drm_x[i] = coordinates[i*3] - drm_y[i] = coordinates[i*3+1] - drm_z[i] = coordinates[i*3+2] - elif coordinates.shape[1] == 1: - for i in range(0, n_coord): - drm_x[i] = coordinates[i*3, 0] - drm_y[i] = coordinates[i*3+1, 0] - drm_z[i] = coordinates[i*3+2, 0] - else: # coordinates.shape[1] == 3 - drm_x[:] = coordinates[:,0] - drm_y[:] = coordinates[:,1] - drm_z[:] = coordinates[:,2] - - return drm_x, drm_y, drm_z, n_coord, isboundary - -def read_coord_h5(h5_filename, verbose): - if verbose: - print('Reading coordinates from input file [%s]' % h5_filename) - - # Get the coordinates from h5 file - with h5py.File(h5_filename, 'r') as h5_file: - coordinates = h5_file['coordinate'] - if coordinates.ndim == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('coordinate dataset length must be divisible by 3 for 1D input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('coordinate dataset length must be divisible by 3 for single-column input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 3: - n_coord = coordinates.shape[0] - else: - raise ValueError('coordinate dataset must have shape (3*n,), (3*n, 1), or (n, 3)') - - h5_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on - h5_y = np.zeros(n_coord, dtype='f8') - h5_z = np.zeros(n_coord, dtype='f8') - nodeTags = h5_file['nodeTag'][:] - # print('h5_x.dtype: ', h5_x.dtype) - - if coordinates.ndim == 1: - for i in range(0, n_coord): - h5_x[i] = coordinates[i*3] - h5_y[i] = coordinates[i*3+1] - h5_z[i] = coordinates[i*3+2] - elif coordinates.shape[1] == 1: - for i in range(0, n_coord): - h5_x[i] = coordinates[i*3, 0] - h5_y[i] = coordinates[i*3+1, 0] - h5_z[i] = coordinates[i*3+2, 0] - else: # coordinates.shape[1] == 3 - # h5_x = coordinates[:,0] # use float32 - # h5_y = coordinates[:,1] - # h5_z = coordinates[:,2] - h5_x[:] = coordinates[:,0] # use float64 - h5_y[:] = coordinates[:,1] - h5_z[:] = coordinates[:,2] - - # print('h5_x.dtype: ', h5_x.dtype) - return h5_x, h5_y, h5_z, n_coord, nodeTags +def read_coord_drm(drm_filename, verbose): + if verbose: + print('Reading coordinates from input file [%s]' % drm_filename) + + # Get the coordinates from DRM file + with h5py.File(drm_filename, 'r') as drm_file: + coordinates = drm_file['Coordinates'] + isboundary = drm_file['Is Boundary Node'][:] + + if coordinates.ndim == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('Coordinates dataset length must be divisible by 3 for 1D input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('Coordinates dataset length must be divisible by 3 for single-column input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 3: + n_coord = coordinates.shape[0] + else: + raise ValueError('Coordinates dataset must have shape (3*n,), (3*n, 1), or (n, 3)') + + drm_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on + drm_y = np.zeros(n_coord, dtype='f8') + drm_z = np.zeros(n_coord, dtype='f8') + + if coordinates.ndim == 1: + for i in range(0, n_coord): + drm_x[i] = coordinates[i*3] + drm_y[i] = coordinates[i*3+1] + drm_z[i] = coordinates[i*3+2] + elif coordinates.shape[1] == 1: + for i in range(0, n_coord): + drm_x[i] = coordinates[i*3, 0] + drm_y[i] = coordinates[i*3+1, 0] + drm_z[i] = coordinates[i*3+2, 0] + else: # coordinates.shape[1] == 3 + drm_x[:] = coordinates[:,0] + drm_y[:] = coordinates[:,1] + drm_z[:] = coordinates[:,2] + + return drm_x, drm_y, drm_z, n_coord, isboundary + +def read_coord_h5(h5_filename, verbose): + if verbose: + print('Reading coordinates from input file [%s]' % h5_filename) + + # Get the coordinates from h5 file + with h5py.File(h5_filename, 'r') as h5_file: + coordinates = h5_file['coordinate'] + if coordinates.ndim == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('coordinate dataset length must be divisible by 3 for 1D input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('coordinate dataset length must be divisible by 3 for single-column input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 3: + n_coord = coordinates.shape[0] + else: + raise ValueError('coordinate dataset must have shape (3*n,), (3*n, 1), or (n, 3)') + + h5_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on + h5_y = np.zeros(n_coord, dtype='f8') + h5_z = np.zeros(n_coord, dtype='f8') + nodeTags = h5_file['nodeTag'][:] + # print('h5_x.dtype: ', h5_x.dtype) + + if coordinates.ndim == 1: + for i in range(0, n_coord): + h5_x[i] = coordinates[i*3] + h5_y[i] = coordinates[i*3+1] + h5_z[i] = coordinates[i*3+2] + elif coordinates.shape[1] == 1: + for i in range(0, n_coord): + h5_x[i] = coordinates[i*3, 0] + h5_y[i] = coordinates[i*3+1, 0] + h5_z[i] = coordinates[i*3+2, 0] + else: # coordinates.shape[1] == 3 + # h5_x = coordinates[:,0] # use float32 + # h5_y = coordinates[:,1] + # h5_z = coordinates[:,2] + h5_x[:] = coordinates[:,0] # use float64 + h5_y[:] = coordinates[:,1] + h5_z[:] = coordinates[:,2] + + # print('h5_x.dtype: ', h5_x.dtype) + return h5_x, h5_y, h5_z, n_coord, nodeTags # changed ref coord as just offsets def convert_to_essi_coord(coord_sys, from_x, from_y, from_z, ref_essi_xyz): @@ -462,111 +470,111 @@ def rotate_coords_ops_xyplane(x, y, z, rotate_angle, ref_coord=[0,0,0]): return rotated_xyz[:,0], rotated_xyz[:,1], rotated_xyz[:,2] -def get_essi_meta(ssi_fname, verbose): - # Get parameter values from HDF5 data - with h5py.File(ssi_fname, 'r') as essiout: - h = essiout['ESSI xyz grid spacing'][0] - x0 = essiout['ESSI xyz origin'][0] - y0 = essiout['ESSI xyz origin'][1] - z0 = essiout['ESSI xyz origin'][2] - t0 = essiout['time start'][0] - dt = essiout['timestep'][0] - nt = essiout['vel_0 ijk layout'].shape[0] - nx = essiout['vel_0 ijk layout'].shape[1] - ny = essiout['vel_0 ijk layout'].shape[2] - nz = essiout['vel_0 ijk layout'].shape[3] - t1 = t0 + dt*(nt-1) - timeseq = np.linspace(t0, t1, nt) - # print('dt, t0, t1, timeseq =', dt, t0, t1, timeseq) - - bTopo = False - zmin, zmax = z0, z0+(nz-1)*h - if 'z coordinates' in essiout: - bTopo = True - zcoords = essiout['z coordinates'] - chunk_x = zcoords.chunks[0] if zcoords.chunks else nx - chunk_y = zcoords.chunks[1] if zcoords.chunks else ny - zmin = np.inf - zmax = -np.inf - for ix in range(0, nx, chunk_x): - for iy in range(0, ny, chunk_y): - upper = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), 0] - lower = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), -1] - zmin = min(zmin, float(np.min(upper))) - zmax = max(zmax, float(np.max(lower))) - - return x0, y0, z0, h, nx, ny, nz, nt, dt, timeseq, bTopo, zmin, zmax - - -def get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose): - ''' when topography is considered, get z array location based on its physical - z coordinates in the interpolated vertical profile between the upper and lower z - interfaces - ''' - with h5py.File(ssi_fname, 'r') as essiout: - h = essiout['ESSI xyz grid spacing'][0] - _, nx, ny, nz = essiout['vel_0 ijk layout'].shape - # print("essiout['vel_0 ijk layout'].shape: ", essiout['vel_0 ijk layout'].shape) - - if nx < 2 or ny < 2 or nz < 2: - raise ValueError('Topography-aware z interpolation requires at least 2 grid points in x, y, and z.') - - # *) x, y indices - # start - coord_x0 = np.floor(coord_x).astype(int) - coord_y0 = np.floor(coord_y).astype(int) - coord_x0 = np.clip(coord_x0, 0, nx-2) - coord_y0 = np.clip(coord_y0, 0, ny-2) - # end - coord_x1 = coord_x0 + 1 - coord_y1 = coord_y0 + 1 - coord_x1 = np.minimum(coord_x1, nx-1) - coord_y1 = np.minimum(coord_y1, ny-1) - # get distinct x, y values - x = np.array(sorted(set(coord_x0) | set(coord_x1)), dtype=int) - y = np.array(sorted(set(coord_y0) | set(coord_y1)), dtype=int) - x_idx = x - x[0] - y_idx = y - y[0] - - # *) interpolated z profile upper and lower z physical coordinates - zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,0] # upper interface - z_upper = zcrds[np.ix_(x_idx, y_idx)] - z_upper_interp = interpn((x, y), z_upper, np.c_[coord_x, coord_y]) - - zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,-1] # lower interface - z_lower = zcrds[np.ix_(x_idx, y_idx)] - z_lower_interp = interpn((x, y), z_lower, np.c_[coord_x, coord_y]) - - z_span = z_lower_interp - z_upper_interp - zero_span = np.isclose(z_span, 0.0) - if np.any(zero_span): - iz = np.flatnonzero(zero_span)[0] - raise ValueError( - f'Error getting z array location: ' - f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' - f'has zero vertical span in the topography profile!' - ) - - z_min = np.minimum(z_upper_interp, z_lower_interp) - z_max = np.maximum(z_upper_interp, z_lower_interp) - tol = max(abs(h) * 1e-6, 1e-8) - out_of_range = (user_essi_z < z_min - tol) | (user_essi_z > z_max + tol) - if np.any(out_of_range): - iz = np.flatnonzero(out_of_range)[0] - raise ValueError( - f'Error getting z array location: ' - f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' - f'not within SW4 domain!' - ) - - # The vertical profile is linear between the interpolated upper/lower interfaces. - coord_z = (user_essi_z - z_upper_interp) * (nz - 1) / z_span - coord_z = np.clip(coord_z, 0.0, nz - 1) - - if verbose: - print('coord_z consider topography:', coord_z) - - return coord_z +def get_essi_meta(ssi_fname, verbose): + # Get parameter values from HDF5 data + with h5py.File(ssi_fname, 'r') as essiout: + h = essiout['ESSI xyz grid spacing'][0] + x0 = essiout['ESSI xyz origin'][0] + y0 = essiout['ESSI xyz origin'][1] + z0 = essiout['ESSI xyz origin'][2] + t0 = essiout['time start'][0] + dt = essiout['timestep'][0] + nt = essiout['vel_0 ijk layout'].shape[0] + nx = essiout['vel_0 ijk layout'].shape[1] + ny = essiout['vel_0 ijk layout'].shape[2] + nz = essiout['vel_0 ijk layout'].shape[3] + t1 = t0 + dt*(nt-1) + timeseq = np.linspace(t0, t1, nt) + # print('dt, t0, t1, timeseq =', dt, t0, t1, timeseq) + + bTopo = False + zmin, zmax = z0, z0+(nz-1)*h + if 'z coordinates' in essiout: + bTopo = True + zcoords = essiout['z coordinates'] + chunk_x = zcoords.chunks[0] if zcoords.chunks else nx + chunk_y = zcoords.chunks[1] if zcoords.chunks else ny + zmin = np.inf + zmax = -np.inf + for ix in range(0, nx, chunk_x): + for iy in range(0, ny, chunk_y): + upper = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), 0] + lower = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), -1] + zmin = min(zmin, float(np.min(upper))) + zmax = max(zmax, float(np.max(lower))) + + return x0, y0, z0, h, nx, ny, nz, nt, dt, timeseq, bTopo, zmin, zmax + + +def get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose): + ''' when topography is considered, get z array location based on its physical + z coordinates in the interpolated vertical profile between the upper and lower z + interfaces + ''' + with h5py.File(ssi_fname, 'r') as essiout: + h = essiout['ESSI xyz grid spacing'][0] + _, nx, ny, nz = essiout['vel_0 ijk layout'].shape + # print("essiout['vel_0 ijk layout'].shape: ", essiout['vel_0 ijk layout'].shape) + + if nx < 2 or ny < 2 or nz < 2: + raise ValueError('Topography-aware z interpolation requires at least 2 grid points in x, y, and z.') + + # *) x, y indices + # start + coord_x0 = np.floor(coord_x).astype(int) + coord_y0 = np.floor(coord_y).astype(int) + coord_x0 = np.clip(coord_x0, 0, nx-2) + coord_y0 = np.clip(coord_y0, 0, ny-2) + # end + coord_x1 = coord_x0 + 1 + coord_y1 = coord_y0 + 1 + coord_x1 = np.minimum(coord_x1, nx-1) + coord_y1 = np.minimum(coord_y1, ny-1) + # get distinct x, y values + x = np.array(sorted(set(coord_x0) | set(coord_x1)), dtype=int) + y = np.array(sorted(set(coord_y0) | set(coord_y1)), dtype=int) + x_idx = x - x[0] + y_idx = y - y[0] + + # *) interpolated z profile upper and lower z physical coordinates + zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,0] # upper interface + z_upper = zcrds[np.ix_(x_idx, y_idx)] + z_upper_interp = interpn((x, y), z_upper, np.c_[coord_x, coord_y]) + + zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,-1] # lower interface + z_lower = zcrds[np.ix_(x_idx, y_idx)] + z_lower_interp = interpn((x, y), z_lower, np.c_[coord_x, coord_y]) + + z_span = z_lower_interp - z_upper_interp + zero_span = np.isclose(z_span, 0.0) + if np.any(zero_span): + iz = np.flatnonzero(zero_span)[0] + raise ValueError( + f'Error getting z array location: ' + f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' + f'has zero vertical span in the topography profile!' + ) + + z_min = np.minimum(z_upper_interp, z_lower_interp) + z_max = np.maximum(z_upper_interp, z_lower_interp) + tol = max(abs(h) * 1e-6, 1e-8) + out_of_range = (user_essi_z < z_min - tol) | (user_essi_z > z_max + tol) + if np.any(out_of_range): + iz = np.flatnonzero(out_of_range)[0] + raise ValueError( + f'Error getting z array location: ' + f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' + f'not within SW4 domain!' + ) + + # The vertical profile is linear between the interpolated upper/lower interfaces. + coord_z = (user_essi_z - z_upper_interp) * (nz - 1) / z_span + coord_z = np.clip(coord_z, 0.0, nz - 1) + + if verbose: + print('coord_z consider topography:', coord_z) + + return coord_z def get_essi_data_btw_step(ssi_fname, start, end, verbose): @@ -716,28 +724,28 @@ def write_to_hdf5_range_1d(h5_fname, gname, dname, data, mystart, myend): dset[mystart:myend] = data[:] h5file.close() -def write_to_hdf5_range_2d(h5_fname, gname, dname, data, mystart, myend): - h5file = h5py.File(h5_fname, 'r+') - if gname == '/': - dset = h5file[dname] +def write_to_hdf5_range_2d(h5_fname, gname, dname, data, mystart, myend): + h5file = h5py.File(h5_fname, 'r+') + if gname == '/': + dset = h5file[dname] else: grp = h5file[gname] dset = grp[dname] - #print('mystart=%d, myend=%d' %(mystart, myend)) - dset[mystart:myend,:] = data[:,:] - h5file.close() - -def get_flat_coord_range(coord_offset, ncoord, width=3): - return coord_offset * width, (coord_offset + ncoord) * width - -def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): - tstep = tsteprange.step - nstep = len(tsteprange) - dt = tstep*essi_dt - tend = max(nstep - 1, 0) * dt - - h5file = h5py.File(h5_fname, 'w') + #print('mystart=%d, myend=%d' %(mystart, myend)) + dset[mystart:myend,:] = data[:,:] + h5file.close() + +def get_flat_coord_range(coord_offset, ncoord, width=3): + return coord_offset * width, (coord_offset + ncoord) * width + +def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): + tstep = tsteprange.step + nstep = len(tsteprange) + dt = tstep*essi_dt + tend = max(nstep - 1, 0) * dt + + h5file = h5py.File(h5_fname, 'w') data_grp = h5file.create_group('DRM_Data') data_location = np.zeros(ncoord, dtype='i4') for i in range(0, ncoord): @@ -751,22 +759,22 @@ def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc dset = data_grp.create_dataset('displacement', (ncoord*3, nstep), dtype='f4') dset = data_grp.create_dataset('data_location', data=data_location, dtype='i4') dset = data_grp.create_dataset(extra_dname, (ncoord,), dtype='i4') - dset = data_grp.create_dataset('xyz', (ncoord, 3), dtype='f4') - - data_grp = h5file.create_group('DRM_Metadata') - dset = data_grp.create_dataset('dt', data=dt, dtype='f8') - dset = data_grp.create_dataset('tstart', data=0.0, dtype='f8') - dset = data_grp.create_dataset('tend', data=tend, dtype='f8') + dset = data_grp.create_dataset('xyz', (ncoord, 3), dtype='f4') + + data_grp = h5file.create_group('DRM_Metadata') + dset = data_grp.create_dataset('dt', data=dt, dtype='f8') + dset = data_grp.create_dataset('tstart', data=0.0, dtype='f8') + dset = data_grp.create_dataset('tend', data=tend, dtype='f8') h5file.close() -def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): - - h5file = h5py.File(h5_fname, 'w') - tstep = tsteprange.step - nstep = len(tsteprange) - dt = tstep*essi_dt - tend = max(nstep - 1, 0) * dt +def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): + + h5file = h5py.File(h5_fname, 'w') + tstep = tsteprange.step + nstep = len(tsteprange) + dt = tstep*essi_dt + tend = max(nstep - 1, 0) * dt print('Create HDF5 file with ', ncoord, ' coordinates, ', nstep, ' steps') if gen_vel: @@ -777,42 +785,42 @@ def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen dset = h5file.create_dataset('displacement', (ncoord*3, nstep), dtype='f4') dset = h5file.create_dataset(extra_dname, (ncoord,), dtype='i4') - dset = h5file.create_dataset('xyz', (ncoord, 3), dtype='f4') - - dset = h5file.create_dataset('dt', data=dt, dtype='f8') - dset = h5file.create_dataset('tstart', data=0.0, dtype='f8') - dset = h5file.create_dataset('tend', data=tend, dtype='f8') + dset = h5file.create_dataset('xyz', (ncoord, 3), dtype='f4') + + dset = h5file.create_dataset('dt', data=dt, dtype='f8') + dset = h5file.create_dataset('tstart', data=0.0, dtype='f8') + dset = h5file.create_dataset('tend', data=tend, dtype='f8') h5file.close() -def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, extra_dname): - h5file = h5py.File(h5_fname, 'r+') - - for dname in ('Velocity', 'Accelerations', 'Displacements', 'Time'): - if dname in h5file: - del h5file[dname] - - if gen_vel: - dset = h5file.create_dataset('Velocity', (ncoord*3, nstep), dtype='f4') +def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, extra_dname): + h5file = h5py.File(h5_fname, 'r+') + + for dname in ('Velocity', 'Accelerations', 'Displacements', 'Time'): + if dname in h5file: + del h5file[dname] + + if gen_vel: + dset = h5file.create_dataset('Velocity', (ncoord*3, nstep), dtype='f4') if gen_acc: dset = h5file.create_dataset('Accelerations', (ncoord*3, nstep), dtype='f4') if gen_dis: dset = h5file.create_dataset('Displacements', (ncoord*3, nstep), dtype='f4') - timeseq = np.arange(nstep, dtype='f8') * dt - - h5file.create_dataset('Time', data=timeseq, dtype='f8') + timeseq = np.arange(nstep, dtype='f8') * dt + + h5file.create_dataset('Time', data=timeseq, dtype='f8') h5file.close() -def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx, ny, nz): - # Due to allocation of neighbor nodes, the indices may exceed the original dataset - # boundary and thus need to be corrected so that the correct nearest chunk can be fetched later - x, y, z = min(x0, nx-1), min(y0, ny-1), min(z0, nz-1) - - val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) - #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) - return val +def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx, ny, nz): + # Due to allocation of neighbor nodes, the indices may exceed the original dataset + # boundary and thus need to be corrected so that the correct nearest chunk can be fetched later + x, y, z = min(x0, nx-1), min(y0, ny-1), min(z0, nz-1) + + val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) + #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) + return val def chunkid_to_start(cid, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z): #print('cid2:', cid, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z) @@ -947,24 +955,24 @@ def linear_interp(data_dict, x, y, z): result = ((c000 * (1-xd) + c100 * xd) * (1-yd) + (c010 * (1-xd) + c110 * xd) * yd) * (1-zd) + ((c001 * (1-xd) + c101 * xd) * (1-yd) + (c011 * (1-xd) + c111 * xd) * yd) * zd return result -def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_data, extra_dname, output_format): - if zeroMotionDir is None or pd.isna(zeroMotionDir): - zeroMotionDir = 'None' - elif isinstance(zeroMotionDir, str): - zeroMotionDir = zeroMotionDir.strip() - if zeroMotionDir == '' or zeroMotionDir.lower() in ('none', 'no'): - zeroMotionDir = 'None' - else: - zeroMotionDir = str(zeroMotionDir).strip() - - if zeroMotionDir != 'None' and zeroMotionDir.upper() not in ('X', 'Y', 'Z'): - raise ValueError(f'Unsupported zeroMotionDir value: {zeroMotionDir}') - - # Read ESSI metadata - essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, essi_nt, essi_dt, essi_timeseq, bTopo, zmin, zmax = get_essi_meta(ssi_fname, verbose) - essi_x_len_max = (essi_nx-1) * essi_h - essi_y_len_max = (essi_ny-1) * essi_h - essi_z_len_max = zmax - zmin +def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_data, extra_dname, output_format): + if zeroMotionDir is None or pd.isna(zeroMotionDir): + zeroMotionDir = 'None' + elif isinstance(zeroMotionDir, str): + zeroMotionDir = zeroMotionDir.strip() + if zeroMotionDir == '' or zeroMotionDir.lower() in ('none', 'no'): + zeroMotionDir = 'None' + else: + zeroMotionDir = str(zeroMotionDir).strip() + + if zeroMotionDir != 'None' and zeroMotionDir.upper() not in ('X', 'Y', 'Z'): + raise ValueError(f'Unsupported zeroMotionDir value: {zeroMotionDir}') + + # Read ESSI metadata + essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, essi_nt, essi_dt, essi_timeseq, bTopo, zmin, zmax = get_essi_meta(ssi_fname, verbose) + essi_x_len_max = (essi_nx-1) * essi_h + essi_y_len_max = (essi_ny-1) * essi_h + essi_z_len_max = zmax - zmin # Start and end time step if start_t > -1e-6 and end_t > -1e-6: @@ -1003,13 +1011,13 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # to get 2D motions for 2D models, modify input node crds to enforce same motion across the width # Note: this should be done before rotation, motion zero-out in the out-of-plane direction will be done later # TODO: here we use the middle crd in that direction by default - user_x, user_y, user_z = user_x0, user_y0, user_z0 - # print(f'zeroMotionDir={zeroMotionDir}') - if zeroMotionDir != 'None': - middleCrd = None - if zeroMotionDir.upper() == 'X': - middleCrd = (np.amin(user_x0) + np.amax(user_x0)) / 2. - user_x = np.full(user_x0.shape, middleCrd) + user_x, user_y, user_z = user_x0, user_y0, user_z0 + # print(f'zeroMotionDir={zeroMotionDir}') + if zeroMotionDir != 'None': + middleCrd = None + if zeroMotionDir.upper() == 'X': + middleCrd = (np.amin(user_x0) + np.amax(user_x0)) / 2. + user_x = np.full(user_x0.shape, middleCrd) elif zeroMotionDir.upper() == 'Y': middleCrd = (np.amin(user_y0) + np.amax(user_y0)) / 2. user_y = np.full(user_y0.shape, middleCrd) @@ -1051,26 +1059,26 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use print('Only generate the plots of input nodes') exit(0) - # Check if all node coordinates are within the sw4 domain - if np.min(user_essi_x) < essi_x0 or np.max(user_essi_x) > essi_x0+essi_x_len_max or \ - np.min(user_essi_y) < essi_y0 or np.max(user_essi_y) > essi_y0+essi_y_len_max or \ - np.min(user_essi_z) < zmin or np.max(user_essi_z) > zmax: - if mpi_rank == 0: - print('Error: all node coordinates (after rotation) should be within the sw4 domain for extracting the motion') - print('while:') - print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) - print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) - print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) - - debugfile = save_path + '/user_essi_xyz.npy' - print('\tcheck user_essi_xyz (after rotation) in file \'{}\''.format(debugfile)) - np.save(debugfile, np.c_[user_essi_x, user_essi_y, user_essi_z]) - raise SystemExit(1) + # Check if all node coordinates are within the sw4 domain + if np.min(user_essi_x) < essi_x0 or np.max(user_essi_x) > essi_x0+essi_x_len_max or \ + np.min(user_essi_y) < essi_y0 or np.max(user_essi_y) > essi_y0+essi_y_len_max or \ + np.min(user_essi_z) < zmin or np.max(user_essi_z) > zmax: + if mpi_rank == 0: + print('Error: all node coordinates (after rotation) should be within the sw4 domain for extracting the motion') + print('while:') + print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) + print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) + print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) + + debugfile = save_path + '/user_essi_xyz.npy' + print('\tcheck user_essi_xyz (after rotation) in file \'{}\''.format(debugfile)) + np.save(debugfile, np.c_[user_essi_x, user_essi_y, user_essi_z]) + raise SystemExit(1) - if verbose and mpi_rank == 0: - print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) - print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) - print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) + if verbose and mpi_rank == 0: + print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) + print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) + print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) # if mpi_rank == 0: # print('while user_essi_xyz (after rotation) is:\n', np.c_[user_essi_x, user_essi_y, user_essi_z]) @@ -1078,28 +1086,28 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Convert to array location (spacing is 1), floating-point coord_x = (user_essi_x - essi_x0) / essi_h - coord_y = (user_essi_y - essi_y0) / essi_h - if not bTopo: - coord_z = (user_essi_z - essi_z0) / essi_h - else: - coord_z = get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose) + coord_y = (user_essi_y - essi_y0) / essi_h + if not bTopo: + coord_z = (user_essi_z - essi_z0) / essi_h + else: + coord_z = get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose) - # Check if we actually need interpolation - # ghost_cell = 0 - # do_interp = True - do_interp = False - coord_x_rounded = np.rint(coord_x) - coord_y_rounded = np.rint(coord_y) - coord_z_rounded = np.rint(coord_z) - for nid in range(0, n_coord): - # if user_essi_x[nid] % essi_h != 0 or user_essi_y[nid] % essi_h != 0 or user_essi_z[nid] % essi_h != 0: - if not math.isclose(coord_x[nid], coord_x_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_y[nid], coord_y_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_z[nid], coord_z_rounded[nid], abs_tol=1e-8): - do_interp = True - # ghost_cell = 1 - if verbose and mpi_rank == 0: - print(f'user_essi_x[{nid}], user_essi_y[{nid}], user_essi_z[{nid}]: {user_essi_x[nid]:.4f}, {user_essi_y[nid]:.4f}, {user_essi_z[nid]:.4f}') - # print(f'coord_x[{nid}], coord_y[{nid}], coord_z[{nid}]: {coord_x[nid]:.4e}, {coord_y[nid]:.4e}, {coord_z[nid]:.4e}') - break + # Check if we actually need interpolation + # ghost_cell = 0 + # do_interp = True + do_interp = False + coord_x_rounded = np.rint(coord_x) + coord_y_rounded = np.rint(coord_y) + coord_z_rounded = np.rint(coord_z) + for nid in range(0, n_coord): + # if user_essi_x[nid] % essi_h != 0 or user_essi_y[nid] % essi_h != 0 or user_essi_z[nid] % essi_h != 0: + if not math.isclose(coord_x[nid], coord_x_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_y[nid], coord_y_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_z[nid], coord_z_rounded[nid], abs_tol=1e-8): + do_interp = True + # ghost_cell = 1 + if verbose and mpi_rank == 0: + print(f'user_essi_x[{nid}], user_essi_y[{nid}], user_essi_z[{nid}]: {user_essi_x[nid]:.4f}, {user_essi_y[nid]:.4f}, {user_essi_z[nid]:.4f}') + # print(f'coord_x[{nid}], coord_y[{nid}], coord_z[{nid}]: {coord_x[nid]:.4e}, {coord_y[nid]:.4e}, {coord_z[nid]:.4e}') + break if mpi_rank == 0: if do_interp: print('Use spline interpolation.') @@ -1446,7 +1454,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Write coordinates and boundary nodes (file created previously), in serial with baton passing comm.Barrier() - if output_format == "opensees": + if output_format == "opensees": if mpi_rank == 0: create_hdf5_opensees(output_fname, n_coord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname) @@ -1478,7 +1486,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if mpi_rank != mpi_size-1: comm.send(my_ncoord, dest=mpi_rank+1, tag=11) - elif output_format == "point": + elif output_format == "point": if mpi_rank == 0: create_hdf5_csv(output_fname, n_coord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname) @@ -1510,30 +1518,30 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if mpi_rank != mpi_size-1: comm.send(my_ncoord, dest=mpi_rank+1, tag=11) - elif output_format == "essi": - coords_start, coords_end = get_flat_coord_range(my_offset, my_ncoord[0]) - if mpi_rank == 0: - create_hdf5_essi(output_fname, n_coord, nsteps, dt, gen_vel, gen_acc, gen_dis, extra_dname) - # Write to the template file - if my_ncoord[0] > 0: - write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) - if gen_acc: - write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) - if gen_dis: + elif output_format == "essi": + coords_start, coords_end = get_flat_coord_range(my_offset, my_ncoord[0]) + if mpi_rank == 0: + create_hdf5_essi(output_fname, n_coord, nsteps, dt, gen_vel, gen_acc, gen_dis, extra_dname) + # Write to the template file + if my_ncoord[0] > 0: + write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + if gen_acc: + write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) + if gen_dis: write_to_hdf5_range(output_fname, '/', 'Displacements', output_dis_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_vel: write_to_hdf5_range(output_fname, '/', 'Velocity', output_vel_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if mpi_size > 1: comm.send(my_ncoord, dest=1, tag=111) - else: - data = comm.recv(source=mpi_rank-1, tag=111) - if my_ncoord[0] > 0: - write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) - if gen_acc: - write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) + else: + data = comm.recv(source=mpi_rank-1, tag=111) + if my_ncoord[0] > 0: + write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + if gen_acc: + write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: write_to_hdf5_range(output_fname, '/', 'Displacements', output_dis_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_vel: @@ -1550,81 +1558,81 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use print('Rank', mpi_rank, 'Finished writing data') return -def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input DRM [%s]' %drm_fname) - print('Input ESSI [%s]' %ssi_fname) - - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('drm', requested_output_mode) - - # original unrotated node coordinates - user_x0, user_y0, user_z0, n_coord, isboundary = read_coord_drm(drm_fname, verbose) - isboundary = np.asarray(isboundary).reshape(-1) - if verbose and mpi_rank == 0: - print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x0[0], user_y0[0], user_z0[0], user_x0[-1], user_y0[-1], user_z0[-1])) - print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x0), np.max(user_x0), np.min(user_y0), np.max(user_y0), np.min(user_z0), np.max(user_z0)) ) - - if output_format == 'opensees': - gen_vel = False - gen_dis = True - gen_acc = True - extra_dname = 'internal' - else: - gen_vel = True - gen_dis = True - gen_acc = True - extra_dname = 'Is Boundary Node' - - output_fname = get_output_filename( - 'drm', output_format, save_path, drm_fname, explicit_output_mode=requested_output_mode is not None - ) - if output_format == 'essi': - if mpi_rank == 0: - prepare_essi_output_file(drm_fname, output_fname) - MPI.COMM_WORLD.Barrier() - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, isboundary, extra_dname, output_format) - - return - -def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input h5 [%s]' %h5_fname) - print('Input ESSI [%s]' %ssi_fname) - - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('h5', requested_output_mode) - gen_vel = True - gen_dis = True - gen_acc = True +def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input DRM [%s]' %drm_fname) + print('Input ESSI [%s]' %ssi_fname) + + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('drm', requested_output_mode) + + # original unrotated node coordinates + user_x0, user_y0, user_z0, n_coord, isboundary = read_coord_drm(drm_fname, verbose) + isboundary = np.asarray(isboundary).reshape(-1) + if verbose and mpi_rank == 0: + print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x0[0], user_y0[0], user_z0[0], user_x0[-1], user_y0[-1], user_z0[-1])) + print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x0), np.max(user_x0), np.min(user_y0), np.max(user_y0), np.min(user_z0), np.max(user_z0)) ) + + if output_format == 'opensees': + gen_vel = False + gen_dis = True + gen_acc = True + extra_dname = 'internal' + else: + gen_vel = True + gen_dis = True + gen_acc = True + extra_dname = 'Is Boundary Node' + + output_fname = get_output_filename( + 'drm', output_format, save_path, drm_fname, explicit_output_mode=requested_output_mode is not None + ) + if output_format == 'essi': + if mpi_rank == 0: + prepare_essi_output_file(drm_fname, output_fname) + MPI.COMM_WORLD.Barrier() + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, isboundary, extra_dname, output_format) + + return + +def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input h5 [%s]' %h5_fname) + print('Input ESSI [%s]' %ssi_fname) + + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('h5', requested_output_mode) + gen_vel = True + gen_dis = True + gen_acc = True extra_dname = 'nodeTag' # original unrotated node coordinates user_x0, user_y0, user_z0, n_coord, node_tags = read_coord_h5(h5_fname, verbose) n_coord = len(node_tags) - if mpi_rank == 0: - print('Generating motions for %i nodes...' % (n_coord)) - - output_fname = get_output_filename( - 'h5', output_format, save_path, h5_fname, explicit_output_mode=requested_output_mode is not None - ) - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) - - return - -def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input CSV [%s]' %csv_fname) - print('Input ESSI [%s]' %ssi_fname) - - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('csv', requested_output_mode) - gen_vel = True - gen_dis = True - gen_acc = True + if mpi_rank == 0: + print('Generating motions for %i nodes...' % (n_coord)) + + output_fname = get_output_filename( + 'h5', output_format, save_path, h5_fname, explicit_output_mode=requested_output_mode is not None + ) + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) + + return + +def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input CSV [%s]' %csv_fname) + print('Input ESSI [%s]' %ssi_fname) + + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('csv', requested_output_mode) + gen_vel = True + gen_dis = True + gen_acc = True extra_dname = 'nodeTag' # original unrotated node coordinates @@ -1639,32 +1647,32 @@ def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste user_y0[i] = df.loc[i, 'y'] user_z0[i] = df.loc[i, 'z'] - if mpi_rank == 0: - print('Generating motions for %i nodes...' % (n_coord)) - - output_fname = get_output_filename( - 'csv', output_format, save_path, csv_fname, explicit_output_mode=requested_output_mode is not None - ) - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) - + if mpi_rank == 0: + print('Generating motions for %i nodes...' % (n_coord)) + + output_fname = get_output_filename( + 'csv', output_format, save_path, csv_fname, explicit_output_mode=requested_output_mode is not None + ) + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) + return def dframeToDict(dFrame): dFrame = list(dFrame.iterrows()) return {i[1].to_list()[0] : i[1].to_list()[1] for i in dFrame} -def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, ref_coord=None, requested_output_mode=None): - if mpi_rank == 0: - print('Input CSV [%s]' %csv_fname) - print('Input ESSI [%s]' %ssi_fname) - - if ref_coord is None: - ref_coord = np.zeros(3) - - output_format = resolve_output_mode('template', requested_output_mode) - - sw4ToESSI_params = dframeToDict(pd.read_csv(csv_fname)) +def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, ref_coord=None, requested_output_mode=None): + if mpi_rank == 0: + print('Input CSV [%s]' %csv_fname) + print('Input ESSI [%s]' %ssi_fname) + + if ref_coord is None: + ref_coord = np.zeros(3) + + output_format = resolve_output_mode('template', requested_output_mode) + + sw4ToESSI_params = dframeToDict(pd.read_csv(csv_fname)) sw4_i_start = sw4ToESSI_params["sw4_i_start"] sw4_i_end = sw4ToESSI_params["sw4_i_end"] sw4_j_start = sw4ToESSI_params["sw4_j_start"] @@ -1681,29 +1689,29 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep # reference point, which is the ESSI or OpenSees origin in the SW4 coordinate system ref_coord[0] = essi_x_start ref_coord[1] = essi_y_start - ref_coord[2] = essi_z_start + ref_coord[2] = essi_z_start - coord_sys = ['y', 'x', '-z'] - if output_format == 'opensees': - gen_vel = False - gen_dis = True - gen_acc = True - extra_dname = 'internal' - else: - gen_vel = True - gen_dis = True - gen_acc = True - extra_dname = 'Is Boundary Node' - - # original unrotated node coordinates - template_file = h5py.File(template_fname) - coordinates = template_file['Coordinates'][:] - - is_boundary = np.asarray(template_file['Is Boundary Node'][:]).reshape(-1) - n_coord = len(is_boundary) - user_x = np.zeros(n_coord) - user_y = np.zeros(n_coord) - user_z = np.zeros(n_coord) + coord_sys = ['y', 'x', '-z'] + if output_format == 'opensees': + gen_vel = False + gen_dis = True + gen_acc = True + extra_dname = 'internal' + else: + gen_vel = True + gen_dis = True + gen_acc = True + extra_dname = 'Is Boundary Node' + + # original unrotated node coordinates + template_file = h5py.File(template_fname) + coordinates = template_file['Coordinates'][:] + + is_boundary = np.asarray(template_file['Is Boundary Node'][:]).reshape(-1) + n_coord = len(is_boundary) + user_x = np.zeros(n_coord) + user_y = np.zeros(n_coord) + user_z = np.zeros(n_coord) for i in range(0, n_coord): user_x[i] = coordinates[i*3] user_y[i] = coordinates[i*3+1] @@ -1711,21 +1719,21 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep template_file.close() - if verbose and mpi_rank == 0: - print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x[0], user_y[0], user_z[0], user_x[-1], user_y[-1], user_z[-1])) - print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x), np.max(user_x), np.min(user_y), np.max(user_y), np.min(user_z), np.max(user_z)) ) - print('Start/end time', start_t, end_t) - - output_fname = get_output_filename( - 'template', output_format, save_path, template_fname, explicit_output_mode=requested_output_mode is not None - ) - if output_format == 'essi': - if mpi_rank == 0: - prepare_essi_output_file(template_fname, output_fname) - MPI.COMM_WORLD.Barrier() - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x, user_y, user_z, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, is_boundary, extra_dname, output_format) - return + if verbose and mpi_rank == 0: + print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x[0], user_y[0], user_z[0], user_x[-1], user_y[-1], user_z[-1])) + print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x), np.max(user_x), np.min(user_y), np.max(user_y), np.min(user_z), np.max(user_z)) ) + print('Start/end time', start_t, end_t) + + output_fname = get_output_filename( + 'template', output_format, save_path, template_fname, explicit_output_mode=requested_output_mode is not None + ) + if output_format == 'essi': + if mpi_rank == 0: + prepare_essi_output_file(template_fname, output_fname) + MPI.COMM_WORLD.Barrier() + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x, user_y, user_z, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, is_boundary, extra_dname, output_format) + return if __name__ == "__main__": @@ -1747,12 +1755,12 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep start_t=0 end_t=0 tstep = 1 - rotate_angle = 0 - zeroMotionDir = 'None' - requested_output_mode = None + rotate_angle = 0 + zeroMotionDir = 'None' + requested_output_mode = None - parser = build_arg_parser() - args = parser.parse_args() + parser = build_arg_parser() + args = parser.parse_args() if args.verbose: verbose=True @@ -1775,8 +1783,8 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep if args.template: template_fname=args.template use_template=True - if args.ssi: - ssi_fname=args.ssi + if args.ssi: + ssi_fname=args.ssi if args.reference: ref_coord[0]=args.reference[0] ref_coord[1]=args.reference[1] @@ -1787,12 +1795,12 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep tstep = args.timerange[2] # time step interval if args.rotateanlge: rotate_angle = args.rotateanlge - if args.savepath: - save_path = args.savepath - if args.zeroMotionDir: - zeroMotionDir = args.zeroMotionDir - if args.output_format: - requested_output_mode = args.output_format + if args.savepath: + save_path = args.savepath + if args.zeroMotionDir: + zeroMotionDir = args.zeroMotionDir + if args.output_format: + requested_output_mode = args.output_format comm = MPI.COMM_WORLD mpi_size = comm.Get_size() @@ -1806,31 +1814,31 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep startTimeStr = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime(startTime)) print('Start time:', startTimeStr) - if drm_fname == '' and h5_fname == '' and csv_fname == '' and template_fname == '': - print('Error, no node coordinate input file is provided, exit...') - exit(0) + if drm_fname == '' and h5_fname == '' and csv_fname == '' and template_fname == '': + print('Error, no node coordinate input file is provided, exit...') + exit(0) if ssi_fname == '': - print('Error, no SW4 SSI output file is provided, exit...') + print('Error, no SW4 SSI output file is provided, exit...') exit(0) - if verbose and mpi_rank == 0: - print('Using ref_coord={}, start_t={}, end_t={}, tstep={}, rotate_angle={}, output_mode={} to extract motions'.format(ref_coord, start_t, end_t, tstep, rotate_angle, requested_output_mode or 'default')) - - try: - if use_drm: - convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_h5: - convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_csv and not use_template: - convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_csv and use_template: - convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, ref_coord, requested_output_mode) - except ValueError as exc: - if mpi_rank == 0: - print(f'Error: {exc}') - raise SystemExit(1) + if verbose and mpi_rank == 0: + print('Using ref_coord={}, start_t={}, end_t={}, tstep={}, rotate_angle={}, output_mode={} to extract motions'.format(ref_coord, start_t, end_t, tstep, rotate_angle, requested_output_mode or 'default')) + + try: + if use_drm: + convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_h5: + convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_csv and not use_template: + convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_csv and use_template: + convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, ref_coord, requested_output_mode) + except ValueError as exc: + if mpi_rank == 0: + print(f'Error: {exc}') + raise SystemExit(1) if mpi_rank == 0: endTime = time.time() endTimeStr = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime(endTime)) - print(f'End time: {endTimeStr} (time spent: {endTime-startTime:.2f} s)') + print(f'End time: {endTimeStr} (time spent: {endTime-startTime:.2f} s)') diff --git a/template/motion_setting.csv b/template/motion_setting.csv index cf3078c..e807fc6 100644 --- a/template/motion_setting.csv +++ b/template/motion_setting.csv @@ -1,5 +1,5 @@ nodeTag,x,y,z,essiXstart,essiYstart,essiZstart,startTime,endTime,tstep,rotationAngle,zeroMotionDir -1,0,0,0,12000,40000,24,0,0,5,0,None +1,0,0,0,12000,40000,24,0,0,5,0,None 2,8,0,0,,,,,,,, 3,8,8,0,,,,,,,, 4,0,8,0,,,,,,,, From b29152edfe5395ce4c112db314a3a2251d1f2641 Mon Sep 17 00:00:00 2001 From: jfhuang Date: Wed, 29 Jan 2025 22:16:55 -0800 Subject: [PATCH 2/6] Sync to include recent updates --- convert.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/convert.py b/convert.py index 980323c..64292d2 100755 --- a/convert.py +++ b/convert.py @@ -820,6 +820,11 @@ def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) + + str_crd = coord_to_str_3d(x0, y0, z0) + if str_crd == '2900, 7100, 1': + print(f'{str_crd}: cid: {val}') + return val def chunkid_to_start(cid, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z): From 51762380464b2be251bb9be4607805f800dd9625 Mon Sep 17 00:00:00 2001 From: Junfei Huang Date: Mon, 20 Apr 2026 15:19:18 -0700 Subject: [PATCH 3/6] Switch on opensees & essi output format for csv & h5 input kinds --- convert.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/convert.py b/convert.py index 64292d2..62326e2 100755 --- a/convert.py +++ b/convert.py @@ -75,8 +75,8 @@ def resolve_output_mode(input_kind, requested_output_mode): } supported = { 'drm': {'point', 'opensees', 'essi'}, - 'h5': {'point'}, - 'csv': {'point'}, + 'h5': {'point', 'opensees', 'essi'}, + 'csv': {'point', 'opensees', 'essi'}, 'template': {'point', 'opensees', 'essi'}, } @@ -1037,6 +1037,8 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # rotate/transform only when rotate_angle is other than 0 (default min difference is 1e-2) b_rotate = np.where(abs(rotate_angle) > 1e-2, True, False) if b_rotate: + if verbose and mpi_rank == 0: + print(f'Rotate coordinates by {rotate_angle} degree around vertical axis (z) in the xy plane') user_x, user_y, user_z = rotate_coords_ops_xyplane(user_x, user_y, user_z, rotate_angle) # Convert user coordinate to sw4 coordinate, relative to ESSI domain (subset of SW4 domain) @@ -1053,7 +1055,8 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # tmpfile.write('(%d, %d, %d)' % (user_essi_x[i], user_essi_y[i], user_essi_z[i])) # print('user_essi_x.dtype: ', user_essi_x.dtype) # print('user_x0.dtype: ', user_x0.dtype) - # np.savetxt('tmp.txt', np.c_[user_essi_x, user_essi_y, user_essi_z], fmt='%.6f') + # np.savetxt('tmp_user.txt', np.c_[user_x, user_y, user_z], fmt='%.6f') + # np.savetxt('tmp_essi.txt', np.c_[user_essi_x, user_essi_y, user_essi_z], fmt='%.6f') # Plot if mpi_rank == 0: From 988c39b1c9c1bfe79470726d5b699d394b1aa9da Mon Sep 17 00:00:00 2001 From: Houjun Tang Date: Tue, 21 Apr 2026 11:58:44 -0700 Subject: [PATCH 4/6] Fix CI for expanded output format support --- README.md | 4 +- convert.py | 42 +++++--- tests/test_convert_fixture_data.py | 156 +++++++++++++++++++++++++---- tests/test_convert_helpers.py | 14 ++- 4 files changed, 179 insertions(+), 37 deletions(-) diff --git a/README.md b/README.md index 1c928ab..aef1ad1 100644 --- a/README.md +++ b/README.md @@ -53,8 +53,8 @@ When `--output-format` is omitted, the script keeps the current default for each Notes: - `--timerange start end step` is in SSI time units. - If `start == end`, the script keeps all steps from `start` to the end of the SSI record. -- `--output-format point` is supported for DRM, HDF5, CSV, and template-driven inputs. -- `--output-format opensees` and `--output-format essi` require input that carries boundary-node metadata, so they are only supported for DRM or template-driven inputs. +- `--output-format point`, `--output-format opensees`, and `--output-format essi` are supported for DRM, HDF5, CSV, and template-driven inputs. +- HDF5 and CSV inputs keep their `nodeTag` metadata when writing OpenSees or ESSI-style outputs. ## Examples diff --git a/convert.py b/convert.py index 62326e2..50c8203 100755 --- a/convert.py +++ b/convert.py @@ -109,10 +109,13 @@ def get_output_filename(input_kind, output_mode, save_path, source_fname, explic output_stem = os.path.splitext(os.path.basename(source_fname))[0] return os.path.join(save_path, f'{output_stem}_motion.h5') - if output_mode == 'essi': - if input_kind == 'template' and not explicit_output_mode: - return source_fname - return os.path.join(save_path, os.path.basename(source_fname)) + if output_mode == 'essi': + if input_kind == 'template' and not explicit_output_mode: + return source_fname + output_name = os.path.basename(source_fname) + if input_kind == 'csv': + output_name = f'{os.path.splitext(output_name)[0]}.h5' + return os.path.join(save_path, output_name) raise ValueError(f'Unsupported output format: {output_mode}') @@ -793,12 +796,25 @@ def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen h5file.close() -def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, extra_dname): - h5file = h5py.File(h5_fname, 'r+') - - for dname in ('Velocity', 'Accelerations', 'Displacements', 'Time'): - if dname in h5file: - del h5file[dname] +def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, extra_dname): + h5file = h5py.File(h5_fname, 'a') + + coord_shape = (ncoord * 3,) + meta_shape = (ncoord,) + + if 'Coordinates' in h5file and h5file['Coordinates'].shape != coord_shape: + del h5file['Coordinates'] + if 'Coordinates' not in h5file: + h5file.create_dataset('Coordinates', coord_shape, dtype='f8') + + if extra_dname in h5file and h5file[extra_dname].shape != meta_shape: + del h5file[extra_dname] + if extra_dname not in h5file: + h5file.create_dataset(extra_dname, meta_shape, dtype='i4') + + for dname in ('Velocity', 'Accelerations', 'Displacements', 'Time'): + if dname in h5file: + del h5file[dname] if gen_vel: dset = h5file.create_dataset('Velocity', (ncoord*3, nstep), dtype='f4') @@ -1647,9 +1663,9 @@ def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste df = pd.read_csv(csv_fname) node_tags = df['nodeTag'][:].tolist() n_coord = len(node_tags) - user_x0 = np.zeros(n_coord, dtype='f4') - user_y0 = np.zeros(n_coord, dtype='f4') - user_z0 = np.zeros(n_coord, dtype='f4') + user_x0 = np.zeros(n_coord, dtype='f8') + user_y0 = np.zeros(n_coord, dtype='f8') + user_z0 = np.zeros(n_coord, dtype='f8') for i in range(0, n_coord): user_x0[i] = df.loc[i, 'x'] user_y0[i] = df.loc[i, 'y'] diff --git a/tests/test_convert_fixture_data.py b/tests/test_convert_fixture_data.py index 547688c..9a6c1fa 100644 --- a/tests/test_convert_fixture_data.py +++ b/tests/test_convert_fixture_data.py @@ -1,3 +1,4 @@ +import csv import os import sys import tempfile @@ -48,6 +49,45 @@ def setUp(self): self.sample_csv = REPO_ROOT / "template" / "motion_setting.csv" self.reference_output = REPO_ROOT / "tests" / "data" / "h5NodeMotion.h5" + def _write_h5_based_csv_fixture(self, output_path): + ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta( + str(self.sample_csv) + ) + + with h5py.File(self.sample_h5, "r") as input_h5: + coordinates = input_h5["coordinate"][:] + node_tags = input_h5["nodeTag"][:] + + if coordinates.ndim == 1: + coordinates = coordinates.reshape(-1, 3) + elif coordinates.ndim == 2 and coordinates.shape[1] == 1: + coordinates = coordinates.reshape(-1, 3) + + with output_path.open("w", newline="") as csvfile: + writer = csv.writer(csvfile) + writer.writerow( + [ + "nodeTag", + "x", + "y", + "z", + "essiXstart", + "essiYstart", + "essiZstart", + "startTime", + "endTime", + "tstep", + "rotationAngle", + "zeroMotionDir", + ] + ) + for index, (node_tag, coord) in enumerate(zip(node_tags, coordinates)): + row = [int(node_tag), float(coord[0]), float(coord[1]), float(coord[2]), "", "", "", "", "", "", "", ""] + if index == 0: + row[4:12] = [ref_coord[0], ref_coord[1], ref_coord[2], start_t, end_t, tstep, rotate_angle, zero_motion_dir] + writer.writerow(row) + + def test_convert_h5_matches_checked_in_fixture_data(self): with tempfile.TemporaryDirectory() as tmpdir: ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta( @@ -116,30 +156,112 @@ def test_convert_h5_accepts_explicit_point_output_mode(self): np.testing.assert_allclose(output_h5["xyz"][:], reference_h5["xyz"][:]) np.testing.assert_array_equal(output_h5["nodeTag"][:], reference_h5["nodeTag"][:]) - def test_convert_h5_rejects_unsupported_essi_output_mode(self): + def test_convert_h5_accepts_explicit_essi_output_mode(self): with tempfile.TemporaryDirectory() as tmpdir: ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta( str(self.sample_csv) ) - with self.assertRaisesRegex(ValueError, 'not supported for h5 input'): - convert.convert_h5( - str(self.sample_h5), - str(self.sample_ssi), - tmpdir, - ref_coord, - start_t, - end_t, - tstep, - rotate_angle, - zero_motion_dir, - False, - 0, - 1, - False, - requested_output_mode="essi", + convert.convert_h5( + str(self.sample_h5), + str(self.sample_ssi), + tmpdir, + ref_coord, + start_t, + end_t, + tstep, + rotate_angle, + zero_motion_dir, + False, + 0, + 1, + False, + requested_output_mode="essi", + ) + + output_path = Path(tmpdir) / "h5NodeCrds.h5" + self.assertTrue(output_path.exists(), f"Missing output file: {output_path}") + + with h5py.File(output_path, "r") as output_h5, h5py.File(self.reference_output, "r") as reference_h5: + self.assertEqual(output_h5["Coordinates"].shape, (108,)) + self.assertEqual(output_h5["Velocity"].shape, (108, 600)) + self.assertEqual(output_h5["Accelerations"].shape, (108, 600)) + self.assertEqual(output_h5["Displacements"].shape, (108, 600)) + self.assertEqual(output_h5["Time"].shape, (600,)) + np.testing.assert_array_equal(output_h5["nodeTag"][:], reference_h5["nodeTag"][:]) + timestep = float(output_h5["Time"][1] - output_h5["Time"][0]) + np.testing.assert_allclose(output_h5["Time"][:], np.arange(600, dtype="f8") * timestep) + + def test_convert_h5_accepts_explicit_opensees_output_mode(self): + with tempfile.TemporaryDirectory() as tmpdir: + ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta( + str(self.sample_csv) + ) + + convert.convert_h5( + str(self.sample_h5), + str(self.sample_ssi), + tmpdir, + ref_coord, + start_t, + end_t, + tstep, + rotate_angle, + zero_motion_dir, + False, + 0, + 1, + False, + requested_output_mode="opensees", + ) + + output_path = Path(tmpdir) / "OpenSeesDRMinput.h5drm" + self.assertTrue(output_path.exists(), f"Missing output file: {output_path}") + + with h5py.File(output_path, "r") as output_h5, h5py.File(self.reference_output, "r") as reference_h5: + self.assertEqual(output_h5["DRM_Data"]["acceleration"].shape, (108, 600)) + self.assertEqual(output_h5["DRM_Data"]["displacement"].shape, (108, 600)) + self.assertEqual(output_h5["DRM_Data"]["nodeTag"].shape, (36,)) + self.assertEqual(output_h5["DRM_Data"]["xyz"].shape, (36, 3)) + self.assertEqual( + float(output_h5["DRM_Metadata"]["dt"][()]), + float(reference_h5["dt"][()]), ) + def test_convert_csv_accepts_explicit_essi_output_mode(self): + with tempfile.TemporaryDirectory() as tmpdir: + csv_path = Path(tmpdir) / "input.csv" + self._write_h5_based_csv_fixture(csv_path) + ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta(str(csv_path)) + + convert.convert_csv( + str(csv_path), + str(self.sample_ssi), + tmpdir, + ref_coord, + start_t, + end_t, + tstep, + rotate_angle, + zero_motion_dir, + False, + 0, + 1, + False, + requested_output_mode="essi", + ) + + output_path = Path(tmpdir) / "input.csv" + self.assertTrue(output_path.exists(), f"Missing source CSV fixture: {output_path}") + essi_output = Path(tmpdir) / "input.h5" + self.assertTrue(essi_output.exists(), f"Missing output file: {essi_output}") + + with h5py.File(essi_output, "r") as output_h5: + self.assertEqual(output_h5["Coordinates"].shape, (108,)) + self.assertEqual(output_h5["nodeTag"].shape, (36,)) + self.assertEqual(output_h5["Velocity"].shape, (108, 600)) + self.assertEqual(output_h5["Time"].shape, (600,)) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_convert_helpers.py b/tests/test_convert_helpers.py index ba5c56c..cef61c8 100644 --- a/tests/test_convert_helpers.py +++ b/tests/test_convert_helpers.py @@ -26,11 +26,11 @@ def test_resolve_output_mode_keeps_existing_defaults(self): self.assertEqual(convert.resolve_output_mode("csv", None), "point") self.assertEqual(convert.resolve_output_mode("template", None), "essi") - def test_resolve_output_mode_rejects_unsupported_input_output_pairs(self): - with self.assertRaisesRegex(ValueError, 'not supported for h5 input'): - convert.resolve_output_mode("h5", "essi") - with self.assertRaisesRegex(ValueError, 'not supported for csv input'): - convert.resolve_output_mode("csv", "opensees") + def test_resolve_output_mode_accepts_explicit_overrides_for_all_input_kinds(self): + self.assertEqual(convert.resolve_output_mode("h5", "essi"), "essi") + self.assertEqual(convert.resolve_output_mode("csv", "opensees"), "opensees") + with self.assertRaisesRegex(ValueError, "Supported formats"): + convert.resolve_output_mode("h5", "unsupported") def test_get_flat_coord_range_scales_node_offsets_by_xyz_width(self): self.assertEqual(convert.get_flat_coord_range(0, 2), (0, 6)) @@ -49,6 +49,10 @@ def test_get_output_filename_uses_requested_mode_and_input_kind(self): convert.get_output_filename("template", "essi", "/tmp/out", "/tmp/in/template.h5drm", explicit_output_mode=True), "/tmp/out/template.h5drm", ) + self.assertEqual( + convert.get_output_filename("csv", "essi", "/tmp/out", "/tmp/in/input.csv", explicit_output_mode=True), + "/tmp/out/input.h5", + ) def test_create_hdf5_essi_writes_float_time_axis_matching_motion_length(self): with tempfile.TemporaryDirectory() as tmpdir: From 99bacc0c36e008d1a426e074c8a0bae5cc95809b Mon Sep 17 00:00:00 2001 From: Houjun Tang Date: Tue, 21 Apr 2026 12:09:14 -0700 Subject: [PATCH 5/6] Address PR review feedback --- convert.py | 65 ++++++++++++++---------------- tests/test_convert_fixture_data.py | 38 +++++++++++++++++ 2 files changed, 68 insertions(+), 35 deletions(-) diff --git a/convert.py b/convert.py index 50c8203..b156d3f 100755 --- a/convert.py +++ b/convert.py @@ -829,19 +829,14 @@ def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, ext h5file.close() -def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx, ny, nz): +def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx, ny, nz): # Due to allocation of neighbor nodes, the indices may exceed the original dataset # boundary and thus need to be corrected so that the correct nearest chunk can be fetched later x, y, z = min(x0, nx-1), min(y0, ny-1), min(z0, nz-1) - val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) - #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) - - str_crd = coord_to_str_3d(x0, y0, z0) - if str_crd == '2900, 7100, 1': - print(f'{str_crd}: cid: {val}') - - return val + val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) + #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) + return val def chunkid_to_start(cid, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z): #print('cid2:', cid, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z) @@ -976,7 +971,7 @@ def linear_interp(data_dict, x, y, z): result = ((c000 * (1-xd) + c100 * xd) * (1-yd) + (c010 * (1-xd) + c110 * xd) * yd) * (1-zd) + ((c001 * (1-xd) + c101 * xd) * (1-yd) + (c011 * (1-xd) + c111 * xd) * yd) * zd return result -def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_data, extra_dname, output_format): +def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values_input, extra_dname, output_format): if zeroMotionDir is None or pd.isna(zeroMotionDir): zeroMotionDir = 'None' elif isinstance(zeroMotionDir, str): @@ -1198,7 +1193,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use read_coords_vel_1 = {} read_coords_vel_2 = {} # coords_str_dict = {} - is_boundary = np.zeros(n_coord, dtype='i4') + extra_values = np.zeros(n_coord, dtype='i4') my_ncoord = np.zeros(1, dtype='int') my_user_coordinates = np.zeros((n_coord,3), dtype='f8') my_converted_coordinates = np.zeros((n_coord,3), dtype='f8') @@ -1256,7 +1251,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use else: my_cids_dict[cid] = {coord_str} - is_boundary[my_ncoord[0]] = extra_data[i] + extra_values[my_ncoord[0]] = extra_values_input[i] my_ncoord[0] += 1 #end if assigned to my rank #end for i in all coordinates @@ -1267,7 +1262,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Allocated more than needed previously, resize here my_user_coordinates.resize(my_ncoord[0], 3) my_converted_coordinates.resize(my_ncoord[0], 3) - is_boundary.resize(my_ncoord[0]) + extra_values.resize(my_ncoord[0]) # if mpi_rank == 0: # # print('read_coords_vel_0 =', read_coords_vel_0) @@ -1484,7 +1479,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if my_ncoord[0] > 0: write_to_hdf5_range_2d(output_fname, 'DRM_Data', 'xyz', my_user_coordinates, my_offset, (my_offset+my_ncoord[0])) - write_to_hdf5_range_1d(output_fname, 'DRM_Data', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, 'DRM_Data', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, 'DRM_Data', 'acceleration', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: @@ -1499,7 +1494,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if verbose: print('Rank', mpi_rank, 'start to write data') write_to_hdf5_range_2d(output_fname, 'DRM_Data', 'xyz', my_user_coordinates, my_offset, (my_offset+my_ncoord[0])) - write_to_hdf5_range_1d(output_fname, 'DRM_Data', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, 'DRM_Data', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, 'DRM_Data', 'acceleration', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) @@ -1516,7 +1511,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if my_ncoord[0] > 0: write_to_hdf5_range_2d(output_fname, '/', 'xyz', my_user_coordinates, my_offset, (my_offset+my_ncoord[0])) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, '/', 'acceleration', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: @@ -1531,7 +1526,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if verbose: print('Rank', mpi_rank, 'start to write data') write_to_hdf5_range_2d(output_fname, '/', 'xyz', my_user_coordinates, my_offset, (my_offset+my_ncoord[0])) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, '/', 'acceleration', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) @@ -1549,7 +1544,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Write to the template file if my_ncoord[0] > 0: write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: @@ -1563,7 +1558,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use data = comm.recv(source=mpi_rank-1, tag=111) if my_ncoord[0] > 0: write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) - write_to_hdf5_range_1d(output_fname, '/', extra_dname, is_boundary, my_offset, my_offset+my_ncoord[0]) + write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) if gen_acc: write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: @@ -1591,8 +1586,8 @@ def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste output_format = resolve_output_mode('drm', requested_output_mode) # original unrotated node coordinates - user_x0, user_y0, user_z0, n_coord, isboundary = read_coord_drm(drm_fname, verbose) - isboundary = np.asarray(isboundary).reshape(-1) + user_x0, user_y0, user_z0, n_coord, boundary_flags = read_coord_drm(drm_fname, verbose) + extra_values = np.asarray(boundary_flags).reshape(-1) if verbose and mpi_rank == 0: print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x0[0], user_y0[0], user_z0[0], user_x0[-1], user_y0[-1], user_z0[-1])) print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x0), np.max(user_x0), np.min(user_y0), np.max(user_y0), np.min(user_z0), np.max(user_z0)) ) @@ -1616,7 +1611,7 @@ def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste prepare_essi_output_file(drm_fname, output_fname) MPI.COMM_WORLD.Barrier() - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, isboundary, extra_dname, output_format) + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values, extra_dname, output_format) return @@ -1625,11 +1620,11 @@ def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, print('Input h5 [%s]' %h5_fname) print('Input ESSI [%s]' %ssi_fname) - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('h5', requested_output_mode) - gen_vel = True - gen_dis = True - gen_acc = True + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('h5', requested_output_mode) + gen_vel = output_format != 'opensees' + gen_dis = True + gen_acc = True extra_dname = 'nodeTag' # original unrotated node coordinates @@ -1652,11 +1647,11 @@ def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste print('Input CSV [%s]' %csv_fname) print('Input ESSI [%s]' %ssi_fname) - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('csv', requested_output_mode) - gen_vel = True - gen_dis = True - gen_acc = True + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('csv', requested_output_mode) + gen_vel = output_format != 'opensees' + gen_dis = True + gen_acc = True extra_dname = 'nodeTag' # original unrotated node coordinates @@ -1731,8 +1726,8 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep template_file = h5py.File(template_fname) coordinates = template_file['Coordinates'][:] - is_boundary = np.asarray(template_file['Is Boundary Node'][:]).reshape(-1) - n_coord = len(is_boundary) + extra_values = np.asarray(template_file['Is Boundary Node'][:]).reshape(-1) + n_coord = len(extra_values) user_x = np.zeros(n_coord) user_y = np.zeros(n_coord) user_z = np.zeros(n_coord) @@ -1756,7 +1751,7 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep prepare_essi_output_file(template_fname, output_fname) MPI.COMM_WORLD.Barrier() - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x, user_y, user_z, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, is_boundary, extra_dname, output_format) + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x, user_y, user_z, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values, extra_dname, output_format) return if __name__ == "__main__": diff --git a/tests/test_convert_fixture_data.py b/tests/test_convert_fixture_data.py index 9a6c1fa..a37ea73 100644 --- a/tests/test_convert_fixture_data.py +++ b/tests/test_convert_fixture_data.py @@ -221,6 +221,7 @@ def test_convert_h5_accepts_explicit_opensees_output_mode(self): with h5py.File(output_path, "r") as output_h5, h5py.File(self.reference_output, "r") as reference_h5: self.assertEqual(output_h5["DRM_Data"]["acceleration"].shape, (108, 600)) self.assertEqual(output_h5["DRM_Data"]["displacement"].shape, (108, 600)) + self.assertNotIn("velocity", output_h5["DRM_Data"]) self.assertEqual(output_h5["DRM_Data"]["nodeTag"].shape, (36,)) self.assertEqual(output_h5["DRM_Data"]["xyz"].shape, (36, 3)) self.assertEqual( @@ -262,6 +263,43 @@ def test_convert_csv_accepts_explicit_essi_output_mode(self): self.assertEqual(output_h5["Velocity"].shape, (108, 600)) self.assertEqual(output_h5["Time"].shape, (600,)) + def test_convert_csv_accepts_explicit_opensees_output_mode(self): + with tempfile.TemporaryDirectory() as tmpdir: + csv_path = Path(tmpdir) / "input.csv" + self._write_h5_based_csv_fixture(csv_path) + ref_coord, start_t, end_t, tstep, rotate_angle, zero_motion_dir = convert.get_csv_meta(str(csv_path)) + + convert.convert_csv( + str(csv_path), + str(self.sample_ssi), + tmpdir, + ref_coord, + start_t, + end_t, + tstep, + rotate_angle, + zero_motion_dir, + False, + 0, + 1, + False, + requested_output_mode="opensees", + ) + + output_path = Path(tmpdir) / "OpenSeesDRMinput.h5drm" + self.assertTrue(output_path.exists(), f"Missing output file: {output_path}") + + with h5py.File(output_path, "r") as output_h5, h5py.File(self.reference_output, "r") as reference_h5: + self.assertEqual(output_h5["DRM_Data"]["acceleration"].shape, (108, 600)) + self.assertEqual(output_h5["DRM_Data"]["displacement"].shape, (108, 600)) + self.assertNotIn("velocity", output_h5["DRM_Data"]) + self.assertEqual(output_h5["DRM_Data"]["nodeTag"].shape, (36,)) + self.assertEqual(output_h5["DRM_Data"]["xyz"].shape, (36, 3)) + self.assertEqual( + float(output_h5["DRM_Metadata"]["dt"][()]), + float(reference_h5["dt"][()]), + ) + if __name__ == "__main__": unittest.main() From 76d9ffcc8d2663cd9fca22284fdaf074b7aee5cd Mon Sep 17 00:00:00 2001 From: Houjun Tang Date: Tue, 21 Apr 2026 12:14:40 -0700 Subject: [PATCH 6/6] Reduce formatting-only PR diff noise --- convert.py | 1191 +++++++++++++++++------------------ template/motion_setting.csv | 2 +- 2 files changed, 591 insertions(+), 602 deletions(-) diff --git a/convert.py b/convert.py index b156d3f..e8f5dd2 100755 --- a/convert.py +++ b/convert.py @@ -1,19 +1,19 @@ #!/usr/bin/env python3 # from genericpath import exists -import os -import shutil -# import sys -import argparse -import h5py -import math - -import scipy -try: - cumtrapz = scipy.integrate.cumulative_trapezoid -except AttributeError: - cumtrapz = scipy.integrate.cumtrapz -from scipy.interpolate import interpn +import os +import shutil +# import sys +import argparse +import h5py +import math + +import scipy +try: + cumtrapz = scipy.integrate.cumulative_trapezoid +except AttributeError: + cumtrapz = scipy.integrate.cumtrapz +from scipy.interpolate import interpn import numpy as np @@ -27,88 +27,88 @@ from mpl_toolkits.mplot3d import Axes3D from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection from mpi4py import MPI -import functools -print = functools.partial(print, flush=True) # Don't buffer print -import hdf5plugin # Use this when SW4 output uses ZFP compression, can be installed with "pip install hdf5plugin" - - -def build_arg_parser(): - parser = argparse.ArgumentParser() - parser.add_argument("-c", "--csv", help="full path to the CSV setting file", default="") - parser.add_argument("-d", "--drm", help="full path to the DRM file with node coordinates", default="") - parser.add_argument("-h5", "--hdf5", help="full path to the hdf5 file with node coordinates", default="") - parser.add_argument("-e", "--ssi", "--essi", dest="ssi", - help="full path to the SW4 SSI output file", default="") - parser.add_argument("-t", "--template", - help="full path to the SSI template file with node coordinates", default="") - parser.add_argument("-p", "--plotonly", help="only generate plots of the input nodes", action="store_true") - parser.add_argument("-r", "--reference", help="reference node coordinate offset, default 0 0 0", - nargs='+', type=float) - # parser.add_argument("-s", "--steprange", help="timestep range, default 0 total_steps", nargs='+', type=int) - parser.add_argument("-s", "--timerange", - help="time range, will return all steps after the lower limit for equal upper and lower limit", - nargs='+', type=int) - parser.add_argument("-R", "--rotateanlge", help="rotate angle for node coordinate and motion: [0, 360)", - type=float) - parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") - parser.add_argument("-P", "--savepath", help="full path for saving the result files", default="") - parser.add_argument( - "-o", - "--output-format", - dest="output_format", - choices=("point", "opensees", "essi"), - help="output format override: point motions, OpenSees DRM, or ESSI template output", - default=None, - ) - parser.add_argument("-z", "--zeroMotionDir", - help="direction for zeroing out motion and enforce same motion across nodes in that direction: None(default), x, y, z", - default="") - return parser - - -def resolve_output_mode(input_kind, requested_output_mode): - defaults = { - 'drm': 'opensees', - 'h5': 'point', - 'csv': 'point', - 'template': 'essi', - } - supported = { - 'drm': {'point', 'opensees', 'essi'}, - 'h5': {'point', 'opensees', 'essi'}, - 'csv': {'point', 'opensees', 'essi'}, - 'template': {'point', 'opensees', 'essi'}, - } - - if requested_output_mode is None: - return defaults[input_kind] - - output_mode = requested_output_mode.lower() - if output_mode not in supported[input_kind]: - supported_modes = ', '.join(sorted(supported[input_kind])) - raise ValueError( - f'Output format "{output_mode}" is not supported for {input_kind} input. ' - f'Supported formats: {supported_modes}.' - ) - return output_mode - - -def get_output_filename(input_kind, output_mode, save_path, source_fname, explicit_output_mode=False): - if output_mode == 'opensees': - return os.path.join(save_path, 'OpenSeesDRMinput.h5drm') - - if output_mode == 'point': - if input_kind == 'h5': - output_stem = os.path.splitext(os.path.basename(source_fname))[0] - return os.path.join(save_path, f'{output_stem}_motion.h5') - if input_kind == 'csv': - return os.path.join(save_path, 'csvNodeMotion.h5') - if input_kind == 'drm': - return os.path.join(save_path, 'drmNodeMotion.h5') - if input_kind == 'template': - output_stem = os.path.splitext(os.path.basename(source_fname))[0] - return os.path.join(save_path, f'{output_stem}_motion.h5') - +import functools +print = functools.partial(print, flush=True) # Don't buffer print +import hdf5plugin # Use this when SW4 output uses ZFP compression, can be installed with "pip install hdf5plugin" + + +def build_arg_parser(): + parser = argparse.ArgumentParser() + parser.add_argument("-c", "--csv", help="full path to the CSV setting file", default="") + parser.add_argument("-d", "--drm", help="full path to the DRM file with node coordinates", default="") + parser.add_argument("-h5", "--hdf5", help="full path to the hdf5 file with node coordinates", default="") + parser.add_argument("-e", "--ssi", "--essi", dest="ssi", + help="full path to the SW4 SSI output file", default="") + parser.add_argument("-t", "--template", + help="full path to the SSI template file with node coordinates", default="") + parser.add_argument("-p", "--plotonly", help="only generate plots of the input nodes", action="store_true") + parser.add_argument("-r", "--reference", help="reference node coordinate offset, default 0 0 0", + nargs='+', type=float) + # parser.add_argument("-s", "--steprange", help="timestep range, default 0 total_steps", nargs='+', type=int) + parser.add_argument("-s", "--timerange", + help="time range, will return all steps after the lower limit for equal upper and lower limit", + nargs='+', type=int) + parser.add_argument("-R", "--rotateanlge", help="rotate angle for node coordinate and motion: [0, 360)", + type=float) + parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") + parser.add_argument("-P", "--savepath", help="full path for saving the result files", default="") + parser.add_argument( + "-o", + "--output-format", + dest="output_format", + choices=("point", "opensees", "essi"), + help="output format override: point motions, OpenSees DRM, or ESSI template output", + default=None, + ) + parser.add_argument("-z", "--zeroMotionDir", + help="direction for zeroing out motion and enforce same motion across nodes in that direction: None(default), x, y, z", + default="") + return parser + + +def resolve_output_mode(input_kind, requested_output_mode): + defaults = { + 'drm': 'opensees', + 'h5': 'point', + 'csv': 'point', + 'template': 'essi', + } + supported = { + 'drm': {'point', 'opensees', 'essi'}, + 'h5': {'point', 'opensees', 'essi'}, + 'csv': {'point', 'opensees', 'essi'}, + 'template': {'point', 'opensees', 'essi'}, + } + + if requested_output_mode is None: + return defaults[input_kind] + + output_mode = requested_output_mode.lower() + if output_mode not in supported[input_kind]: + supported_modes = ', '.join(sorted(supported[input_kind])) + raise ValueError( + f'Output format "{output_mode}" is not supported for {input_kind} input. ' + f'Supported formats: {supported_modes}.' + ) + return output_mode + + +def get_output_filename(input_kind, output_mode, save_path, source_fname, explicit_output_mode=False): + if output_mode == 'opensees': + return os.path.join(save_path, 'OpenSeesDRMinput.h5drm') + + if output_mode == 'point': + if input_kind == 'h5': + output_stem = os.path.splitext(os.path.basename(source_fname))[0] + return os.path.join(save_path, f'{output_stem}_motion.h5') + if input_kind == 'csv': + return os.path.join(save_path, 'csvNodeMotion.h5') + if input_kind == 'drm': + return os.path.join(save_path, 'drmNodeMotion.h5') + if input_kind == 'template': + output_stem = os.path.splitext(os.path.basename(source_fname))[0] + return os.path.join(save_path, f'{output_stem}_motion.h5') + if output_mode == 'essi': if input_kind == 'template' and not explicit_output_mode: return source_fname @@ -116,25 +116,17 @@ def get_output_filename(input_kind, output_mode, save_path, source_fname, explic if input_kind == 'csv': output_name = f'{os.path.splitext(output_name)[0]}.h5' return os.path.join(save_path, output_name) - - raise ValueError(f'Unsupported output format: {output_mode}') - - -def prepare_essi_output_file(source_fname, output_fname): - if os.path.abspath(source_fname) != os.path.abspath(output_fname): - shutil.copyfile(source_fname, output_fname) - -# from scipy.signal import butter,filtfilt -# def butter_lowpass_filter(data, cutoff, nyq, order): -# normal_cutoff = cutoff / nyq -# # Get the filter coefficients -# b, a = butter(order, normal_cutoff, btype='low', analog=False) -# y = filtfilt(b, a, data, axis=0, method="gust") -# return y - -# from scipy.signal import butter,filtfilt -# def butter_lowpass_filter(data, cutoff, nyq, order): -# normal_cutoff = cutoff / nyq + + raise ValueError(f'Unsupported output format: {output_mode}') + + +def prepare_essi_output_file(source_fname, output_fname): + if os.path.abspath(source_fname) != os.path.abspath(output_fname): + shutil.copyfile(source_fname, output_fname) + +# from scipy.signal import butter,filtfilt +# def butter_lowpass_filter(data, cutoff, nyq, order): +# normal_cutoff = cutoff / nyq # # Get the filter coefficients # b, a = butter(order, normal_cutoff, btype='low', analog=False) # y = filtfilt(b, a, data, axis=0, method="gust") @@ -229,23 +221,23 @@ def plot_cube(save_path, cube_definition, x, y, z, view): ax.yaxis.set_tick_params(labelsize=lblsize) ax.xaxis.set_tick_params(labelsize=lblsize) - # ax.dist = 12 - #ax.set_aspect('equal') - if hasattr(ax, 'set_box_aspect'): - box_aspect = ( - max(x_plot_max - x_plot_min, 1e-12), - max(y_plot_max - y_plot_min, 1e-12), - max(z_plot_max - z_plot_min, 1e-12), - ) - try: - ax.set_box_aspect(box_aspect, zoom=0.5) - except TypeError: - try: - ax.set_box_aspect(box_aspect) - except TypeError: - ax.dist = 12 - else: - ax.dist = 12 + # ax.dist = 12 + #ax.set_aspect('equal') + if hasattr(ax, 'set_box_aspect'): + box_aspect = ( + max(x_plot_max - x_plot_min, 1e-12), + max(y_plot_max - y_plot_min, 1e-12), + max(z_plot_max - z_plot_min, 1e-12), + ) + try: + ax.set_box_aspect(box_aspect, zoom=0.5) + except TypeError: + try: + ax.set_box_aspect(box_aspect) + except TypeError: + ax.dist = 12 + else: + ax.dist = 12 ax.text(cube_definition[2][0], cube_definition[2][1], cube_definition[2][2]-z_len*.05, 'SW4-ESSI domain', fontsize=7) @@ -301,95 +293,95 @@ def plot_coords(essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, us plot_cube(save_path, cube_definition, user_essi_x, user_essi_y, user_essi_z, 'XZ') plot_cube(save_path, cube_definition, user_essi_x, user_essi_y, user_essi_z, 'XY') -def read_coord_drm(drm_filename, verbose): - if verbose: - print('Reading coordinates from input file [%s]' % drm_filename) - - # Get the coordinates from DRM file - with h5py.File(drm_filename, 'r') as drm_file: - coordinates = drm_file['Coordinates'] - isboundary = drm_file['Is Boundary Node'][:] - - if coordinates.ndim == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('Coordinates dataset length must be divisible by 3 for 1D input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('Coordinates dataset length must be divisible by 3 for single-column input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 3: - n_coord = coordinates.shape[0] - else: - raise ValueError('Coordinates dataset must have shape (3*n,), (3*n, 1), or (n, 3)') - - drm_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on - drm_y = np.zeros(n_coord, dtype='f8') - drm_z = np.zeros(n_coord, dtype='f8') - - if coordinates.ndim == 1: - for i in range(0, n_coord): - drm_x[i] = coordinates[i*3] - drm_y[i] = coordinates[i*3+1] - drm_z[i] = coordinates[i*3+2] - elif coordinates.shape[1] == 1: - for i in range(0, n_coord): - drm_x[i] = coordinates[i*3, 0] - drm_y[i] = coordinates[i*3+1, 0] - drm_z[i] = coordinates[i*3+2, 0] - else: # coordinates.shape[1] == 3 - drm_x[:] = coordinates[:,0] - drm_y[:] = coordinates[:,1] - drm_z[:] = coordinates[:,2] - - return drm_x, drm_y, drm_z, n_coord, isboundary - -def read_coord_h5(h5_filename, verbose): - if verbose: - print('Reading coordinates from input file [%s]' % h5_filename) - - # Get the coordinates from h5 file - with h5py.File(h5_filename, 'r') as h5_file: - coordinates = h5_file['coordinate'] - if coordinates.ndim == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('coordinate dataset length must be divisible by 3 for 1D input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 1: - if coordinates.shape[0] % 3 != 0: - raise ValueError('coordinate dataset length must be divisible by 3 for single-column input') - n_coord = coordinates.shape[0] // 3 - elif coordinates.shape[1] == 3: - n_coord = coordinates.shape[0] - else: - raise ValueError('coordinate dataset must have shape (3*n,), (3*n, 1), or (n, 3)') - - h5_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on - h5_y = np.zeros(n_coord, dtype='f8') - h5_z = np.zeros(n_coord, dtype='f8') - nodeTags = h5_file['nodeTag'][:] - # print('h5_x.dtype: ', h5_x.dtype) - - if coordinates.ndim == 1: - for i in range(0, n_coord): - h5_x[i] = coordinates[i*3] - h5_y[i] = coordinates[i*3+1] - h5_z[i] = coordinates[i*3+2] - elif coordinates.shape[1] == 1: - for i in range(0, n_coord): - h5_x[i] = coordinates[i*3, 0] - h5_y[i] = coordinates[i*3+1, 0] - h5_z[i] = coordinates[i*3+2, 0] - else: # coordinates.shape[1] == 3 - # h5_x = coordinates[:,0] # use float32 - # h5_y = coordinates[:,1] - # h5_z = coordinates[:,2] - h5_x[:] = coordinates[:,0] # use float64 - h5_y[:] = coordinates[:,1] - h5_z[:] = coordinates[:,2] - - # print('h5_x.dtype: ', h5_x.dtype) - return h5_x, h5_y, h5_z, n_coord, nodeTags +def read_coord_drm(drm_filename, verbose): + if verbose: + print('Reading coordinates from input file [%s]' % drm_filename) + + # Get the coordinates from DRM file + with h5py.File(drm_filename, 'r') as drm_file: + coordinates = drm_file['Coordinates'] + isboundary = drm_file['Is Boundary Node'][:] + + if coordinates.ndim == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('Coordinates dataset length must be divisible by 3 for 1D input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('Coordinates dataset length must be divisible by 3 for single-column input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 3: + n_coord = coordinates.shape[0] + else: + raise ValueError('Coordinates dataset must have shape (3*n,), (3*n, 1), or (n, 3)') + + drm_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on + drm_y = np.zeros(n_coord, dtype='f8') + drm_z = np.zeros(n_coord, dtype='f8') + + if coordinates.ndim == 1: + for i in range(0, n_coord): + drm_x[i] = coordinates[i*3] + drm_y[i] = coordinates[i*3+1] + drm_z[i] = coordinates[i*3+2] + elif coordinates.shape[1] == 1: + for i in range(0, n_coord): + drm_x[i] = coordinates[i*3, 0] + drm_y[i] = coordinates[i*3+1, 0] + drm_z[i] = coordinates[i*3+2, 0] + else: # coordinates.shape[1] == 3 + drm_x[:] = coordinates[:,0] + drm_y[:] = coordinates[:,1] + drm_z[:] = coordinates[:,2] + + return drm_x, drm_y, drm_z, n_coord, isboundary + +def read_coord_h5(h5_filename, verbose): + if verbose: + print('Reading coordinates from input file [%s]' % h5_filename) + + # Get the coordinates from h5 file + with h5py.File(h5_filename, 'r') as h5_file: + coordinates = h5_file['coordinate'] + if coordinates.ndim == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('coordinate dataset length must be divisible by 3 for 1D input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 1: + if coordinates.shape[0] % 3 != 0: + raise ValueError('coordinate dataset length must be divisible by 3 for single-column input') + n_coord = coordinates.shape[0] // 3 + elif coordinates.shape[1] == 3: + n_coord = coordinates.shape[0] + else: + raise ValueError('coordinate dataset must have shape (3*n,), (3*n, 1), or (n, 3)') + + h5_x = np.zeros(n_coord, dtype='f8') # use float64 to avoid rounding-off error during crd manipulation (rotation, etc) later on + h5_y = np.zeros(n_coord, dtype='f8') + h5_z = np.zeros(n_coord, dtype='f8') + nodeTags = h5_file['nodeTag'][:] + # print('h5_x.dtype: ', h5_x.dtype) + + if coordinates.ndim == 1: + for i in range(0, n_coord): + h5_x[i] = coordinates[i*3] + h5_y[i] = coordinates[i*3+1] + h5_z[i] = coordinates[i*3+2] + elif coordinates.shape[1] == 1: + for i in range(0, n_coord): + h5_x[i] = coordinates[i*3, 0] + h5_y[i] = coordinates[i*3+1, 0] + h5_z[i] = coordinates[i*3+2, 0] + else: # coordinates.shape[1] == 3 + # h5_x = coordinates[:,0] # use float32 + # h5_y = coordinates[:,1] + # h5_z = coordinates[:,2] + h5_x[:] = coordinates[:,0] # use float64 + h5_y[:] = coordinates[:,1] + h5_z[:] = coordinates[:,2] + + # print('h5_x.dtype: ', h5_x.dtype) + return h5_x, h5_y, h5_z, n_coord, nodeTags # changed ref coord as just offsets def convert_to_essi_coord(coord_sys, from_x, from_y, from_z, ref_essi_xyz): @@ -473,111 +465,111 @@ def rotate_coords_ops_xyplane(x, y, z, rotate_angle, ref_coord=[0,0,0]): return rotated_xyz[:,0], rotated_xyz[:,1], rotated_xyz[:,2] -def get_essi_meta(ssi_fname, verbose): - # Get parameter values from HDF5 data - with h5py.File(ssi_fname, 'r') as essiout: - h = essiout['ESSI xyz grid spacing'][0] - x0 = essiout['ESSI xyz origin'][0] - y0 = essiout['ESSI xyz origin'][1] - z0 = essiout['ESSI xyz origin'][2] - t0 = essiout['time start'][0] - dt = essiout['timestep'][0] - nt = essiout['vel_0 ijk layout'].shape[0] - nx = essiout['vel_0 ijk layout'].shape[1] - ny = essiout['vel_0 ijk layout'].shape[2] - nz = essiout['vel_0 ijk layout'].shape[3] - t1 = t0 + dt*(nt-1) - timeseq = np.linspace(t0, t1, nt) - # print('dt, t0, t1, timeseq =', dt, t0, t1, timeseq) - - bTopo = False - zmin, zmax = z0, z0+(nz-1)*h - if 'z coordinates' in essiout: - bTopo = True - zcoords = essiout['z coordinates'] - chunk_x = zcoords.chunks[0] if zcoords.chunks else nx - chunk_y = zcoords.chunks[1] if zcoords.chunks else ny - zmin = np.inf - zmax = -np.inf - for ix in range(0, nx, chunk_x): - for iy in range(0, ny, chunk_y): - upper = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), 0] - lower = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), -1] - zmin = min(zmin, float(np.min(upper))) - zmax = max(zmax, float(np.max(lower))) - - return x0, y0, z0, h, nx, ny, nz, nt, dt, timeseq, bTopo, zmin, zmax - - -def get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose): - ''' when topography is considered, get z array location based on its physical - z coordinates in the interpolated vertical profile between the upper and lower z - interfaces - ''' - with h5py.File(ssi_fname, 'r') as essiout: - h = essiout['ESSI xyz grid spacing'][0] - _, nx, ny, nz = essiout['vel_0 ijk layout'].shape - # print("essiout['vel_0 ijk layout'].shape: ", essiout['vel_0 ijk layout'].shape) - - if nx < 2 or ny < 2 or nz < 2: - raise ValueError('Topography-aware z interpolation requires at least 2 grid points in x, y, and z.') - - # *) x, y indices - # start - coord_x0 = np.floor(coord_x).astype(int) - coord_y0 = np.floor(coord_y).astype(int) - coord_x0 = np.clip(coord_x0, 0, nx-2) - coord_y0 = np.clip(coord_y0, 0, ny-2) - # end - coord_x1 = coord_x0 + 1 - coord_y1 = coord_y0 + 1 - coord_x1 = np.minimum(coord_x1, nx-1) - coord_y1 = np.minimum(coord_y1, ny-1) - # get distinct x, y values - x = np.array(sorted(set(coord_x0) | set(coord_x1)), dtype=int) - y = np.array(sorted(set(coord_y0) | set(coord_y1)), dtype=int) - x_idx = x - x[0] - y_idx = y - y[0] - - # *) interpolated z profile upper and lower z physical coordinates - zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,0] # upper interface - z_upper = zcrds[np.ix_(x_idx, y_idx)] - z_upper_interp = interpn((x, y), z_upper, np.c_[coord_x, coord_y]) - - zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,-1] # lower interface - z_lower = zcrds[np.ix_(x_idx, y_idx)] - z_lower_interp = interpn((x, y), z_lower, np.c_[coord_x, coord_y]) - - z_span = z_lower_interp - z_upper_interp - zero_span = np.isclose(z_span, 0.0) - if np.any(zero_span): - iz = np.flatnonzero(zero_span)[0] - raise ValueError( - f'Error getting z array location: ' - f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' - f'has zero vertical span in the topography profile!' - ) - - z_min = np.minimum(z_upper_interp, z_lower_interp) - z_max = np.maximum(z_upper_interp, z_lower_interp) - tol = max(abs(h) * 1e-6, 1e-8) - out_of_range = (user_essi_z < z_min - tol) | (user_essi_z > z_max + tol) - if np.any(out_of_range): - iz = np.flatnonzero(out_of_range)[0] - raise ValueError( - f'Error getting z array location: ' - f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' - f'not within SW4 domain!' - ) - - # The vertical profile is linear between the interpolated upper/lower interfaces. - coord_z = (user_essi_z - z_upper_interp) * (nz - 1) / z_span - coord_z = np.clip(coord_z, 0.0, nz - 1) - - if verbose: - print('coord_z consider topography:', coord_z) - - return coord_z +def get_essi_meta(ssi_fname, verbose): + # Get parameter values from HDF5 data + with h5py.File(ssi_fname, 'r') as essiout: + h = essiout['ESSI xyz grid spacing'][0] + x0 = essiout['ESSI xyz origin'][0] + y0 = essiout['ESSI xyz origin'][1] + z0 = essiout['ESSI xyz origin'][2] + t0 = essiout['time start'][0] + dt = essiout['timestep'][0] + nt = essiout['vel_0 ijk layout'].shape[0] + nx = essiout['vel_0 ijk layout'].shape[1] + ny = essiout['vel_0 ijk layout'].shape[2] + nz = essiout['vel_0 ijk layout'].shape[3] + t1 = t0 + dt*(nt-1) + timeseq = np.linspace(t0, t1, nt) + # print('dt, t0, t1, timeseq =', dt, t0, t1, timeseq) + + bTopo = False + zmin, zmax = z0, z0+(nz-1)*h + if 'z coordinates' in essiout: + bTopo = True + zcoords = essiout['z coordinates'] + chunk_x = zcoords.chunks[0] if zcoords.chunks else nx + chunk_y = zcoords.chunks[1] if zcoords.chunks else ny + zmin = np.inf + zmax = -np.inf + for ix in range(0, nx, chunk_x): + for iy in range(0, ny, chunk_y): + upper = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), 0] + lower = zcoords[ix:min(ix + chunk_x, nx), iy:min(iy + chunk_y, ny), -1] + zmin = min(zmin, float(np.min(upper))) + zmax = max(zmax, float(np.max(lower))) + + return x0, y0, z0, h, nx, ny, nz, nt, dt, timeseq, bTopo, zmin, zmax + + +def get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose): + ''' when topography is considered, get z array location based on its physical + z coordinates in the interpolated vertical profile between the upper and lower z + interfaces + ''' + with h5py.File(ssi_fname, 'r') as essiout: + h = essiout['ESSI xyz grid spacing'][0] + _, nx, ny, nz = essiout['vel_0 ijk layout'].shape + # print("essiout['vel_0 ijk layout'].shape: ", essiout['vel_0 ijk layout'].shape) + + if nx < 2 or ny < 2 or nz < 2: + raise ValueError('Topography-aware z interpolation requires at least 2 grid points in x, y, and z.') + + # *) x, y indices + # start + coord_x0 = np.floor(coord_x).astype(int) + coord_y0 = np.floor(coord_y).astype(int) + coord_x0 = np.clip(coord_x0, 0, nx-2) + coord_y0 = np.clip(coord_y0, 0, ny-2) + # end + coord_x1 = coord_x0 + 1 + coord_y1 = coord_y0 + 1 + coord_x1 = np.minimum(coord_x1, nx-1) + coord_y1 = np.minimum(coord_y1, ny-1) + # get distinct x, y values + x = np.array(sorted(set(coord_x0) | set(coord_x1)), dtype=int) + y = np.array(sorted(set(coord_y0) | set(coord_y1)), dtype=int) + x_idx = x - x[0] + y_idx = y - y[0] + + # *) interpolated z profile upper and lower z physical coordinates + zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,0] # upper interface + z_upper = zcrds[np.ix_(x_idx, y_idx)] + z_upper_interp = interpn((x, y), z_upper, np.c_[coord_x, coord_y]) + + zcrds = essiout['z coordinates'][x[0]:x[-1]+1,y[0]:y[-1]+1,-1] # lower interface + z_lower = zcrds[np.ix_(x_idx, y_idx)] + z_lower_interp = interpn((x, y), z_lower, np.c_[coord_x, coord_y]) + + z_span = z_lower_interp - z_upper_interp + zero_span = np.isclose(z_span, 0.0) + if np.any(zero_span): + iz = np.flatnonzero(zero_span)[0] + raise ValueError( + f'Error getting z array location: ' + f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' + f'has zero vertical span in the topography profile!' + ) + + z_min = np.minimum(z_upper_interp, z_lower_interp) + z_max = np.maximum(z_upper_interp, z_lower_interp) + tol = max(abs(h) * 1e-6, 1e-8) + out_of_range = (user_essi_z < z_min - tol) | (user_essi_z > z_max + tol) + if np.any(out_of_range): + iz = np.flatnonzero(out_of_range)[0] + raise ValueError( + f'Error getting z array location: ' + f'({h*coord_x[iz]:.2f}, {h*coord_y[iz]:.2f}, {user_essi_z[iz]:.2f}) ' + f'not within SW4 domain!' + ) + + # The vertical profile is linear between the interpolated upper/lower interfaces. + coord_z = (user_essi_z - z_upper_interp) * (nz - 1) / z_span + coord_z = np.clip(coord_z, 0.0, nz - 1) + + if verbose: + print('coord_z consider topography:', coord_z) + + return coord_z def get_essi_data_btw_step(ssi_fname, start, end, verbose): @@ -727,28 +719,28 @@ def write_to_hdf5_range_1d(h5_fname, gname, dname, data, mystart, myend): dset[mystart:myend] = data[:] h5file.close() -def write_to_hdf5_range_2d(h5_fname, gname, dname, data, mystart, myend): - h5file = h5py.File(h5_fname, 'r+') - if gname == '/': - dset = h5file[dname] +def write_to_hdf5_range_2d(h5_fname, gname, dname, data, mystart, myend): + h5file = h5py.File(h5_fname, 'r+') + if gname == '/': + dset = h5file[dname] else: grp = h5file[gname] dset = grp[dname] - #print('mystart=%d, myend=%d' %(mystart, myend)) - dset[mystart:myend,:] = data[:,:] - h5file.close() - -def get_flat_coord_range(coord_offset, ncoord, width=3): - return coord_offset * width, (coord_offset + ncoord) * width - -def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): - tstep = tsteprange.step - nstep = len(tsteprange) - dt = tstep*essi_dt - tend = max(nstep - 1, 0) * dt - - h5file = h5py.File(h5_fname, 'w') + #print('mystart=%d, myend=%d' %(mystart, myend)) + dset[mystart:myend,:] = data[:,:] + h5file.close() + +def get_flat_coord_range(coord_offset, ncoord, width=3): + return coord_offset * width, (coord_offset + ncoord) * width + +def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): + tstep = tsteprange.step + nstep = len(tsteprange) + dt = tstep*essi_dt + tend = max(nstep - 1, 0) * dt + + h5file = h5py.File(h5_fname, 'w') data_grp = h5file.create_group('DRM_Data') data_location = np.zeros(ncoord, dtype='i4') for i in range(0, ncoord): @@ -762,22 +754,22 @@ def create_hdf5_opensees(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc dset = data_grp.create_dataset('displacement', (ncoord*3, nstep), dtype='f4') dset = data_grp.create_dataset('data_location', data=data_location, dtype='i4') dset = data_grp.create_dataset(extra_dname, (ncoord,), dtype='i4') - dset = data_grp.create_dataset('xyz', (ncoord, 3), dtype='f4') - - data_grp = h5file.create_group('DRM_Metadata') - dset = data_grp.create_dataset('dt', data=dt, dtype='f8') - dset = data_grp.create_dataset('tstart', data=0.0, dtype='f8') - dset = data_grp.create_dataset('tend', data=tend, dtype='f8') + dset = data_grp.create_dataset('xyz', (ncoord, 3), dtype='f4') + + data_grp = h5file.create_group('DRM_Metadata') + dset = data_grp.create_dataset('dt', data=dt, dtype='f8') + dset = data_grp.create_dataset('tstart', data=0.0, dtype='f8') + dset = data_grp.create_dataset('tend', data=tend, dtype='f8') h5file.close() -def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): - - h5file = h5py.File(h5_fname, 'w') - tstep = tsteprange.step - nstep = len(tsteprange) - dt = tstep*essi_dt - tend = max(nstep - 1, 0) * dt +def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname): + + h5file = h5py.File(h5_fname, 'w') + tstep = tsteprange.step + nstep = len(tsteprange) + dt = tstep*essi_dt + tend = max(nstep - 1, 0) * dt print('Create HDF5 file with ', ncoord, ' coordinates, ', nstep, ' steps') if gen_vel: @@ -788,11 +780,11 @@ def create_hdf5_csv(h5_fname, ncoord, tsteprange, essi_dt, gen_vel, gen_acc, gen dset = h5file.create_dataset('displacement', (ncoord*3, nstep), dtype='f4') dset = h5file.create_dataset(extra_dname, (ncoord,), dtype='i4') - dset = h5file.create_dataset('xyz', (ncoord, 3), dtype='f4') - - dset = h5file.create_dataset('dt', data=dt, dtype='f8') - dset = h5file.create_dataset('tstart', data=0.0, dtype='f8') - dset = h5file.create_dataset('tend', data=tend, dtype='f8') + dset = h5file.create_dataset('xyz', (ncoord, 3), dtype='f4') + + dset = h5file.create_dataset('dt', data=dt, dtype='f8') + dset = h5file.create_dataset('tstart', data=0.0, dtype='f8') + dset = h5file.create_dataset('tend', data=tend, dtype='f8') h5file.close() @@ -815,25 +807,25 @@ def create_hdf5_essi(h5_fname, ncoord, nstep, dt, gen_vel, gen_acc, gen_dis, ext for dname in ('Velocity', 'Accelerations', 'Displacements', 'Time'): if dname in h5file: del h5file[dname] - - if gen_vel: - dset = h5file.create_dataset('Velocity', (ncoord*3, nstep), dtype='f4') + + if gen_vel: + dset = h5file.create_dataset('Velocity', (ncoord*3, nstep), dtype='f4') if gen_acc: dset = h5file.create_dataset('Accelerations', (ncoord*3, nstep), dtype='f4') if gen_dis: dset = h5file.create_dataset('Displacements', (ncoord*3, nstep), dtype='f4') - timeseq = np.arange(nstep, dtype='f8') * dt - - h5file.create_dataset('Time', data=timeseq, dtype='f8') + timeseq = np.arange(nstep, dtype='f8') * dt + + h5file.create_dataset('Time', data=timeseq, dtype='f8') h5file.close() def coord_to_chunkid(x0, y0, z0, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, nx, ny, nz): - # Due to allocation of neighbor nodes, the indices may exceed the original dataset - # boundary and thus need to be corrected so that the correct nearest chunk can be fetched later - x, y, z = min(x0, nx-1), min(y0, ny-1), min(z0, nz-1) - + # Due to allocation of neighbor nodes, the indices may exceed the original dataset + # boundary and thus need to be corrected so that the correct nearest chunk can be fetched later + x, y, z = min(x0, nx-1), min(y0, ny-1), min(z0, nz-1) + val = int(np.floor(x/chk_x)*nchk_y*nchk_z + np.floor(y/chk_y)*nchk_z + np.floor(z/chk_z)) #print('coord_to_chunkid:', x, y, z, chk_x, chk_y, chk_z, nchk_x, nchk_y, nchk_z, val) return val @@ -972,23 +964,23 @@ def linear_interp(data_dict, x, y, z): return result def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values_input, extra_dname, output_format): - if zeroMotionDir is None or pd.isna(zeroMotionDir): - zeroMotionDir = 'None' - elif isinstance(zeroMotionDir, str): - zeroMotionDir = zeroMotionDir.strip() - if zeroMotionDir == '' or zeroMotionDir.lower() in ('none', 'no'): - zeroMotionDir = 'None' - else: - zeroMotionDir = str(zeroMotionDir).strip() - - if zeroMotionDir != 'None' and zeroMotionDir.upper() not in ('X', 'Y', 'Z'): - raise ValueError(f'Unsupported zeroMotionDir value: {zeroMotionDir}') - - # Read ESSI metadata - essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, essi_nt, essi_dt, essi_timeseq, bTopo, zmin, zmax = get_essi_meta(ssi_fname, verbose) - essi_x_len_max = (essi_nx-1) * essi_h - essi_y_len_max = (essi_ny-1) * essi_h - essi_z_len_max = zmax - zmin + if zeroMotionDir is None or pd.isna(zeroMotionDir): + zeroMotionDir = 'None' + elif isinstance(zeroMotionDir, str): + zeroMotionDir = zeroMotionDir.strip() + if zeroMotionDir == '' or zeroMotionDir.lower() in ('none', 'no'): + zeroMotionDir = 'None' + else: + zeroMotionDir = str(zeroMotionDir).strip() + + if zeroMotionDir != 'None' and zeroMotionDir.upper() not in ('X', 'Y', 'Z'): + raise ValueError(f'Unsupported zeroMotionDir value: {zeroMotionDir}') + + # Read ESSI metadata + essi_x0, essi_y0, essi_z0, essi_h, essi_nx, essi_ny, essi_nz, essi_nt, essi_dt, essi_timeseq, bTopo, zmin, zmax = get_essi_meta(ssi_fname, verbose) + essi_x_len_max = (essi_nx-1) * essi_h + essi_y_len_max = (essi_ny-1) * essi_h + essi_z_len_max = zmax - zmin # Start and end time step if start_t > -1e-6 and end_t > -1e-6: @@ -1027,13 +1019,13 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # to get 2D motions for 2D models, modify input node crds to enforce same motion across the width # Note: this should be done before rotation, motion zero-out in the out-of-plane direction will be done later # TODO: here we use the middle crd in that direction by default - user_x, user_y, user_z = user_x0, user_y0, user_z0 - # print(f'zeroMotionDir={zeroMotionDir}') - if zeroMotionDir != 'None': - middleCrd = None - if zeroMotionDir.upper() == 'X': - middleCrd = (np.amin(user_x0) + np.amax(user_x0)) / 2. - user_x = np.full(user_x0.shape, middleCrd) + user_x, user_y, user_z = user_x0, user_y0, user_z0 + # print(f'zeroMotionDir={zeroMotionDir}') + if zeroMotionDir != 'None': + middleCrd = None + if zeroMotionDir.upper() == 'X': + middleCrd = (np.amin(user_x0) + np.amax(user_x0)) / 2. + user_x = np.full(user_x0.shape, middleCrd) elif zeroMotionDir.upper() == 'Y': middleCrd = (np.amin(user_y0) + np.amax(user_y0)) / 2. user_y = np.full(user_y0.shape, middleCrd) @@ -1048,8 +1040,6 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # rotate/transform only when rotate_angle is other than 0 (default min difference is 1e-2) b_rotate = np.where(abs(rotate_angle) > 1e-2, True, False) if b_rotate: - if verbose and mpi_rank == 0: - print(f'Rotate coordinates by {rotate_angle} degree around vertical axis (z) in the xy plane') user_x, user_y, user_z = rotate_coords_ops_xyplane(user_x, user_y, user_z, rotate_angle) # Convert user coordinate to sw4 coordinate, relative to ESSI domain (subset of SW4 domain) @@ -1066,8 +1056,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # tmpfile.write('(%d, %d, %d)' % (user_essi_x[i], user_essi_y[i], user_essi_z[i])) # print('user_essi_x.dtype: ', user_essi_x.dtype) # print('user_x0.dtype: ', user_x0.dtype) - # np.savetxt('tmp_user.txt', np.c_[user_x, user_y, user_z], fmt='%.6f') - # np.savetxt('tmp_essi.txt', np.c_[user_essi_x, user_essi_y, user_essi_z], fmt='%.6f') + # np.savetxt('tmp.txt', np.c_[user_essi_x, user_essi_y, user_essi_z], fmt='%.6f') # Plot if mpi_rank == 0: @@ -1078,26 +1067,26 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use print('Only generate the plots of input nodes') exit(0) - # Check if all node coordinates are within the sw4 domain - if np.min(user_essi_x) < essi_x0 or np.max(user_essi_x) > essi_x0+essi_x_len_max or \ - np.min(user_essi_y) < essi_y0 or np.max(user_essi_y) > essi_y0+essi_y_len_max or \ - np.min(user_essi_z) < zmin or np.max(user_essi_z) > zmax: - if mpi_rank == 0: - print('Error: all node coordinates (after rotation) should be within the sw4 domain for extracting the motion') - print('while:') - print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) - print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) - print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) - - debugfile = save_path + '/user_essi_xyz.npy' - print('\tcheck user_essi_xyz (after rotation) in file \'{}\''.format(debugfile)) - np.save(debugfile, np.c_[user_essi_x, user_essi_y, user_essi_z]) - raise SystemExit(1) + # Check if all node coordinates are within the sw4 domain + if np.min(user_essi_x) < essi_x0 or np.max(user_essi_x) > essi_x0+essi_x_len_max or \ + np.min(user_essi_y) < essi_y0 or np.max(user_essi_y) > essi_y0+essi_y_len_max or \ + np.min(user_essi_z) < zmin or np.max(user_essi_z) > zmax: + if mpi_rank == 0: + print('Error: all node coordinates (after rotation) should be within the sw4 domain for extracting the motion') + print('while:') + print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) + print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) + print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) + + debugfile = save_path + '/user_essi_xyz.npy' + print('\tcheck user_essi_xyz (after rotation) in file \'{}\''.format(debugfile)) + np.save(debugfile, np.c_[user_essi_x, user_essi_y, user_essi_z]) + raise SystemExit(1) - if verbose and mpi_rank == 0: - print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) - print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) - print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) + if verbose and mpi_rank == 0: + print('\t','Min/Max SW4 x:',essi_x0,essi_x0+essi_x_len_max,'Min/Max user x:',np.min(user_essi_x),np.max(user_essi_x)) + print('\t','Min/Max SW4 y:',essi_y0,essi_y0+essi_y_len_max,'Min/Max user y:',np.min(user_essi_y),np.max(user_essi_y)) + print('\t','Min/Max SW4 z:',zmin,zmax,'Min/Max user z:',np.min(user_essi_z),np.max(user_essi_z)) # if mpi_rank == 0: # print('while user_essi_xyz (after rotation) is:\n', np.c_[user_essi_x, user_essi_y, user_essi_z]) @@ -1105,28 +1094,28 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Convert to array location (spacing is 1), floating-point coord_x = (user_essi_x - essi_x0) / essi_h - coord_y = (user_essi_y - essi_y0) / essi_h - if not bTopo: - coord_z = (user_essi_z - essi_z0) / essi_h - else: - coord_z = get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose) + coord_y = (user_essi_y - essi_y0) / essi_h + if not bTopo: + coord_z = (user_essi_z - essi_z0) / essi_h + else: + coord_z = get_coordz_topo(ssi_fname, coord_x, coord_y, user_essi_z, verbose) - # Check if we actually need interpolation - # ghost_cell = 0 - # do_interp = True - do_interp = False - coord_x_rounded = np.rint(coord_x) - coord_y_rounded = np.rint(coord_y) - coord_z_rounded = np.rint(coord_z) - for nid in range(0, n_coord): - # if user_essi_x[nid] % essi_h != 0 or user_essi_y[nid] % essi_h != 0 or user_essi_z[nid] % essi_h != 0: - if not math.isclose(coord_x[nid], coord_x_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_y[nid], coord_y_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_z[nid], coord_z_rounded[nid], abs_tol=1e-8): - do_interp = True - # ghost_cell = 1 - if verbose and mpi_rank == 0: - print(f'user_essi_x[{nid}], user_essi_y[{nid}], user_essi_z[{nid}]: {user_essi_x[nid]:.4f}, {user_essi_y[nid]:.4f}, {user_essi_z[nid]:.4f}') - # print(f'coord_x[{nid}], coord_y[{nid}], coord_z[{nid}]: {coord_x[nid]:.4e}, {coord_y[nid]:.4e}, {coord_z[nid]:.4e}') - break + # Check if we actually need interpolation + # ghost_cell = 0 + # do_interp = True + do_interp = False + coord_x_rounded = np.rint(coord_x) + coord_y_rounded = np.rint(coord_y) + coord_z_rounded = np.rint(coord_z) + for nid in range(0, n_coord): + # if user_essi_x[nid] % essi_h != 0 or user_essi_y[nid] % essi_h != 0 or user_essi_z[nid] % essi_h != 0: + if not math.isclose(coord_x[nid], coord_x_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_y[nid], coord_y_rounded[nid], abs_tol=1e-8) or not math.isclose(coord_z[nid], coord_z_rounded[nid], abs_tol=1e-8): + do_interp = True + # ghost_cell = 1 + if verbose and mpi_rank == 0: + print(f'user_essi_x[{nid}], user_essi_y[{nid}], user_essi_z[{nid}]: {user_essi_x[nid]:.4f}, {user_essi_y[nid]:.4f}, {user_essi_z[nid]:.4f}') + # print(f'coord_x[{nid}], coord_y[{nid}], coord_z[{nid}]: {coord_x[nid]:.4e}, {coord_y[nid]:.4e}, {coord_z[nid]:.4e}') + break if mpi_rank == 0: if do_interp: print('Use spline interpolation.') @@ -1473,7 +1462,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use # Write coordinates and boundary nodes (file created previously), in serial with baton passing comm.Barrier() - if output_format == "opensees": + if output_format == "opensees": if mpi_rank == 0: create_hdf5_opensees(output_fname, n_coord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname) @@ -1505,7 +1494,7 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if mpi_rank != mpi_size-1: comm.send(my_ncoord, dest=mpi_rank+1, tag=11) - elif output_format == "point": + elif output_format == "point": if mpi_rank == 0: create_hdf5_csv(output_fname, n_coord, tsteprange, essi_dt, gen_vel, gen_acc, gen_dis, extra_dname) @@ -1537,30 +1526,30 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use if mpi_rank != mpi_size-1: comm.send(my_ncoord, dest=mpi_rank+1, tag=11) - elif output_format == "essi": - coords_start, coords_end = get_flat_coord_range(my_offset, my_ncoord[0]) - if mpi_rank == 0: - create_hdf5_essi(output_fname, n_coord, nsteps, dt, gen_vel, gen_acc, gen_dis, extra_dname) - # Write to the template file - if my_ncoord[0] > 0: - write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) + elif output_format == "essi": + coords_start, coords_end = get_flat_coord_range(my_offset, my_ncoord[0]) + if mpi_rank == 0: + create_hdf5_essi(output_fname, n_coord, nsteps, dt, gen_vel, gen_acc, gen_dis, extra_dname) + # Write to the template file + if my_ncoord[0] > 0: + write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) - if gen_acc: - write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) - if gen_dis: + if gen_acc: + write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) + if gen_dis: write_to_hdf5_range(output_fname, '/', 'Displacements', output_dis_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_vel: write_to_hdf5_range(output_fname, '/', 'Velocity', output_vel_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if mpi_size > 1: comm.send(my_ncoord, dest=1, tag=111) - else: - data = comm.recv(source=mpi_rank-1, tag=111) - if my_ncoord[0] > 0: - write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) + else: + data = comm.recv(source=mpi_rank-1, tag=111) + if my_ncoord[0] > 0: + write_to_hdf5_range_1d(output_fname, '/', 'Coordinates', my_user_coordinates.reshape(my_ncoord[0]*3), coords_start, coords_end) write_to_hdf5_range_1d(output_fname, '/', extra_dname, extra_values, my_offset, my_offset+my_ncoord[0]) - if gen_acc: - write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) + if gen_acc: + write_to_hdf5_range(output_fname, '/', 'Accelerations', output_acc_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_dis: write_to_hdf5_range(output_fname, '/', 'Displacements', output_dis_all[:,tsteprange], my_offset*3, (my_offset+my_ncoord[0])*3) if gen_vel: @@ -1577,49 +1566,49 @@ def generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, use print('Rank', mpi_rank, 'Finished writing data') return -def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input DRM [%s]' %drm_fname) - print('Input ESSI [%s]' %ssi_fname) - - coord_sys = ['y', 'x', '-z'] - output_format = resolve_output_mode('drm', requested_output_mode) - - # original unrotated node coordinates +def convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input DRM [%s]' %drm_fname) + print('Input ESSI [%s]' %ssi_fname) + + coord_sys = ['y', 'x', '-z'] + output_format = resolve_output_mode('drm', requested_output_mode) + + # original unrotated node coordinates user_x0, user_y0, user_z0, n_coord, boundary_flags = read_coord_drm(drm_fname, verbose) extra_values = np.asarray(boundary_flags).reshape(-1) - if verbose and mpi_rank == 0: - print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x0[0], user_y0[0], user_z0[0], user_x0[-1], user_y0[-1], user_z0[-1])) - print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x0), np.max(user_x0), np.min(user_y0), np.max(user_y0), np.min(user_z0), np.max(user_z0)) ) - - if output_format == 'opensees': - gen_vel = False - gen_dis = True - gen_acc = True - extra_dname = 'internal' - else: - gen_vel = True - gen_dis = True - gen_acc = True - extra_dname = 'Is Boundary Node' - - output_fname = get_output_filename( - 'drm', output_format, save_path, drm_fname, explicit_output_mode=requested_output_mode is not None - ) - if output_format == 'essi': - if mpi_rank == 0: - prepare_essi_output_file(drm_fname, output_fname) - MPI.COMM_WORLD.Barrier() - + if verbose and mpi_rank == 0: + print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x0[0], user_y0[0], user_z0[0], user_x0[-1], user_y0[-1], user_z0[-1])) + print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x0), np.max(user_x0), np.min(user_y0), np.max(user_y0), np.min(user_z0), np.max(user_z0)) ) + + if output_format == 'opensees': + gen_vel = False + gen_dis = True + gen_acc = True + extra_dname = 'internal' + else: + gen_vel = True + gen_dis = True + gen_acc = True + extra_dname = 'Is Boundary Node' + + output_fname = get_output_filename( + 'drm', output_format, save_path, drm_fname, explicit_output_mode=requested_output_mode is not None + ) + if output_format == 'essi': + if mpi_rank == 0: + prepare_essi_output_file(drm_fname, output_fname) + MPI.COMM_WORLD.Barrier() + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values, extra_dname, output_format) - - return - -def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input h5 [%s]' %h5_fname) - print('Input ESSI [%s]' %ssi_fname) - + + return + +def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input h5 [%s]' %h5_fname) + print('Input ESSI [%s]' %ssi_fname) + coord_sys = ['y', 'x', '-z'] output_format = resolve_output_mode('h5', requested_output_mode) gen_vel = output_format != 'opensees' @@ -1631,22 +1620,22 @@ def convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, user_x0, user_y0, user_z0, n_coord, node_tags = read_coord_h5(h5_fname, verbose) n_coord = len(node_tags) - if mpi_rank == 0: - print('Generating motions for %i nodes...' % (n_coord)) - - output_fname = get_output_filename( - 'h5', output_format, save_path, h5_fname, explicit_output_mode=requested_output_mode is not None - ) - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) - - return - -def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): - if mpi_rank == 0: - print('Input CSV [%s]' %csv_fname) - print('Input ESSI [%s]' %ssi_fname) - + if mpi_rank == 0: + print('Generating motions for %i nodes...' % (n_coord)) + + output_fname = get_output_filename( + 'h5', output_format, save_path, h5_fname, explicit_output_mode=requested_output_mode is not None + ) + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) + + return + +def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, requested_output_mode=None): + if mpi_rank == 0: + print('Input CSV [%s]' %csv_fname) + print('Input ESSI [%s]' %ssi_fname) + coord_sys = ['y', 'x', '-z'] output_format = resolve_output_mode('csv', requested_output_mode) gen_vel = output_format != 'opensees' @@ -1666,32 +1655,32 @@ def convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tste user_y0[i] = df.loc[i, 'y'] user_z0[i] = df.loc[i, 'z'] - if mpi_rank == 0: - print('Generating motions for %i nodes...' % (n_coord)) - - output_fname = get_output_filename( - 'csv', output_format, save_path, csv_fname, explicit_output_mode=requested_output_mode is not None - ) - - generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) - + if mpi_rank == 0: + print('Generating motions for %i nodes...' % (n_coord)) + + output_fname = get_output_filename( + 'csv', output_format, save_path, csv_fname, explicit_output_mode=requested_output_mode is not None + ) + + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x0, user_y0, user_z0, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir,gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, node_tags, extra_dname, output_format) + return def dframeToDict(dFrame): dFrame = list(dFrame.iterrows()) return {i[1].to_list()[0] : i[1].to_list()[1] for i in dFrame} -def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, ref_coord=None, requested_output_mode=None): - if mpi_rank == 0: - print('Input CSV [%s]' %csv_fname) - print('Input ESSI [%s]' %ssi_fname) - - if ref_coord is None: - ref_coord = np.zeros(3) - - output_format = resolve_output_mode('template', requested_output_mode) - - sw4ToESSI_params = dframeToDict(pd.read_csv(csv_fname)) +def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plot_only, mpi_rank, mpi_size, verbose, ref_coord=None, requested_output_mode=None): + if mpi_rank == 0: + print('Input CSV [%s]' %csv_fname) + print('Input ESSI [%s]' %ssi_fname) + + if ref_coord is None: + ref_coord = np.zeros(3) + + output_format = resolve_output_mode('template', requested_output_mode) + + sw4ToESSI_params = dframeToDict(pd.read_csv(csv_fname)) sw4_i_start = sw4ToESSI_params["sw4_i_start"] sw4_i_end = sw4ToESSI_params["sw4_i_end"] sw4_j_start = sw4ToESSI_params["sw4_j_start"] @@ -1708,29 +1697,29 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep # reference point, which is the ESSI or OpenSees origin in the SW4 coordinate system ref_coord[0] = essi_x_start ref_coord[1] = essi_y_start - ref_coord[2] = essi_z_start + ref_coord[2] = essi_z_start - coord_sys = ['y', 'x', '-z'] - if output_format == 'opensees': - gen_vel = False - gen_dis = True - gen_acc = True - extra_dname = 'internal' - else: - gen_vel = True - gen_dis = True - gen_acc = True - extra_dname = 'Is Boundary Node' - - # original unrotated node coordinates - template_file = h5py.File(template_fname) - coordinates = template_file['Coordinates'][:] - + coord_sys = ['y', 'x', '-z'] + if output_format == 'opensees': + gen_vel = False + gen_dis = True + gen_acc = True + extra_dname = 'internal' + else: + gen_vel = True + gen_dis = True + gen_acc = True + extra_dname = 'Is Boundary Node' + + # original unrotated node coordinates + template_file = h5py.File(template_fname) + coordinates = template_file['Coordinates'][:] + extra_values = np.asarray(template_file['Is Boundary Node'][:]).reshape(-1) n_coord = len(extra_values) - user_x = np.zeros(n_coord) - user_y = np.zeros(n_coord) - user_z = np.zeros(n_coord) + user_x = np.zeros(n_coord) + user_y = np.zeros(n_coord) + user_z = np.zeros(n_coord) for i in range(0, n_coord): user_x[i] = coordinates[i*3] user_y[i] = coordinates[i*3+1] @@ -1738,21 +1727,21 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep template_file.close() - if verbose and mpi_rank == 0: - print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x[0], user_y[0], user_z[0], user_x[-1], user_y[-1], user_z[-1])) - print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x), np.max(user_x), np.min(user_y), np.max(user_y), np.min(user_z), np.max(user_z)) ) - print('Start/end time', start_t, end_t) - - output_fname = get_output_filename( - 'template', output_format, save_path, template_fname, explicit_output_mode=requested_output_mode is not None - ) - if output_format == 'essi': - if mpi_rank == 0: - prepare_essi_output_file(template_fname, output_fname) - MPI.COMM_WORLD.Barrier() - + if verbose and mpi_rank == 0: + print('Done read %d coordinates, first is (%d, %d, %d), last is (%d, %d, %d)' % (n_coord, user_x[0], user_y[0], user_z[0], user_x[-1], user_y[-1], user_z[-1])) + print('x, y, z (min/max): (%.0f, %.0f), (%.0f, %.0f), (%.0f, %.0f)' % (np.min(user_x), np.max(user_x), np.min(user_y), np.max(user_y), np.min(user_z), np.max(user_z)) ) + print('Start/end time', start_t, end_t) + + output_fname = get_output_filename( + 'template', output_format, save_path, template_fname, explicit_output_mode=requested_output_mode is not None + ) + if output_format == 'essi': + if mpi_rank == 0: + prepare_essi_output_file(template_fname, output_fname) + MPI.COMM_WORLD.Barrier() + generate_acc_dis_time(ssi_fname, coord_sys, ref_coord, user_x, user_y, user_z, n_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, gen_vel, gen_acc, gen_dis, verbose, plot_only, output_fname, mpi_rank, mpi_size, extra_values, extra_dname, output_format) - return + return if __name__ == "__main__": @@ -1774,12 +1763,12 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep start_t=0 end_t=0 tstep = 1 - rotate_angle = 0 - zeroMotionDir = 'None' - requested_output_mode = None + rotate_angle = 0 + zeroMotionDir = 'None' + requested_output_mode = None - parser = build_arg_parser() - args = parser.parse_args() + parser = build_arg_parser() + args = parser.parse_args() if args.verbose: verbose=True @@ -1802,8 +1791,8 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep if args.template: template_fname=args.template use_template=True - if args.ssi: - ssi_fname=args.ssi + if args.ssi: + ssi_fname=args.ssi if args.reference: ref_coord[0]=args.reference[0] ref_coord[1]=args.reference[1] @@ -1814,12 +1803,12 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep tstep = args.timerange[2] # time step interval if args.rotateanlge: rotate_angle = args.rotateanlge - if args.savepath: - save_path = args.savepath - if args.zeroMotionDir: - zeroMotionDir = args.zeroMotionDir - if args.output_format: - requested_output_mode = args.output_format + if args.savepath: + save_path = args.savepath + if args.zeroMotionDir: + zeroMotionDir = args.zeroMotionDir + if args.output_format: + requested_output_mode = args.output_format comm = MPI.COMM_WORLD mpi_size = comm.Get_size() @@ -1833,31 +1822,31 @@ def convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep startTimeStr = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime(startTime)) print('Start time:', startTimeStr) - if drm_fname == '' and h5_fname == '' and csv_fname == '' and template_fname == '': - print('Error, no node coordinate input file is provided, exit...') - exit(0) + if drm_fname == '' and h5_fname == '' and csv_fname == '' and template_fname == '': + print('Error, no node coordinate input file is provided, exit...') + exit(0) if ssi_fname == '': - print('Error, no SW4 SSI output file is provided, exit...') + print('Error, no SW4 SSI output file is provided, exit...') exit(0) - if verbose and mpi_rank == 0: - print('Using ref_coord={}, start_t={}, end_t={}, tstep={}, rotate_angle={}, output_mode={} to extract motions'.format(ref_coord, start_t, end_t, tstep, rotate_angle, requested_output_mode or 'default')) - - try: - if use_drm: - convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_h5: - convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_csv and not use_template: - convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) - elif use_csv and use_template: - convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, ref_coord, requested_output_mode) - except ValueError as exc: - if mpi_rank == 0: - print(f'Error: {exc}') - raise SystemExit(1) + if verbose and mpi_rank == 0: + print('Using ref_coord={}, start_t={}, end_t={}, tstep={}, rotate_angle={}, output_mode={} to extract motions'.format(ref_coord, start_t, end_t, tstep, rotate_angle, requested_output_mode or 'default')) + + try: + if use_drm: + convert_drm(drm_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_h5: + convert_h5(h5_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_csv and not use_template: + convert_csv(csv_fname, ssi_fname, save_path, ref_coord, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, requested_output_mode) + elif use_csv and use_template: + convert_template(csv_fname, template_fname, ssi_fname, start_t, end_t, tstep, rotate_angle, zeroMotionDir, plotonly, mpi_rank, mpi_size, verbose, ref_coord, requested_output_mode) + except ValueError as exc: + if mpi_rank == 0: + print(f'Error: {exc}') + raise SystemExit(1) if mpi_rank == 0: endTime = time.time() endTimeStr = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime(endTime)) - print(f'End time: {endTimeStr} (time spent: {endTime-startTime:.2f} s)') + print(f'End time: {endTimeStr} (time spent: {endTime-startTime:.2f} s)') diff --git a/template/motion_setting.csv b/template/motion_setting.csv index e807fc6..cf3078c 100644 --- a/template/motion_setting.csv +++ b/template/motion_setting.csv @@ -1,5 +1,5 @@ nodeTag,x,y,z,essiXstart,essiYstart,essiZstart,startTime,endTime,tstep,rotationAngle,zeroMotionDir -1,0,0,0,12000,40000,24,0,0,5,0,None +1,0,0,0,12000,40000,24,0,0,5,0,None 2,8,0,0,,,,,,,, 3,8,8,0,,,,,,,, 4,0,8,0,,,,,,,,