Source code for pyNastran.converters.shabp.shabp

"""
dfeines:
 - read_shabp(shabp_filename, log=None, debug=False)
 - SHABP(log=None, debug=False)

"""
from numpy import array, zeros, arange, ones, cross
from numpy.linalg import norm  # type: ignore
from cpylog import get_logger2

from pyNastran.converters.shabp.shabp_results import ShabpOut
#from pyNastran.converters.shabp.parse_trailer import parse_trailer

[docs] def read_shabp(shabp_filename, read_special_routines=False, log=None, debug=False): """reads an S/HABP file""" model = SHABP(log=log, debug=debug) model.read_shabp(shabp_filename, read_special_routines=read_special_routines) return model
[docs] class SHABP(ShabpOut): """defines the SHABP class""" def __init__(self, log=None, debug=False): """ Initializes the SHABP object Parameters ---------- debug : bool/None; default=True used to set the logger if no logger is passed in True: logs debug/info/error messages False: logs info/error messages None: logs error messages log : logging module object / None if log is set, debug is ignored and uses the settings the logging object has """ #self.xyz = {} self.X = {} self.Y = {} self.Z = {} self.trailer = None self.component_name_to_patch = {} self.patch_to_component_num = {} self.component_to_params = {} self.component_num_to_name = {} self.component_name_to_num = {} self.log = get_logger2(log, debug=debug) self.title = '' self.header = '' self.shabp_cases = {} #def write_shabp(self, out_filename): #"""writes an S/HABP file""" #pass
[docs] def get_area_xlength_by_component(self, components=None): """gets the area and length of a set of components""" if components is None: components = self.component_name_to_patch.keys() #ncomponents = len(components) # we're not using a vector because component name is # probably more useful areas = {} lengths = {} xmax = None xmin = None for name in components: patches = self.component_name_to_patch[name] area = 0. for unused_i, ipatch in enumerate(patches): X = self.X[ipatch-1] nrows, ncols = X.shape npoints = nrows * ncols xyz = zeros((npoints, 3), dtype='float32') xyz[:, 0] = self.X[ipatch-1].ravel() xyz[:, 1] = self.Y[ipatch-1].ravel() xyz[:, 2] = self.Z[ipatch-1].ravel() if xmin is None: xmin = xyz[:, 0].min() xmax = xyz[:, 0].max() else: xmin = min(xmin, xyz[:, 0].min()) xmax = max(xmax, xyz[:, 0].max()) # TODO: can we vectorize this efficiently? for irow in range(nrows-1): for jcol in range(ncols-1): i1 = irow*ncols +jcol, i2 = (irow+1)*ncols + jcol, i3 = (irow+1)*ncols + (jcol+1), i4 = irow*ncols +(jcol+1), a = xyz[i3, :] - xyz[i1, :] b = xyz[i4, :] - xyz[i2, :] area += 0.5 * norm(cross(a, b)) areas[name] = area lengths[name] = xmax - xmin return areas, lengths
[docs] def get_area_by_component(self, components=None): """gets the area of a set of components""" if components is None: components = self.component_name_to_patch.keys() #ncomponents = len(components) # we're not using a vector because component name is # probably more useful areas = {} for name in components: patches = self.component_name_to_patch[name] area = self.get_area_by_patch(patches) areas[name] = area.sum() return areas
[docs] def get_area_by_patch(self, ipatches=None): """gets the area of a set of patches""" if ipatches is None: npatches = len(self.X) ipatches = arange(npatches) areas = zeros(len(ipatches), dtype='float32') for i, ipatch in enumerate(ipatches): X = self.X[ipatch-1] nrows, ncols = X.shape npoints = nrows * ncols xyz = zeros((npoints, 3), dtype='float32') xyz[:, 0] = self.X[ipatch-1].ravel() xyz[:, 1] = self.Y[ipatch-1].ravel() xyz[:, 2] = self.Z[ipatch-1].ravel() area = 0. # TODO: can we vectorize this efficiently? for irow in range(nrows-1): for jcol in range(ncols-1): i1 = irow*ncols +jcol, i2 = (irow+1)*ncols +(jcol), i3 = (irow+1)*ncols +(jcol+1), i4 = irow*ncols + (jcol+1), a = xyz[i3, :] - xyz[i1, :] b = xyz[i4, :] - xyz[i2, :] area += 0.5 * norm(cross(a, b)) areas[i] = area return areas
[docs] def read_shabp(self, shabp_filename, read_special_routines=False): """reads an SHABP.INP / SHABP.mk5 file""" with open(shabp_filename) as shabp_file: lines = shabp_file.readlines() if shabp_filename.lower().endswith(".geo"): i = 0 else: # this supports standard .inp and .mk5 files i = 3 self.header = lines[0:i] patches = [] while i < len(lines)-1: patch = [] # C 1001 10 3 patch 1 header = lines[i] #print "%r" % header.strip() if 1: name = header[:7].strip() #num = int(header[4:7]) #nrows = int(header[7:20]) #ncols = int(header[20:29]) #name2 = header[29:].strip() #print lines[i].strip() flag = name[4] self.log.debug(f'name={name} flag={flag}') else: stream_cross = header[5] # heat/tk A-10 symmetry = header[16] # heat/tk A-10 scale_factor = header[17] # heat/tk A-10 #nadj1 = header[18:20] # not used #nadj2 = header[20:22] # not used #nadj3 = header[22:24] # not used #nadj4 = header[24:26] # not used self.log.debug(f'stream_cross={stream_cross} symmetry={symmetry} scale_factor={scale_factor}') vis_type = int(header[26]) #if vis_type == 0: #print("inviscid_A") #if vis_type == 1: #print("neither") #if vis_type == 2: #print("viscous_A") #if vis_type == 3: #print("inviscid_B") name_old = header[:7].strip() self.log.info(f'name={name!r} stream_cross={stream_cross!r} symmetry={symmetry!r} ' f'scale_factor={scale_factor!r} vis_type={vis_type!r}') header2 = lines[i+1] xsc = header2[0:10].strip() ysc = header2[10:20].strip() zsc = header2[20:30].strip() delx = header2[30:40].strip() dely = header2[40:50].strip() delz = header2[50:60].strip() self.log.info('sc=%s,%s,%s del=%r,%r,%r' % (xsc, ysc, zsc, delx, dely, delz)) #print(lines[i].strip()) flag = name_old[4] if flag == '1': is_last_patch = True elif flag == '0': is_last_patch = False else: raise RuntimeError('last patch flag = %r; must be 0 or 1' % flag) #print(f'name={name} nrows={nrows} ncols={ncols:d} name2={name2!r}') i += 2 row = [] while 1: #-1071.0480 77.2500 -66.94202 -987.7440 77.2500 -66.94200 3 line = lines[i] #print(i, line, end='') i += 1 #t1 = int(line[30]) # STAT - status flag - p. 27 x1, y1, z1, t1 = line[:10], line[10:20], line[20:30], line[30:31] t1 = line[30] x1 = float(x1) y1 = float(y1) z1 = float(z1) t1 = int(t1) if t1 == 1: patch.append(row) row = [] row.append([x1, y1, z1]) elif t1 in [0, 2]: row.append([x1, y1, z1]) elif t1 == 3: row.append([x1, y1, z1]) patch.append(row) #assert len(row) == nrows, 'len(row)=%s nrows=%s ncols=%s' % (len(row), nrows, ncols) break else: raise RuntimeError() if t1 == 3: break x2, y2, z2, t2 = line[31:41], line[41:51], line[51:61], line[61:62] t2 = line[61] x2 = float(x2) y2 = float(y2) z2 = float(z2) t2 = int(t2) if t2 == 1: patch.append(row) row = [] row.append([x2, y2, z2]) elif t2 in [0, 2]: row.append([x2, y2, z2]) elif t2 == 3: row.append([x2, y2, z2]) patch.append(row) break else: raise RuntimeError() try: patchi = array(patch, dtype='float32') except Exception: print('patch =', patch) for i, patchi in enumerate(patch): print('i=%s n=%s patch=%s' % (i, len(patchi), patchi)) raise patches.append(patchi) if is_last_patch: #print "***last patch - lines[%i] = %r" % (i, lines[i]) break #print "lines[%i] = %r" % (i, lines[i]) self.trailer = lines[i:] self.build_patches(patches) try: self.parse_trailer(read_special_routines) except (RuntimeError, ValueError): self.log.error('failed parsing trailer')
#raise
[docs] def build_patches(self, patches): X = [] Y = [] Z = [] #XYZ = [] for patch in patches: nrows = len(patch) ncols = len(patch[0]) #xyz = zeros((nrows, ncols, 3), dtype='float32') x = zeros((nrows, ncols), dtype='float32') y = zeros((nrows, ncols), dtype='float32') z = zeros((nrows, ncols), dtype='float32') for irow, row in enumerate(patch): for icol, col in enumerate(row): x[irow, icol] = col[0] y[irow, icol] = col[1] z[irow, icol] = col[2] #xyz[irow, icol, :] = [col[0], col[1], col[2]] X.append(x) Y.append(y) Z.append(z) #XYZ.append(xyz) self.X = X self.Y = Y self.Z = Z
#self.XYZ = XYZ
[docs] def get_points_elements_regions(self): npatches = len(self.X) npoints = 0 nelements = 0 for ipatch in range(npatches): X = self.X[ipatch] nrows, ncols = X.shape npoints += nrows * ncols nelements += (nrows-1) * (ncols-1) ipoint = 0 ielement = 0 xyz = zeros((npoints, 3), dtype='float32') elements2 = zeros((nelements, 4), dtype='int32') components = ones(nelements, dtype='int32') patches = ones(nelements, dtype='int32') impact = ones(nelements, dtype='int32') shadow = ones(nelements, dtype='int32') for ipatch in range(npatches): comp_num, impact_val, shadow_val = self.get_impact_shadow(ipatch) nrows, ncols = self.X[ipatch].shape npointsi = nrows * ncols nelementsi = (nrows-1) * (ncols-1) xyz[ipoint:ipoint+npointsi, 0] = self.X[ipatch].ravel() xyz[ipoint:ipoint+npointsi, 1] = self.Y[ipatch].ravel() xyz[ipoint:ipoint+npointsi, 2] = self.Z[ipatch].ravel() #if comp_num == 0: #continue elements = [] for irow in range(nrows-1): for jcol in range(ncols-1): element = [ irow*ncols + jcol, (irow+1)*ncols + (jcol), (irow+1)*ncols + (jcol+1), irow*ncols + (jcol+1), ] elements.append(element) elements = array(elements, dtype='int32') patches[ielement:ielement+nelementsi] *= (ipatch+1) components[ielement:ielement+nelementsi] *= comp_num impact[ielement:ielement+nelementsi] *= impact_val shadow[ielement:ielement+nelementsi] *= shadow_val elements2[ielement:ielement+nelementsi, :] = elements[:, :] + ipoint #print(" ipatch=%i Cp[%i:%i]" % (ipatch+1, ielement, ielement+nelementsi)) ipoint += npointsi ielement += nelementsi return xyz, elements2, patches, components, impact, shadow
[docs] def get_impact_shadow(self, ipatch): """We may not have inviscid pressure""" comp_num = 0 impact_val = 0 shadow_val = 0 if ipatch not in self.patch_to_component_num: self.log.debug(f'skipping ipatch={ipatch} comp_num={comp_num}') return comp_num, impact_val, shadow_val #raise RuntimeError('ipatch=%s keys=%s' % ( #ipatch, sorted(self.patch_to_component_num))) comp_num = self.patch_to_component_num[ipatch] if comp_num not in self.component_num_to_name: self.log.debug(f'skipping ipatch={ipatch} comp_num={comp_num}') return comp_num, impact_val, shadow_val name = self.component_num_to_name[comp_num] patch = self.component_name_to_patch[name] if comp_num-1 not in self.component_to_params: self.log.debug(f'skipping ipatch={ipatch} comp_num={comp_num} name={name!r} patch={patch}') return comp_num, impact_val, shadow_val data = self.component_to_params[comp_num-1] #except KeyError: #name = self.component_num_to_name[comp_num] #self.log.error(f'ipatch={ipatch} comp_num={comp_num} name={name}') #print(self.component_to_params) #raise impact_val = data[0] shadow_val = data[1] self.log.debug(f'loading ipatch={ipatch} comp_num={comp_num} name={name!r} patch={patch}') return comp_num, impact_val, shadow_val
[docs] def parse_trailer(self, read_special_routines=False): """parses the case information (e.g., number of angles of attack, control surface info)""" #out = parse_trailer(self.trailer) #order, component_names, cases, components = out #print('order = %s' % order) #print('component_names = %s' % component_names) #print('cases = %s' % cases) #print('components = %s' % components) #return out #for line in self.trailer: #print line.rstrip() self.title = self.trailer[0].strip() #print('title = %r' % self.title) line2 = self.trailer[1] #print('line2 = %r' % line2[:20]) methods, summation_flag = _get_methods(line2) #print('methods =', methods) #print('summation_flag =', summation_flag) unused_npatches = line2[:2] #print('npatches =', npatches) mach_line = self.trailer[2].rstrip() mach, alt, pstag, tstag, igas, nalpha_beta, ideriv, dangle = _parse_flight_condition(mach_line) #print(f'mach={mach} alt={alt} pstag={pstag} tstag={tstag} igas={igas} ' #f'nalpha_beta={nalpha_beta} ideriv={ideriv} dangle={dangle}') ref_line = self.trailer[3].rstrip() unused_reference_quantities = _parse_reference_conditions(ref_line) #sref, cref, bref, xcg, ycg, zcg = reference_quantities #print(f'sref={sref} cref={cref} bref={bref} cg=[{xcg},{ycg},{zcg}]') i = 4 for n in range(nalpha_beta): alpha_line = self.trailer[i].rstrip() alpha, beta, roll, cdelta, qi, ri, pi = _read_alpha_beta_line(alpha_line) #print(f'alpha={alpha} beta={beta} roll={roll} cdelta={cdelta} pqr=[{qi},{ri},{pi}]') self.shabp_cases[n] = [mach, alpha, beta] i += 1 #self.getMethods() self.log.debug(f'methods = {methods}') for method in methods: line = self.trailer[i].rstrip() if method == 3: self.log.debug(f'inviscid_pressure: {line}') i, ncomponents, ifsave, params = parse_inviscid_pressure(self.trailer, line, i) self.component_to_params = params elif method == 4: i = parse_viscous(self.trailer, line, i) elif method == 5: i = parse_special_routines(self.trailer, line, i, read_special_routines=read_special_routines) else: raise NotImplementedError(method) #print "ncomps =", len(self.component_to_params) #print "keys =", sorted(self.component_to_params.keys()) #print "**lines[%i] = %s\n" % (i+1, self.trailer[i].rstrip()) #i += 1 # component names 7 comp_names_line = self.trailer[i].rstrip() #print("comp_names_line = %r" % comp_names_line) ncomps = int(comp_names_line.strip().split()[2]) i += 1 for icomp in range(ncomps): line = self.trailer[i] #print("line =", line.strip()) unused_npatches = int(line[:2]) name = line[2:].strip() line = self.trailer[i+1] + ' ' patches = [] for n in range(0, len(line), 2): #print("n =", n) ipatch = line[n:n+2].strip() if len(ipatch) == 0: break int_ipatch = int(ipatch) #print("int_ipatch =", int_ipatch) patches.append(int_ipatch) self.patch_to_component_num[int_ipatch-1] = icomp+1 #print("patches =", patches, '\n') self.component_name_to_patch[name] = patches self.component_num_to_name[icomp] = name self.component_name_to_num[name] = icomp i += 2 # 2noseconeright #3142 self.log.debug('done with trailer')
def _parse_flight_condition(mach_line): """reads a flight condition line""" mach = float(mach_line[0 :10].strip()) alt = float(mach_line[10:20].strip()) pstag = float(mach_line[20:30].strip()) tstag = float(mach_line[30:40].strip()) igas = int(mach_line[40:41].strip()) nalpha_beta = int(mach_line[41:43].strip()) # 0-no derivatives # 1: pitch static stability derivatives # 2: pitch control derivatives # 3: pitch dynamic derivatives # 4: lat/dir static stability derivatives # 5: lat/dir control derivatives # 6: lat/dir dynamic derivatives # 7: ideriv=1,2,3 (all pitch) # 8: ideriv=4,5,6 (all lat/dir) # 9: ideriv=1-6 (all) try: ideriv = int(mach_line[43].strip()) # 0-9 dangle = float(mach_line[45:55].strip()) # used for ideriv=1,2,4,5 except Exception: # old SHABP raise # old SHABP #iDeriv = 0 #dAngle = 0.0 assert igas in [0, 1], igas # 0=air, 1=helium return mach, alt, pstag, tstag, igas, nalpha_beta, ideriv, dangle def _get_methods(line): """ per vecc A-53 1 - Flow Field Analysis 2 - Shielding Analysis 3 - Inviscid Pressures 4 - Viscous Methods 5 - Special Routines """ fields = line[:20] assert len(fields) == 20, fields methods = [] for inti in line[:20]: if inti == '0': continue methodi = int(inti) methods.append(methodi) summation_flag = int(line[20]) return methods, summation_flag def _read_alpha_beta_line(alpha_line): """reads the alpha line""" #print("alpha_line =", alpha_line) alpha = float(alpha_line[0:10]) beta = float(alpha_line[10:20]) roll = float(alpha_line[20:30]) cdelta = float(alpha_line[30:40]) # unused qi = float(alpha_line[40:50]) ri = float(alpha_line[50:60]) #pi = float(alpha_line[60:70]) pi = None return alpha, beta, roll, cdelta, qi, ri, pi def _parse_reference_conditions(ref_line): """reads the reference quantities line""" sref = float(ref_line[0:10].strip()) cref = float(ref_line[10:20].strip()) bref = float(ref_line[20:30].strip()) xcg = float(ref_line[30:40].strip()) ycg = float(ref_line[40:50].strip()) zcg = float(ref_line[50:60].strip()) return sref, cref, bref, xcg, ycg, zcg
[docs] def parse_inviscid_pressure(lines, line, i): """reads inviscid pressure (method=3)""" ncomponents = int(line[0:2]) # 0 : save component force data; use new Unit 9 # 1 : save component force data; use old Unit 9 # 2 : dont save force data ifsave = int(line[2:3]) unused_title = line[6:66].strip() assert ifsave == 0, 'ifsave=%r; %s' % (ifsave, line) #print(line) #print(f'ncomponents={ncomponents} ifsave={ifsave} title={title!r}') i += 1 #return i, ncomponents, ifsave params = {} for icomponent in range(ncomponents): hinge_line = lines[i].rstrip() #print('hinge_line = %r' % hinge_line) #print "icount=",icount # get deflection angle, check for a hinge #n+=1 #print "n = ",n #print "trailer[",n+1,"] = ",trailer[n+1].strip() hinge_method = int(hinge_line[7:8]) unused_dangle = hinge_line[12:18] unused_hinge_patch = int(hinge_line[1:3]) #print(f'hinge_method={hinge_method} dangle={dangle} hinge_patch={hinge_patch}') i += 1 line2 = lines[i].rstrip() #print(line2) impact = int(line2[0:2].strip()) shadow = int(line2[2:4].strip()) iprint = int(line2[4:5].strip()) ipin = int(line2[5:6].strip()) isave = int(line2[6:7].strip()) pdata1 = float(line2[10:20].strip()) pdata2 = float(line2[20:30].strip()) pdata3 = float(line2[30:40].strip()) pdata4 = float(line2[40:50].strip()) pdata5 = float(line2[50:60].strip()) pdata6 = float(line2[60:70].strip()) #print(f'icomponent={icomponent} impact={impact} shadow={shadow} iprint={iprint}') params[icomponent] = [ impact, shadow, iprint, ipin, isave, pdata1, pdata2, pdata3, pdata4, pdata5, pdata6 ] i += 1 if hinge_method == 0: pass elif hinge_method == 3: hinge_line = lines[i].rstrip() #print('hinge_line2 = %r' % hinge_line) x1 = float(hinge_line[0:10]) y1 = float(hinge_line[10:20]) z1 = float(hinge_line[20:30]) x2 = float(hinge_line[30:40]) y2 = float(hinge_line[40:50]) z2 = float(hinge_line[50:60]) unused_hinge = (x1, y1, z1, x2, y2, z2) i += 1 #print('hinge =', hinge) else: print(hinge_line) raise NotImplementedError(hinge_method) return i, ncomponents, ifsave, params
[docs] def parse_viscous(lines, line, i): """ Parses the viscous table (method=4) 140 viscous_run (level 2) 11 1 000 5 5000 0.0000 1.0000 0.0000 1.0000 3.0000 3.0000 2 0 4 1 1 1 1 1 0 0 0 1 1 0 1 1.500000E2 3.000000E2 1.600000E1 3.000000E1 9.000000E1 5.000000E2 8.000000E-1 1.200000E1 0.000000E0 0.000000E0 21 1 ... """ line = lines[i].rstrip() # Skin Friction Base & Title card; A-73 ncomp = int(line[0:2]) unused_ifsave = line[2:3] unused_title = line[6:66].strip() #print(f'ncomp={ncomp} ifsave={ifsave} title={title}') i += 1 #print('***', line) for unused_icompi in range(ncomp): #geometry data source card; A-73 line = lines[i].rstrip() #print(line) unused_a, unused_b = line.split() unused_icomp = int(line[:3]) # 0=skin friction method cards will be read for each alpha/beta # 1=read only one skin friction method card and assume it applies to all alpha/betas for this component # 2=use the same skin friction data as for the prepvious component (skip reading data) unused_isk = line[3] unused_nskin_friction_elements = line[4:7] #print(f'icomp={icomp} isk={isk} nskin_friction_elements={nskin_friction_elements}') i += 6 # skin friction method cards; A-74 # 0=level 2 viscous # 1=level 1 viscous #isfmeth = # 0=dont print # 1=print detailed skin friction intermediate results #iprint = # 0=don't save # 1=write skin friction coefficient for each element #isave = line = lines[i].rstrip() #print(line) return i
[docs] def parse_special_routines(lines, line, i, read_special_routines=False): """ Parses the special routines (method=5) 10000000000000000000 01100 1 3 0 0 0 0 01100 2 4 0 0 0 0 11100 3 1 2 5 6 7 """ #if not read_special_routines: #raise RuntimeError('special routines') line = lines[i].rstrip() #print('***', line) special_lines = [] while line[0:1] in ['0', '1']: special_lines.append(line) i += 1 line = lines[i].rstrip() #print('***', line) #print('**2', lines[i]) return i