"""
Defines:
- Cart3D(log=None, debug=False)
- read_cart3d(self, infilename, result_names=None)
- write_cart3d(self, outfilename, is_binary=False, float_fmt='%6.7f')
- flip_model()
- make_mirror_model(self, nodes, elements, regions, loads, axis='y', tol=0.000001)
- make_half_model(self, axis='y', remap_nodes=True)
- get_free_edges(self, elements)
- get_area(self)
- get_normals(self)
- get_normals_at_nodes(self, cnormals)
- comp2tri(in_filenames, out_filename,
is_binary=False, float_fmt='%6.7f')
"""
import sys
from struct import pack, unpack
from math import ceil
from collections import defaultdict
from typing import Tuple, Union
import numpy as np
from cpylog import get_logger2
from pyNastran.utils import is_binary_file, _filename, b
[docs]def read_cart3d(cart3d_filename, log=None, debug=False, result_names=None):
"""loads a Cart3D file"""
model = Cart3D(log=log, debug=debug)
model.read_cart3d(cart3d_filename, result_names)
return model
[docs]def comp2tri(in_filenames, out_filename,
is_binary=False, float_fmt='%6.7f',
log=None, debug=False):
"""
Combines multiple Cart3d files (binary or ascii) into a single file.
Parameters
----------
in_filenames : List[str]
list of filenames
out_filename : str
output filename
is_binary : bool; default=False
is the output file binary
float_fmt : str; default='%6.7f'
the format string to use for ascii writing
Notes
-----
assumes loads is None
"""
points = []
elements = []
regions = []
npoints = 0
nregions = 0
model = Cart3D(log=log, debug=debug)
for infilename in in_filenames:
model.read_cart3d(infilename)
npointsi = model.nodes.shape[0]
nregionsi = len(np.unique(model.regions))
points.append(model.nodes)
elements.append(model.elements + npoints)
regions.append(model.regions + nregions)
npoints += npointsi
nregions += nregionsi
points = np.vstack(points)
elements = np.vstack(elements)
regions = np.vstack(regions)
model.points = points
model.elements = elements
model.regions = regions
if out_filename:
model.write_cart3d(out_filename, is_binary=is_binary, float_fmt=float_fmt)
return model
[docs]class Cart3dIO:
"""
Cart3d IO class
"""
def __init__(self, log=None, debug=False):
self.log = get_logger2(log, debug=debug)
self._endian = b''
self._encoding = 'latin1'
self.n = 0
self.infile = None
self.infilename = None
self.points = None
self.elements = None
self.regions = None
self.loads = {}
def _write_header(self, outfile, points, elements, is_loads, is_binary=False):
"""
writes the cart3d header
Without results
---------------
npoints nelements
With results
------------
npoints nelements nresults
"""
npoints = points.shape[0]
nelements = elements.shape[0]
if is_binary:
if is_loads:
fmt = self._endian + b('iiiii')
msg = pack(fmt, 3*4, npoints, nelements, 6, 4)
else:
fmt = self._endian + b('iiii')
msg = pack(fmt, 2*4, npoints, nelements, 4)
int_fmt = None
else:
# this is ASCII data
if is_loads:
msg = b("%i %i 6\n" % (npoints, nelements))
else:
msg = b("%i %i\n" % (npoints, nelements))
# take the max value, string it, and length it
# so 123,456 is length 6
int_fmt = b('%%%si' % len(str(nelements)))
outfile.write(msg)
return int_fmt
def _write_points(self, outfile, points, is_binary, float_fmt='%6.6f'):
"""writes the points"""
if is_binary:
four = pack(self._endian + b('i'), 4)
outfile.write(four)
npoints = points.shape[0]
fmt = self._endian + b('%if' % (npoints * 3))
floats = pack(fmt, *np.ravel(points))
outfile.write(floats)
outfile.write(four)
else:
if isinstance(float_fmt, bytes):
fmt_ascii = float_fmt
else:
fmt_ascii = float_fmt.encode('latin1')
np.savetxt(outfile, points, fmt_ascii)
def _write_elements(self, outfile, elements, is_binary, int_fmt='%6i'):
"""writes the triangles"""
min_e = elements.min()
assert min_e == 0, 'min(elements)=%s' % min_e
if is_binary:
fmt = self._endian + b('i')
four = pack(fmt, 4)
outfile.write(four)
nelements = elements.shape[0]
fmt = self._endian + b('%ii' % (nelements * 3))
ints = pack(fmt, *np.ravel(elements+1))
outfile.write(ints)
outfile.write(four)
else:
if isinstance(int_fmt, bytes):
fmt_ascii = int_fmt
else:
fmt_ascii = int_fmt.encode('latin1')
np.savetxt(outfile, elements+1, fmt_ascii)
def _write_regions(self, outfile, regions, is_binary):
"""writes the regions"""
if is_binary:
fmt = self._endian + b('i')
four = pack(fmt, 4)
outfile.write(four)
nregions = len(regions)
fmt = self._endian + b('%ii' % nregions)
ints = pack(fmt, *regions)
outfile.write(ints)
outfile.write(four)
else:
fmt = b'%i'
np.savetxt(outfile, regions, fmt)
def _write_loads(self, outfile, loads, is_binary, float_fmt='%6.6f'):
"""writes the *.triq loads"""
if is_binary:
raise NotImplementedError('is_binary=%s' % is_binary)
else:
Cp = loads['Cp']
rho = loads['rho']
rhoU = loads['rhoU']
rhoV = loads['rhoV']
rhoW = loads['rhoW']
E = loads['E']
npoints = self.points.shape[0]
assert len(Cp) == npoints, 'len(Cp)=%s npoints=%s' % (len(Cp), npoints)
#nrows = len(Cp)
fmt = '%s\n%s %s %s %s %s\n' % (float_fmt, float_fmt, float_fmt,
float_fmt, float_fmt, float_fmt)
for (cpi, rhoi, rhou, rhov, rhoe, e) in zip(Cp, rho, rhoU, rhoV, rhoW, E):
outfile.write(fmt % (cpi, rhoi, rhou, rhov, rhoe, e))
def _read_header_ascii(self):
"""
Reads the header::
npoints nelements # geometry
npoints nelements nresults # results
"""
line = self.infile.readline()
sline = line.strip().split()
if len(sline) == 2:
npoints, nelements = int(sline[0]), int(sline[1])
nresults = 0
elif len(sline) == 3:
npoints = int(sline[0])
nelements = int(sline[1])
nresults = int(sline[2])
else:
raise ValueError('invalid result type')
return npoints, nelements, nresults
@property
def nresults(self):
"""get the number of results"""
if isinstance(self.loads, dict):
return len(self.loads)
return 0
@property
def nnodes(self):
"""alternate way to access number of points"""
return self.npoints
@property
def npoints(self):
"""get the number of points"""
return self.points.shape[0]
@property
def nodes(self):
"""alternate way to access the points"""
return self.points
@nodes.setter
def nodes(self, points):
"""alternate way to access the points"""
self.points = points
@property
def nelements(self):
"""get the number of elements"""
return self.elements.shape[0]
def _read_points_ascii(self, npoints):
"""
A point is defined by x,y,z and the ID is the location in points.
"""
p = 0
data = []
assert npoints > 0, 'npoints=%s' % npoints
points = np.zeros((npoints, 3), dtype='float32')
while p < npoints:
data += self.infile.readline().strip().split()
while len(data) > 2:
x = data.pop(0)
y = data.pop(0)
z = data.pop(0)
points[p] = [x, y, z]
p += 1
return points
def _read_elements_ascii(self, nelements):
"""
An element is defined by n1,n2,n3 and the ID is the location in elements.
"""
assert nelements > 0, 'npoints=%s nelements=%s' % (self.npoints, nelements)
elements = np.zeros((nelements, 3), dtype='int32')
ieid = 0
data = []
while ieid < nelements:
data += self.infile.readline().strip().split()
while len(data) > 2:
n1 = int(data.pop(0))
n2 = int(data.pop(0))
n3 = int(data.pop(0))
elements[ieid] = [n1, n2, n3]
ieid += 1
nid_min = elements.min()
if nid_min != 1:
nid_max = elements.max()
nnodes = self.nodes.shape[0]
if nid_max == nnodes:
msg = (
'Possible Cart3d error due to unused nodes\n'
'min(nids)=%s; expected 1; nid_max=%s nnodes=%s' % (
nid_min, nid_max, nnodes))
self.log.warning(msg)
else:
msg = 'elements:\n%s\nmin(nids)=%s; expected 1; nid_max=%s nnodes=%s' % (
elements, nid_min, nid_max, nnodes, )
raise RuntimeError(msg)
#assert elements.min() == 1, elements.min()
return elements - 1
def _read_regions_ascii(self, nelements):
"""reads the region section"""
regions = np.zeros(nelements, dtype='int32')
iregion = 0
data = []
while iregion < nelements:
data = self.infile.readline().strip().split()
ndata = len(data)
regions[iregion : iregion + ndata] = data
iregion += ndata
return regions
def _read_header_binary(self):
"""
Reads the header::
npoints nelements # geometry
npoints nelements nresults # results
"""
data = self.infile.read(4)
size_little, = unpack(b'<i', data)
size_big, = unpack(b'>i', data)
if size_big in [12, 8]:
self._endian = b'>'
size = size_big
elif size_little in [8, 12]:
self._endian = b'<'
size = size_little
else:
self._rewind()
self.show(100)
raise RuntimeError('unknown endian')
self.n += 4
data = self.infile.read(size)
self.n += size
so4 = size // 4 # size over 4
if so4 == 3:
(npoints, nelements, nresults) = unpack(self._endian + b('iii'), data)
self.log.info("npoints=%s nelements=%s nresults=%s" % (npoints, nelements, nresults))
elif so4 == 2:
(npoints, nelements) = unpack(self._endian + b('ii'), data)
nresults = 0
self.log.info("npoints=%s nelements=%s" % (npoints, nelements))
else:
self._rewind()
self.show(100)
raise RuntimeError('in the wrong spot...endian...size/4=%s' % so4)
self.infile.read(8) # end of first block, start of second block
return (npoints, nelements, nresults)
def _read_points_binary(self, npoints):
"""reads the xyz points"""
size = npoints * 12 # 12=3*4 all the points
data = self.infile.read(size)
dtype = np.dtype(self._endian + b('f4'))
points = np.frombuffer(data, dtype=dtype).reshape((npoints, 3)).copy()
self.infile.read(8) # end of second block, start of third block
return points
def _read_elements_binary(self, nelements):
"""reads the triangles"""
size = nelements * 12 # 12=3*4 all the elements
data = self.infile.read(size)
dtype = np.dtype(self._endian + b('i4'))
elements = np.frombuffer(data, dtype=dtype).reshape((nelements, 3)).copy()
self.infile.read(8) # end of third (element) block, start of regions (fourth) block
assert elements.min() == 1, elements.min()
return elements - 1
def _read_regions_binary(self, nelements):
"""reads the regions"""
size = nelements * 4 # 12=3*4 all the elements
data = self.infile.read(size)
regions = np.zeros(nelements, dtype='int32')
dtype = self._endian + b'i'
regions = np.frombuffer(data, dtype=dtype).copy()
self.infile.read(4) # end of regions (fourth) block
return regions
def _read_results_binary(self, i, infile, result_names=None):
"""binary results are not supported"""
pass
def _rewind(self): # pragma: no cover
"""go back to the beginning of the file"""
self.n = 0
self.infile.seek(self.n)
[docs] def show(self, n, types='ifs', endian=None): # pragma: no cover
assert self.n == self.infile.tell(), 'n=%s tell=%s' % (self.n, self.infile.tell())
#nints = n // 4
data = self.infile.read(4 * n)
strings, ints, floats = self.show_data(data, types=types, endian=endian)
self.infile.seek(self.n)
return strings, ints, floats
[docs] def show_data(self, data, types='ifs', endian=None): # pragma: no cover
return self._write_data(sys.stdout, data, types=types, endian=endian)
def _write_data(self, outfile, data, types='ifs', endian=None): # pragma: no cover
"""
Useful function for seeing what's going on locally when debugging.
"""
n = len(data)
nints = n // 4
ndoubles = n // 8
strings = None
ints = None
floats = None
longs = None
if endian is None:
endian = self._endian
if 's' in types:
strings = unpack(b'%s%is' % (endian, n), data)
outfile.write("strings = %s\n" % str(strings))
if 'i' in types:
ints = unpack(b'%s%ii' % (endian, nints), data)
outfile.write("ints = %s\n" % str(ints))
if 'f' in types:
floats = unpack(b'%s%if' % (endian, nints), data)
outfile.write("floats = %s\n" % str(floats))
if 'l' in types:
longs = unpack(b'%s%il' % (endian, nints), data)
outfile.write("long = %s\n" % str(longs))
if 'I' in types:
ints2 = unpack(b'%s%iI' % (endian, nints), data)
outfile.write("unsigned int = %s\n" % str(ints2))
if 'L' in types:
longs2 = unpack(b'%s%iL' % (endian, nints), data)
outfile.write("unsigned long = %s\n" % str(longs2))
if 'q' in types:
longs = unpack(b'%s%iq' % (endian, ndoubles), data[:ndoubles*8])
outfile.write("long long = %s\n" % str(longs))
return strings, ints, floats
[docs] def show_ndata(self, n, types='ifs'): # pragma: no cover
return self._write_ndata(sys.stdout, n, types=types)
def _write_ndata(self, outfile, n, types='ifs'): # pragma: no cover
"""
Useful function for seeing what's going on locally when debugging.
"""
nold = self.n
data = self.infile.read(n)
self.n = nold
self.infile.seek(self.n)
return self._write_data(outfile, data, types=types)
[docs]class Cart3D(Cart3dIO):
"""Cart3d interface class"""
model_type = 'cart3d'
is_structured = False
is_outward_normals = True
def __init__(self, log=None, debug=False):
Cart3dIO.__init__(self, log=log, debug=debug)
self.loads = {}
self.points = None
self.elements = None
[docs] def flip_model(self):
"""flip the model about the y-axis"""
self.points[:, 1] *= -1.
self.elements = np.hstack([
self.elements[:, 0:1],
self.elements[:, 2:3],
self.elements[:, 1:2],
])
print(self.elements.shape)
[docs] def make_mirror_model(self, nodes, elements, regions, loads, axis='y', tol=0.000001):
"""
Makes a full cart3d model about the given axis.
Parameters
----------
nodes : (nnodes, 3) ndarray
the nodes
elements : (nelements, 3) ndarray
the elmements
regions : (nelements) ndarray
the regions
loads : dict[str] = (nnodes) ndarray
not supported
axis : str; {"x", "y", "z", "-x", "-y", "-z"}
a string of the axis
tol : float; default=0.000001
the tolerance for the centerline points
"""
raise NotImplementedError()
#self.log.info('---starting make_mirror_model---')
#assert tol >= 0, 'tol=%r' % tol # prevents hacks to the axis
#nnodes = nodes.shape[0]
#assert nnodes > 0, 'nnodes=%s' % nnodes
#nelements = elements.shape[0]
#assert nelements > 0, 'nelements=%s' % nelements
#ax, ax0 = self._get_ax(axis)
#if ax in [0, 1, 2]: # positive x, y, z values; mirror to -side
#iy0 = np.where(nodes[:, ax] > tol)[0]
#ax2 = ax
#elif ax in [3, 4, 5]: # negative x, y, z values; mirror to +side
#iy0 = np.where(nodes[:, ax-3] < -tol)[0]
#ax2 = ax - 3 # we create ax2 in order to generalize the data updating
#else:
#raise NotImplementedError(axis)
## the nodes to be duplicated are the nodes that aren't below the tolerance
#nodes_upper = nodes[iy0]
#nodes_upper[:, ax2] *= -1.0 # flip the nodes about the axis
#nodes2 = np.vstack([nodes, nodes_upper])
#nnodes2 = nodes2.shape[0]
#assert nnodes2 > nnodes, 'nnodes2=%s nnodes=%s' % (nnodes2, nnodes)
#nnodes_upper = nodes_upper.shape[0]
#elements_upper = elements.copy()
#nelements = elements.shape[0]
## remap the mirrored nodes with the new node ids
#for eid in range(nelements):
#element = elements_upper[eid, :]
#for i, eidi in enumerate(element):
#if eidi in iy0:
#elements_upper[eid][i] = nnodes_upper + eidi
## we need to reverse the element in order to get
## the proper normal vector
#elements_upper[eid] = elements_upper[eid, ::-1]
#elements2 = np.vstack([elements, elements_upper])
#nelements2 = elements2.shape[0]
#assert nelements2 > nelements, 'nelements2=%s nelements=%s' % (nelements2, nelements)
#nregions = len(np.unique(regions))
#regions_upper = regions.copy() + nregions
#regions2 = np.hstack([regions, regions_upper])
#loads2 = {}
#for key, data in loads.items():
## flip the sign on the flipping axis terms
#if((key in ['U', 'rhoU'] and ax2 == 0) or
#(key in ['V', 'rhoV'] and ax2 == 1) or
#(key in ['W', 'rhoW'] and ax2 == 2)):
#data_upper = -data[iy0]
#else:
#data_upper = data[iy0]
#loads2[key] = np.hstack([data, data_upper])
#self.log.info('---finished make_mirror_model---')
#return (nodes2, elements2, regions2, loads2)
def _get_ax(self, axis: Union[str, int]) -> Tuple[int, int]:
"""helper method to convert an axis_string into an integer"""
if isinstance(axis, str):
axis = axis.lower().strip()
if axis in ['+x', 'x', 0]:
ax = 0
elif axis in ['+y', 'y', 1]:
ax = 1
elif axis in ['+z', 'z', 2]:
ax = 2
elif axis in ['-x', 3]:
ax = 3
elif axis == ['-y', 4]:
ax = 4
elif axis == ['-z', 5]:
ax = 5
else: # pragma: no cover
raise NotImplementedError('axis=%r' % axis)
self.log.debug("axis=%r ax=%s" % (axis, ax))
# shift ax to the actual column index in the nodes array
ax0 = ax if ax in [0, 1, 2] else ax - 3
return ax, ax0
[docs] def make_half_model(self, axis='y', remap_nodes=True):
"""
Makes a half model from a full model
Notes
-----
Cp is really loads['Cp'] and was meant for loads analysis only
"""
nodes = self.nodes
elements = self.elements
regions = self.regions
loads = self.loads
if loads is None:
loads = {}
nnodes = nodes.shape[0]
assert nnodes > 0, 'nnodes=%s' % nnodes
nelements = elements.shape[0]
assert nelements > 0, 'nelements=%s' % nelements
self.log.info('---starting make_half_model---')
ax, ax0 = self._get_ax(axis)
self.log.debug('find elements to remove')
#print(f'axis={axis} ax={ax}')
ynode = nodes[:, ax0]
y1 = ynode[elements[:, 0]]
y2 = ynode[elements[:, 1]]
y3 = ynode[elements[:, 2]]
if ax in [0, 1, 2]: # keep values > 0
ielements_save = np.where((y1 >= 0.0) & (y2 >= 0.0) & (y3 >= 0.0))[0]
elif ax in [3, 4, 5]: # keep values < 0
ielements_save = np.where((y1 <= 0.0) & (y2 <= 0.0) & (y3 <= 0.0))[0]
else:
raise NotImplementedError('axis=%r ax=%s' % (axis, ax))
elements2 = elements[ielements_save, :]
regions2 = regions[ielements_save]
nelements2 = elements2.shape[0]
assert 0 < nelements2 < nelements, 'nelements=%s nelements2=%s' % (nelements, nelements2)
# location of nodes to keep
inodes_save = np.unique(elements2.ravel())
#------------------------------------------
#checks
#inodes_map = np.arange(len(inodes_save))
is_nodes = 0 < len(inodes_save) < nnodes
if not is_nodes:
msg = 'There are no nodes in the model; len(inodes_save)=%s nnodes=%s' % (
len(inodes_save), nnodes)
raise RuntimeError(msg)
nodes2 = nodes[inodes_save, :]
nnodes2 = nodes2.shape[0]
assert 0 < nnodes2 < nnodes, 'no nodes were removed; nnodes=%s nnodes2=%s' % (nnodes, nnodes2)
#------------------------------------------
# renumbers mesh
# node id renumbering
emin0 = elements2.min()
emax0 = elements2.max()
# only renumber if we need to
if emin0 > 0 or emax0 != nelements2:
# we're making a 0-based node map of old -> new node_id
nodes_map = {}
for i, nid in enumerate(inodes_save):
nodes_map[nid] = i
# update the node ids
for ielement in range(nelements2):
element2 = elements2[ielement, :]
elements2[ielement, :] = [nodes_map[nid] for nid in element2]
assert len(elements2) > 0, elements2.shape
min_e = elements2.min()
assert min_e == 0, 'min(elements)=%s' % min_e
#------------------------------------------
loads2 = {} # 'Cp', 'Mach', 'U', etc.
for key, load in loads.items():
loads2[key] = load[inodes_save]
self.log.info('---finished make_half_model---')
self.nodes = nodes2
self.elements = elements2
self.regions = regions2
self.loads = loads2
return (nodes2, elements2, regions2, loads2)
[docs] def get_free_edges(self, elements):
"""
Cart3d must be a closed model with each edge shared by 2 elements
The free edges indicate the problematic areas.
Returns
-------
free edges : (nedges, 2) int ndarray
the free edge node ids
"""
edge_to_eid_map = defaultdict(list)
for i, element in enumerate(elements):
edge1 = tuple(sorted([element[0], element[1]]))
edge2 = tuple(sorted([element[1], element[2]]))
edge3 = tuple(sorted([element[2], element[0]]))
edge_to_eid_map[edge1].append(i)
edge_to_eid_map[edge2].append(i)
edge_to_eid_map[edge3].append(i)
free_edges = []
for edge, eids in sorted(edge_to_eid_map.items()):
if len(eids) != 2:
free_edges.append(edge)
return np.array(free_edges, dtype='int32')
[docs] def read_cart3d(self, infilename, result_names=None):
"""extracts the points, elements, and Cp"""
self.infilename = infilename
self.log.info("---reading cart3d...%r---" % self.infilename)
self.infilename = infilename
if is_binary_file(infilename):
with open(infilename, 'rb') as self.infile:
try:
npoints, nelements, nresults = self._read_header_binary()
self.points = self._read_points_binary(npoints)
self.elements = self._read_elements_binary(nelements)
self.regions = self._read_regions_binary(nelements)
# TODO: loads
except Exception:
msg = 'failed reading %r' % infilename
self.log.error(msg)
raise
else:
with open(_filename(infilename), 'r', encoding=self._encoding) as self.infile:
try:
npoints, nelements, nresults = self._read_header_ascii()
self.points = self._read_points_ascii(npoints)
self.elements = self._read_elements_ascii(nelements)
self.regions = self._read_regions_ascii(nelements)
self._read_results_ascii(0, self.infile, nresults, result_names=result_names)
except Exception:
msg = 'failed reading %r' % infilename
self.log.error(msg)
raise
self.log.debug("npoints=%s nelements=%s" % (self.npoints, self.nelements))
assert self.npoints > 0, 'npoints=%s' % self.npoints
assert self.nelements > 0, 'nelements=%s' % self.nelements
[docs] def write_cart3d(self, outfilename, is_binary=False, float_fmt='%6.7f'):
"""
writes a cart3d file
"""
assert len(self.points) > 0, 'len(self.points)=%s' % len(self.points)
if self.loads is None or self.loads == {}:
#loads = {}
is_loads = False
else:
is_loads = True
self.log.info("---writing cart3d...%r---" % outfilename)
with open(outfilename, 'wb') as outfile:
int_fmt = self._write_header(outfile, self.points, self.elements, is_loads, is_binary)
self._write_points(outfile, self.points, is_binary, float_fmt)
self._write_elements(outfile, self.elements, is_binary, int_fmt)
self._write_regions(outfile, self.regions, is_binary)
if is_loads:
assert is_binary is False, 'is_binary=%r is not supported for loads' % is_binary
self._write_loads(outfile, self.loads, is_binary, float_fmt)
def _read_results_ascii(self, i, infile, nresults, result_names=None):
"""
Reads the Cp results.
Results are read on a nodal basis from the following table:
Cp
rho,rhoU,rhoV,rhoW,rhoE
With the following definitions:
Cp = (p - 1/gamma) / (0.5*M_inf*M_inf)
rhoVel^2 = rhoU^2+rhoV^2+rhoW^2
M^2 = rhoVel^2/rho^2
Thus:
p = (gamma-1)*(e- (rhoU**2+rhoV**2+rhoW**2)/(2.*rho))
p_dimensional = qInf * Cp + pInf
# ???
rho,rhoU,rhoV,rhoW,rhoE
Parameters
----------
result_names : List[str]; default=None (All)
result_names = ['Cp', 'rho', 'rhoU', 'rhoV', 'rhoW', 'rhoE',
'Mach', 'U', 'V', 'W', 'E']
"""
if nresults == 0:
return
#loads = {}
if result_names is None:
result_names = ['Cp', 'rho', 'rhoU', 'rhoV', 'rhoW', 'rhoE',
'Mach', 'U', 'V', 'W', 'E', 'a', 'T', 'Pressure', 'q']
self.log.debug('---starting read_results---')
results = np.zeros((self.npoints, 6), dtype='float32')
nresult_lines = int(ceil(nresults / 5.)) - 1
for ipoint in range(self.npoints):
# rho rhoU,rhoV,rhoW,pressure/rhoE/E
sline = infile.readline().strip().split()
i += 1
for unused_n in range(nresult_lines):
sline += infile.readline().strip().split() # Cp
i += 1
#gamma = 1.4
#else:
# p=0.
sline = _get_list(sline)
# Cp
# rho rhoU rhoV rhoW E
# 0.416594
# 1.095611 0.435676 0.003920 0.011579 0.856058
results[ipoint, :] = sline
#p=0
#cp = sline[0]
#rho = float(sline[1])
#if(rho > abs(0.000001)):
#rhoU = float(sline[2])
#rhoV = float(sline[3])
#rhoW = float(sline[4])
#rhoE = float(sline[5])
#mach2 = (rhoU) ** 2 + (rhoV) ** 2 + (rhoW) ** 2 / rho ** 2
#mach = sqrt(mach2)
#if mach > 10:
#print("nid=%s Cp=%s mach=%s rho=%s rhoU=%s rhoV=%s rhoW=%s" % (
#pointNum, cp, mach, rho, rhoU, rhoV, rhoW))
#print("pt=%s i=%s Cp=%s p=%s" %(pointNum,i,sline[0],p))
del sline
self.loads = self._calculate_results(result_names, results)
def _calculate_results(self, result_names, results, loads=None):
"""
Takes the Cart3d variables and calculates additional variables
Parameters
----------
result_names : List[str]
the variables to calculate
results : (n,6) ndarray
the non-dimensional primitive flow variables
loads : dict; default=None -> {}
key : ???
value : ???
"""
if loads is None:
loads = {}
Cp = results[:, 0]
rho = results[:, 1]
rho_u = results[:, 2]
rho_v = results[:, 3]
rho_w = results[:, 4]
E = results[:, 5]
ibad = np.where(rho <= 0.000001)[0]
if len(ibad) > 0:
if 'Mach' in result_names:
Mach = np.sqrt(rho_u**2 + rho_v**2 + rho_w**2)# / rho
Mach[ibad] = 0.0
if 'U' in result_names:
U = rho_u / rho
U[ibad] = 0.0
if 'U' in result_names:
V = rho_v / rho
V[ibad] = 0.0
if 'W' in result_names:
W = rho_w / rho
W[ibad] = 0.0
#if 'rhoE' in result_names:
#rho_e = rhoE / rho
#e[ibad] = 0.0
is_bad = True
#n = 0
#for i in ibad:
#print("nid=%s Cp=%s mach=%s rho=%s rhoU=%s rhoV=%s rhoW=%s" % (
#i, Cp[i], Mach[i], rho[i], rho_u[i], rho_v[i], rho_w[i]))
#Mach[i] = 0.0
#n += 1
#if n > 10:
# break
else:
is_bad = False
#loc = locals()
if 'Cp' in result_names:
loads['Cp'] = Cp
if 'rhoU' in result_names:
loads['rhoU'] = rho_u
if 'rhoV' in result_names:
loads['rhoV'] = rho_v
if 'rhoW' in result_names:
loads['rhoW'] = rho_w
#if 'rhoE' in result_names:
#loads['rhoE'] = rho_e
if 'rho' in result_names:
loads['rho'] = rho
if 'Mach' in result_names:
if not is_bad:
#Mach = np.sqrt(rho_u**2 + rho_v**2 + rho_w**2) / rho
Mach = np.sqrt(rho_u**2 + rho_v**2 + rho_w**2)
loads['Mach'] = Mach
if 'U' in result_names:
if not is_bad:
U = rho_u / rho
loads['U'] = U
if 'V' in result_names:
if not is_bad:
V = rho_v / rho
loads['V'] = V
if 'W' in result_names:
if not is_bad:
W = rho_w / rho
loads['W'] = W
if 'E' in result_names:
#if not is_bad:
#E = rhoE / rho
loads['E'] = E
gamma = 1.4
qinf = 1.0
pinf = 1. / gamma
Tinf = 1.0
#Cp = (p - pinf) / qinf
p = Cp * qinf + pinf
T = (Tinf * gamma) * p / rho
q = 0.5 * rho * Mach ** 2
if 'a' in result_names:
#print('T: min=%s max=%s' % (T.min(), T.max()))
loads['a'] = np.sqrt(T)
if 'T' in result_names:
loads['T'] = T
if 'Pressure' in result_names:
loads['Pressure'] = p
if 'q' in result_names:
loads['q'] = q
# dynamic pressure
# speed of sound
# total pressure = p0/rhoi*ainf**2
# total density
# entropy
# kinetic energy
# enthalpy
# energy, E
# total energy
# total enthalpy
#i = where(Mach == max(Mach))[0][0]
#self.log.info("i=%s Cp=%s rho=%s rho_u=%s rho_v=%s rho_w=%s Mach=%s" % (
#i, Cp[i], rho[i], rho_u[i], rho_v[i], rho_w[i], Mach[i]))
self.log.debug('---finished read_results---')
return loads
def _get_area_vector(self):
"""
Gets the area vector (unnormalized normal vector)
Returns
-------
normals : (n, 3) ndarray
unnormalized centroidal normal vectors
"""
elements = self.elements
nodes = self.nodes
p1 = nodes[elements[:, 0], :]
p2 = nodes[elements[:, 1], :]
p3 = nodes[elements[:, 2], :]
ne = elements.shape[0]
avec = p2 - p1
bvec = p3 - p1
n = np.cross(avec, bvec)
assert len(n) == ne, 'len(n)=%s ne=%s' % (len(n), ne)
return n
[docs] def get_area(self):
"""
Gets the element area
Returns
-------
area : (n, 3) ndarray
the element areas
"""
ne = self.elements.shape[0]
n = self._get_area_vector()
ni = np.linalg.norm(n, axis=1)
assert len(ni) == ne, 'len(ni)=%s ne=%s' % (len(ni), ne)
return 0.5 * ni
[docs] def get_normals(self):
"""
Gets the centroidal normals
Returns
-------
cnormals : (n, 3) ndarray
normalized centroidal normal vectors
"""
ne = self.elements.shape[0]
n = self._get_area_vector()
ni = np.linalg.norm(n, axis=1)
assert len(ni) == ne, 'len(ni)=%s ne=%s' % (len(ni), ne)
assert ni.min() > 0.0, ni[np.where(ni <= 0.0)[0]]
n /= ni[:, None] # normal vector
return n
[docs] def get_normals_at_nodes(self, cnormals):
"""
Gets the nodal normals
Parameters
----------
cnormals : (n, 3) ndarray
normalized centroidal normal vectors
Returns
-------
nnormals : (n, 3) ndarray
normalized nodal normal vectors
"""
elements = self.elements
#nodes = self.nodes
nnodes = self.nnodes
nid_to_eids = defaultdict(list)
# find the elements to consider for each node
for eid, element in enumerate(elements):
n1, n2, n3 = element
nid_to_eids[n1].append(eid)
nid_to_eids[n2].append(eid)
nid_to_eids[n3].append(eid)
nnormals = np.zeros((nnodes, 3), dtype='float64')
for nid in range(nnodes):
eids = nid_to_eids[nid]
if len(eids) == 0:
raise RuntimeError('nid=%s is not used' % nid)
#ni_avg = cnormals[eids, :]
nnormals[nid] = cnormals[eids, :].sum(axis=0)
ni = np.linalg.norm(nnormals, axis=1)
assert ni.min() > 0, ni
nnormals /= ni[:, None] # normal vector
return nnormals
[docs]def convert_to_float(svalues):
"""Takes a list of strings and converts them to floats."""
values = []
for value in svalues:
values.append(float(value))
return values
def _get_list(sline):
"""Takes a list of strings and converts them to floats."""
try:
sline2 = convert_to_float(sline)
except ValueError:
print("sline = %s" % sline)
raise SyntaxError('cannot parse %s' % sline)
return sline2