"""
Support for calculating one or more scalar averages across one or more
regions in a domain.
"""
from math import cos, sin, sqrt
from openmdao.lib.datatypes.domain.flow import CELL_CENTER
from openmdao.lib.datatypes.domain.zone import CYLINDRICAL
from openmdao.lib.datatypes.domain.metrics import get_metric, list_metrics, \
create_scalar_metric
_SCHEMES = ('area', 'mass')
# TODO: account for ghost cells in index calculations.
[docs]def mesh_probe(domain, regions, variables, weighting_scheme='area'):
"""
Calculate metrics on mesh regions.
Currently only supports structured grids.
domain: DomainObj
The domain to be processed.
regions: list
List of ``(zone_name, imin, imax[, jmin, jmax[, kmin, kmax]])``
mesh region specifications to be used for the calculation.
Indices start at 0. Negative indices are relative to the end of the
array.
variables: list
List of ``(metric_name, units)`` tuples. Legal pre-existing metric
names can be obtained from :meth:`list_metrics`. If `units` is
None, then no unit conversion is attempted. `metric_name` may also
be the name of a flow solution scalar variable if `units` is None.
In this case a minimal metric calculator will be auto-generated.
weighting_scheme: string
Specifies how individual values are weighted. Legal values are
'area' for area averaging and 'mass' for mass averaging.
Returns a list of metric values in the order of the `variables` list.
.. note::
The per-item averaging scheme is simplistic. For instance, all four
vertices of a face get equal weight, or the two cells sharing a
face get equal weight. This will lead to inaccuracies for highly
irregular grids.
"""
# Check validity of region specifications.
_regions = _check_regions(domain, regions)
# Check validity of variables.
need_weights = False
for name, units in variables:
if name not in list_metrics():
if units is None:
# See if it's something we can create dynamically.
zone_name = regions[0][0]
zone = getattr(domain, zone_name)
if hasattr(zone.flow_solution, name):
create_scalar_metric(name)
else:
raise ValueError('Unknown variable %r' % name)
else:
raise ValueError('Unsupported variable %r' % name)
cls, integrate, geometry = get_metric(name)
if not integrate:
need_weights = True
# Check validity of weighting scheme.
if weighting_scheme not in _SCHEMES:
raise ValueError('Unknown/unsupported weighting scheme %r'
% weighting_scheme)
# Collect weights.
if need_weights:
weights, weight_total = _calc_weights(weighting_scheme, domain, _regions)
else:
weights, weight_total = {}, 0.
# Collect metric values.
metrics = []
for name, units in variables:
# Compute total for each region.
total = None
for region in _regions:
zone_name = region[0]
zone = getattr(domain, zone_name)
if units is None:
ref = None
else: # Check for a reference_state dictionary.
ref = zone.reference_state or domain.reference_state
if not ref:
raise ValueError('No zone or domain reference_state'
' dictionary supplied for zone %s.'
% zone_name)
value = _calc_metric(name, domain, region, weights, ref)
value *= zone.symmetry_instances # Adjust for symmetry.
if total is None:
total = value # Set initial PhysicalQuantity (or float).
else:
total += value
# If not integrating adjust for overall weighting.
cls, integrate, geometry = get_metric(name)
if not integrate:
total /= weight_total
if units is None:
metrics.append(total)
else:
total.convert_to_unit(units) # Convert to requested units.
metrics.append(total.value)
return metrics
def _check_regions(domain, regions):
""" Check validity of region specifications and normalize. """
_regions = []
for i, region in enumerate(regions, 1):
if len(region) == 7:
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
elif len(region) == 5:
zone_name, imin, imax, jmin, jmax = region
kmin, kmax = None, None
elif len(region) == 3:
zone_name, imin, imax = region
jmin, jmax, kmin, kmax = None, None, None, None
else:
raise ValueError('region specification %d invalid: %r'
% (i, region))
try:
zone = getattr(domain, zone_name)
except AttributeError:
raise ValueError('region %d: Domain does not contain zone %r'
% (i, zone_name))
shape = zone.shape
if len(zone.shape) > 2:
if len(region) != 7:
raise ValueError('region %d: index dimensionality mismatch' % i)
zone_imax, zone_jmax, zone_kmax = shape
elif len(zone.shape) > 1:
if len(region) != 5:
raise ValueError('region %d: index dimensionality mismatch' % i)
zone_imax, zone_jmax = zone.shape
zone_kmax = None
elif len(shape) > 0:
if len(region) != 3:
raise ValueError('region %d: index dimensionality mismatch' % i)
zone_imax, = zone.shape
zone_jmax, zone_kmax = None, None
else:
raise ValueError('region %d: zone is empty' % i)
# Normalize end-relative indexing.
if imin < 0:
imin = zone_imax + imin
if imax < 0:
imax = zone_imax + imax
if imin < 0 or imin >= zone_imax:
raise ValueError('region %d: imin %d invalid (max %d)'
% (i, imin, zone_imax))
if imax < imin or imax >= zone_imax:
raise ValueError('region %d: imax %d invalid (max %d)'
% (i, imax, zone_imax))
if zone_jmax is not None:
if jmin < 0:
jmin = zone_jmax + jmin
if jmax < 0:
jmax = zone_jmax + jmax
if jmin < 0 or jmin >= zone_jmax:
raise ValueError('region %d: jmin %d invalid (max %d)'
% (i, jmin, zone_jmax))
if jmax < jmin or jmax >= zone_jmax:
raise ValueError('region %d: jmax %d invalid (max %d)'
% (i, jmax, zone_jmax))
if zone_kmax is not None:
if kmin < 0:
kmin = zone_kmax + kmin
if kmax < 0:
kmax = zone_kmax + kmax
if kmin < 0 or kmin >= zone_kmax:
raise ValueError('region %d: kmin %d invalid (max %d)'
% (i, kmin, zone_kmax))
if kmax < kmin or kmax >= zone_kmax:
raise ValueError('region %d: kmax %d invalid (max %d)'
% (i, kmax, zone_kmax))
if len(region) == 7:
_regions.append((zone_name, imin, imax, jmin, jmax, kmin, kmax))
elif len(region) == 5:
_regions.append((zone_name, imin, imax, jmin, jmax))
else:
_regions.append((zone_name, imin, imax))
return _regions
def _get_dimension(region):
""" Return dimensionality of `region` (0, 1, 2, or 3). """
if len(region) == 7:
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
elif len(region) == 5:
zone_name, imin, imax, jmin, jmax = region
kmin, kmax = 0, 0
else:
zone_name, imin, imax = region
jmin, jmax, kmin, kmax = 0, 0, 0, 0
dim = 0
if imin != imax:
dim += 1
if jmin != jmax:
dim += 1
if kmin != kmax:
dim += 1
return dim
def _calc_weights(scheme, domain, regions):
"""
Calculate averaging weights, updating `weights` and returning total value.
"""
weights = {}
weight_total = 0.
for region in regions:
dim = _get_dimension(region)
if dim == 3:
zone_weights = _volume_weights(scheme, domain, region)
elif dim == 2:
if len(region) == 7:
zone_weights = _surface_weights_3d(scheme, domain, region)
else:
zone_weights = _surface_weights_2d(scheme, domain, region)
elif dim == 1:
if len(region) == 7:
zone_weights = _curve_weights_3d(scheme, domain, region)
elif len(region) == 5:
zone_weights = _curve_weights_2d(scheme, domain, region)
else:
zone_weights = _curve_weights_1d(scheme, domain, region)
else:
zone_weights = [1.]
zone_name = region[0]
zone = getattr(domain, zone_name)
zone_weights *= zone.symmetry_instances # Adjust for symmetry.
if zone_name in weights:
raise RuntimeError('Zone %r used more than once' % zone_name)
else:
weights[zone_name] = zone_weights
weight_total += sum(zone_weights)
return (weights, weight_total)
def _volume_weights(scheme, domain, region):
""" Returns weights for a mesh volume. """
raise NotImplementedError('_volume_weights')
def _surface_weights_3d(scheme, domain, region):
""" Returns 3D (index space) weights for a mesh surface. """
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
zone = getattr(domain, zone_name)
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = grid.z.item
if scheme == 'mass':
try:
if cylindrical:
mom_c1 = flow.momentum.z.item
mom_c2 = flow.momentum.r.item
mom_c3 = flow.momentum.t.item
else:
mom_c1 = flow.momentum.x.item
mom_c2 = flow.momentum.y.item
mom_c3 = flow.momentum.z.item
except AttributeError:
raise AttributeError("For mass averaging zone %s is missing"
" 'momentum'." % zone_name)
if imin == imax:
imax += 1
face_normal = _iface_normal
face_value = _iface_cell_value if cell_center else _iface_node_value
elif jmin == jmax:
jmax += 1
face_normal = _jface_normal
face_value = _jface_cell_value if cell_center else _jface_node_value
else:
kmax += 1
face_normal = _kface_normal
face_value = _kface_cell_value if cell_center else _kface_node_value
weights = []
for i in range(imin, imax):
for j in range(jmin, jmax):
for k in range(kmin, kmax):
sc1, sc2, sc3 = face_normal(c1, c2, c3, i, j, k, cylindrical)
if scheme == 'mass':
loc = (i, j, k)
rvu = face_value(mom_c1, loc)
rvv = face_value(mom_c2, loc)
rvw = face_value(mom_c3, loc)
weights.append(rvu*sc1 + rvv*sc2 + rvw*sc3)
else:
weights.append(sqrt(sc1*sc1 + sc2*sc2 + sc3*sc3))
return weights
def _surface_weights_2d(scheme, domain, region):
""" Returns 2D (index space) weights for a mesh surface. """
zone_name, imin, imax, jmin, jmax = region
zone = getattr(domain, zone_name)
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = None if grid.z is None else grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = None if grid.z is None else grid.z.item
if scheme == 'mass':
try:
if cylindrical:
mom_c1 = None if flow.momentum.z is None else flow.momentum.z.item
mom_c2 = flow.momentum.r.item
mom_c3 = flow.momentum.t.item
else:
mom_c1 = flow.momentum.x.item
mom_c2 = flow.momentum.y.item
mom_c3 = None if flow.momentum.z is None else flow.momentum.z.item
except AttributeError:
raise AttributeError("For mass averaging zone %s is missing"
" 'momentum'." % zone_name)
weights = []
for i in range(imin, imax):
for j in range(jmin, jmax):
sc1, sc2, sc3 = _cell_normal(c1, c2, c3, i, j, cylindrical)
if scheme == 'mass':
ip1 = i + 1
jp1 = j + 1
if cell_center:
# Cell value is value.
# FIXME: built-in ghosts
rvu = 0. if mom_c1 is None else mom_c1(ip1, jp1)
rvv = mom_c2(ip1, jp1)
rvw = 0. if mom_c1 is None else mom_c3(ip1, jp1)
else:
# Average across vertices.
if mom_c1 is None:
rvu = 0.
else:
rvu = 0.25 * (mom_c1(i, j) + mom_c1(ip1, j) + \
mom_c1(i, jp1) + mom_c1(ip1, jp1))
rvv = 0.25 * (mom_c2(i, j) + mom_c2(ip1, j) + \
mom_c2(i, jp1) + mom_c2(ip1, jp1))
if mom_c3 is None:
rvw = 0.
else:
rvw = 0.25 * (mom_c3(i, j) + mom_c3(ip1, j) + \
mom_c3(i, jp1) + mom_c3(ip1, jp1))
weights.append(rvu*sc1 + rvv*sc2 + rvw*sc3)
else:
weights.append(sqrt(sc1*sc1 + sc2*sc2 + sc3*sc3))
return weights
def _curve_weights_3d(scheme, domain, region):
""" Returns 3D (index space) weights for a mesh curve. """
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
zone = getattr(domain, zone_name)
grid = zone.grid_coordinates
cylindrical = zone.coordinate_system == CYLINDRICAL
if cylindrical:
raise NotImplementedError('curve weights for cylindrical coordinates')
else:
x = grid.x.item
y = grid.y.item
z = grid.z.item
if scheme == 'mass':
raise NotImplementedError('curve mass averaging')
along_i, along_j = False, False
if imin != imax:
along_i = True
jmax += 1
kmax += 1
elif jmin != jmax:
along_j = True
imax += 1
kmax += 1
else:
imax += 1
jmax += 1
weights = []
for i in range(imin, imax):
for j in range(jmin, jmax):
for k in range(kmin, kmax):
if along_i:
dx = x(i+1, j, k) - x(i, j, k)
dy = y(i+1, j, k) - y(i, j, k)
dz = z(i+1, j, k) - z(i, j, k)
elif along_j:
dx = x(i, j+1, k) - x(i, j, k)
dy = y(i, j+1, k) - y(i, j, k)
dz = z(i, j+1, k) - z(i, j, k)
else:
dx = x(i, j, k+1) - x(i, j, k)
dy = y(i, j, k+1) - y(i, j, k)
dz = z(i, j, k+1) - z(i, j, k)
weights.append(sqrt(dx*dx + dy*dy + dz*dz))
return weights
def _curve_weights_2d(scheme, domain, region):
""" Returns 2D (index space) weights for a mesh curve. """
zone_name, imin, imax, jmin, jmax = region
zone = getattr(domain, zone_name)
grid = zone.grid_coordinates
cylindrical = zone.coordinate_system == CYLINDRICAL
if cylindrical:
raise NotImplementedError('curve weights for cylindrical coordinates')
else:
x = grid.x.item
y = grid.y.item
z = None if grid.z is None else grid.z.item
if scheme == 'mass':
raise NotImplementedError('curve mass averaging')
if imin != imax:
along_i = True
jmax += 1
else:
along_i = False
imax += 1
weights = []
for i in range(imin, imax):
for j in range(jmin, jmax):
if along_i:
dx = x(i+1, j) - x(i, j)
dy = y(i+1, j) - y(i, j)
dz = 0. if z is None else z(i+1, j) - z(i, j)
else:
dx = x(i, j+1) - x(i, j)
dy = y(i, j+1) - y(i, j)
dz = 0. if z is None else z(i, j+1) - z(i, j)
weights.append(sqrt(dx*dx + dy*dy + dz*dz))
return weights
def _curve_weights_1d(scheme, domain, region):
""" Returns 1D (index space) weights for a mesh curve. """
zone_name, imin, imax = region
zone = getattr(domain, zone_name)
grid = zone.grid_coordinates
cylindrical = zone.coordinate_system == CYLINDRICAL
if cylindrical:
raise NotImplementedError('curve weights for cylindrical coordinates')
else:
x = grid.x.item
y = None if grid.y is None else grid.y.item
z = None if grid.z is None else grid.z.item
if scheme == 'mass':
raise NotImplementedError('curve mass averaging')
weights = []
for i in range(imin, imax):
dx = x(i+1) - x(i)
dy = 0. if y is None else y(i+1) - y(i)
dz = 0. if z is None else z(i+1) - z(i)
weights.append(sqrt(dx*dx + dy*dy + dz*dz))
return weights
def _calc_metric(name, domain, region, weights, reference_state):
"""
Calculate metric `name` on `region` using `weights` and `reference_state`.
"""
zone_name = region[0]
zone = getattr(domain, zone_name)
cls, integrate, geometry = get_metric(name)
metric = cls(zone, zone_name, reference_state)
weights = weights.get(zone_name)
# Could be volume, surface, curve, or point.
dim = _get_dimension(region)
if dim == 3:
if geometry not in ('volume', 'any'):
raise RuntimeError('metric %r not applicable to volumes')
total = _volume(metric, integrate, zone, region, weights)
elif dim == 2:
if geometry not in ('surface', 'any'):
raise RuntimeError('metric %r not applicable to surfaces')
if len(region) == 7:
total = _surface_3d(metric, integrate, zone, region, weights)
else:
total = _surface_2d(metric, integrate, zone, region, weights)
elif dim == 1:
if geometry not in ('curve', 'any'):
raise RuntimeError('metric %r not applicable to curves')
if len(region) == 7:
total = _curve_3d(metric, integrate, zone, region, weights)
elif len(region) == 5:
total = _curve_2d(metric, integrate, zone, region, weights)
else:
total = _curve_1d(metric, integrate, zone, region, weights)
else:
if geometry != 'any':
raise RuntimeError('metric %r not applicable to points')
total = _point(metric, zone, region)
if reference_state is None:
return total
else:
return metric.dimensionalize(total)
def _volume(metric, integrate, zone, region, weights):
""" Calculate metric on a volume. """
raise NotImplementedError('metric calculation on volume')
''' Not implemented yet
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = grid.z.item
vol = None
weight_index = 0
total = 0.
for i in range(imin, imax):
for j in range(jmin, jmax):
for k in range(kmin, kmax):
# TODO: calculate volume for integrating.
if cell_center:
# FIXME: built-in ghosts
# Cell value is value.
val = metric.calculate((i+1, j+1, k+1), vol)
else:
# Average across vertices.
val = metric.calculate((i, j, k), vol)
val += metric.calculate((i, j+1, k), vol)
val += metric.calculate((i, j+1, k+1), vol)
val += metric.calculate((i, j, k+1), vol)
val += metric.calculate((i+1, j, k), vol)
val += metric.calculate((i+1, j+1, k), vol)
val += metric.calculate((i+1, j+1, k+1), vol)
val += metric.calculate((i+1, j, k+1), vol)
val *= 0.125
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
'''
def _surface_3d(metric, integrate, zone, region, weights):
""" Calculate metric on a 3D (index space) surface. """
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = grid.z.item
if imin == imax:
face = 'i'
imax += 1
get_normal = _iface_normal
elif jmin == jmax:
face = 'j'
jmax += 1
get_normal = _jface_normal
else:
face = 'k'
kmax += 1
get_normal = _kface_normal
normal = None
weight_index = 0
total = 0.
for i in range(imin, imax):
for j in range(jmin, jmax):
for k in range(kmin, kmax):
if integrate:
normal = get_normal(c1, c2, c3, i, j, k, cylindrical)
if cell_center:
# FIXME: built-in ghosts
# Average across cells sharing surface.
val = metric.calculate((i+1, j+1, k+1), normal)
if face == 'i':
val += metric.calculate((i, j+1, k+1), normal)
elif face == 'j':
val += metric.calculate((i+1, j, k+1), normal)
else:
val += metric.calculate((i+1, j+1, k), normal)
val *= 0.5
else:
# Average across vertices.
val = metric.calculate((i, j, k), normal)
if face == 'i':
val += metric.calculate((i, j+1, k), normal)
val += metric.calculate((i, j+1, k+1), normal)
val += metric.calculate((i, j, k+1), normal)
elif face == 'j':
val += metric.calculate((i+1, j, k), normal)
val += metric.calculate((i+1, j, k+1), normal)
val += metric.calculate((i, j, k+1), normal)
else:
val += metric.calculate((i+1, j, k), normal)
val += metric.calculate((i+1, j+1, k), normal)
val += metric.calculate((i, j+1, k), normal)
val *= 0.25
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
def _surface_2d(metric, integrate, zone, region, weights):
""" Calculate metric on a 2D (index space) surface. """
zone_name, imin, imax, jmin, jmax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = None if grid.z is None else grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = None if grid.z is None else grid.z.item
normal = None
weight_index = 0
total = 0.
for i in range(imin, imax):
for j in range(jmin, jmax):
if integrate:
normal = _cell_normal(c1, c2, c3, i, j, cylindrical)
if cell_center:
# FIXME: built-in ghosts
# Cell value is value.
val = metric.calculate((i+1, j+1), normal)
else:
# Average across vertices.
val = metric.calculate((i, j), normal)
val += metric.calculate((i, j+1), normal)
val += metric.calculate((i+1, j+1), normal)
val += metric.calculate((i+1, j), normal)
val *= 0.25
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
def _curve_3d(metric, integrate, zone, region, weights):
""" Calculate metric on a 3D (index space) curve. """
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = grid.z.item
if imin != imax:
edge = 'i'
jmax += 1
kmax += 1
get_length = _iedge_length
elif jmin != jmax:
edge = 'j'
imax += 1
kmax += 1
get_length = _jedge_length
else:
edge = 'k'
imax += 1
jmax += 1
get_length = _kedge_length
length = None
weight_index = 0
total = 0.
for i in range(imin, imax):
for j in range(jmin, jmax):
for k in range(kmin, kmax):
if integrate:
length = get_length(c1, c2, c3, (i, j, k), cylindrical)
if cell_center:
# FIXME: built-in ghosts
# Average across cells sharing edge.
val = metric.calculate((i+1, j+1, k+1), length)
if edge == 'i':
val += metric.calculate((i+1, j, k+1), length)
val += metric.calculate((i+1, j+1, k), length)
val += metric.calculate((i+1, j, k), length)
elif edge == 'j':
val += metric.calculate((i, j+1, k+1), length)
val += metric.calculate((i+1, j+1, k), length)
val += metric.calculate((i, j+1, k), length)
else:
val += metric.calculate((i, j+1, k+1), length)
val += metric.calculate((i+1, j, k+1), length)
val += metric.calculate((i, j, k+1), length)
val *= 0.25
else:
# Average across vertices.
val = metric.calculate((i, j, k), length)
if edge == 'i':
val += metric.calculate((i+1, j, k), length)
elif edge == 'j':
val +=metric.calculate((i, j+1, k), length)
else:
val +=metric.calculate((i, j, k+1), length)
val *= 0.5
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
def _curve_2d(metric, integrate, zone, region, weights):
""" Calculate metric on a 2D (index space) curve. """
zone_name, imin, imax, jmin, jmax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = None if grid.z is None else grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = grid.y.item
c3 = None if grid.z is None else grid.z.item
if imin != imax:
edge = 'i'
jmax += 1
get_length = _iedge_length
else:
edge = 'j'
imax += 1
get_length = _jedge_length
length = None
weight_index = 0
total = 0.
for i in range(imin, imax):
for j in range(jmin, jmax):
if integrate:
length = get_length(c1, c2, c3, (i, j), cylindrical)
if cell_center:
# FIXME: built-in ghosts
# Average across cells sharing edge.
val = metric.calculate((i+1, j+1), length)
if edge == 'i':
val += metric.calculate((i+1, j), length)
else:
val += metric.calculate((i, j+1), length)
val *= 0.5
else:
# Average across vertices.
val = metric.calculate((i, j), length)
if edge == 'i':
val += metric.calculate((i+1, j), length)
else:
val += metric.calculate((i, j+1), length)
val *= 0.5
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
def _curve_1d(metric, integrate, zone, region, weights):
""" Calculate metric on a 1D (index space) curve. """
zone_name, imin, imax = region
grid = zone.grid_coordinates
flow = zone.flow_solution
cylindrical = zone.coordinate_system == CYLINDRICAL
cell_center = flow.grid_location == CELL_CENTER
if cylindrical:
c1 = None if grid.z is None else grid.z.item
c2 = grid.r.item
c3 = grid.t.item
else:
c1 = grid.x.item
c2 = None if grid.y is None else grid.y.item
c3 = None if grid.z is None else grid.z.item
get_length = _iedge_length
length = None
weight_index = 0
total = 0.
for i in range(imin, imax):
if integrate:
length = get_length(c1, c2, c3, (i,), cylindrical)
if cell_center:
# FIXME: built-in ghosts
# Cell value is value.
val = metric.calculate((i+1,), length)
else:
# Average across vertices.
val = metric.calculate((i,), length)
val += metric.calculate((i+1,), length)
val *= 0.5
if integrate:
total += val
else:
weight = weights[weight_index]
weight_index += 1
total += val * weight
return total
def _point(metric, zone, region):
""" Calculate metric at a point. """
zone_name = region[0]
flow = zone.flow_solution
cell_center = flow.grid_location == CELL_CENTER
if cell_center:
# Average across cells sharing point.
# FIXME: built-in ghosts
if len(region) == 7:
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
val = metric.calculate((imin+1, jmin+1, kmin+1), None)
val += metric.calculate((imin+1, jmin+1, kmin), None)
val += metric.calculate((imin+1, jmin, kmin), None)
val += metric.calculate((imin+1, jmin, kmin+1), None)
val += metric.calculate((imin, jmin+1, kmin+1), None)
val += metric.calculate((imin, jmin+1, kmin), None)
val += metric.calculate((imin, jmin, kmin), None)
val += metric.calculate((imin, jmin, kmin+1), None)
return 0.125 * val
elif len(region) == 5:
zone_name, imin, imax, jmin, jmax = region
val = metric.calculate((imin+1, jmin+1), None)
val += metric.calculate((imin, jmin+1), None)
val += metric.calculate((imin, jmin), None)
val += metric.calculate((imin+1, jmin), None)
return 0.25 * val
else:
zone_name, imin, imax = region
val = metric.calculate((imin+1,), None)
val += metric.calculate((imin,), None)
return 0.5 * val
else:
# Vertex value is value.
if len(region) == 7:
zone_name, imin, imax, jmin, jmax, kmin, kmax = region
return metric.calculate((imin, jmin, kmin), None)
elif len(region) == 5:
zone_name, imin, imax, jmin, jmax = region
return metric.calculate((imin, jmin), None)
else:
zone_name, imin, imax = region
return metric.calculate((imin,), None)
def _iface_normal(c1, c2, c3, i, j, k, cylindrical):
"""
Return non-dimensional vector normal to I face with magnitude equal to area.
"""
# FIXME: built-in ghosts
jp1 = j + 1
kp1 = k + 1
# upper-left - lower-right.
diag_c11 = c1(i, jp1, k) - c1(i, j, kp1)
diag_c21 = c2(i, jp1, k) - c2(i, j, kp1)
diag_c31 = c3(i, jp1, k) - c3(i, j, kp1)
# upper-right - lower-left.
diag_c12 = c1(i, jp1, kp1) - c1(i, j, k)
diag_c22 = c2(i, jp1, kp1) - c2(i, j, k)
diag_c32 = c3(i, jp1, kp1) - c3(i, j, k)
if cylindrical:
r1 = (c2(i, j, kp1) + c2(i, jp1, k )) / 2.
r2 = (c2(i, j, k ) + c2(i, jp1, kp1)) / 2.
else:
r1 = 1.
r2 = 1.
sc1 = -0.5 * ( r2 * diag_c21 * diag_c32 - r1 * diag_c22 * diag_c31)
sc2 = -0.5 * (-r2 * diag_c11 * diag_c32 + r1 * diag_c12 * diag_c31)
sc3 = -0.5 * ( diag_c11 * diag_c22 - diag_c12 * diag_c21)
return (sc1, sc2, sc3)
def _jface_normal(c1, c2, c3, i, j, k, cylindrical):
"""
Return non-dimensional vector normal to J face with magnitude equal to area.
"""
# FIXME: built-in ghosts
ip1 = i + 1
kp1 = k + 1
# upper-left - lower-right.
diag_c11 = c1(ip1, j, k) - c1(i, j, kp1)
diag_c21 = c2(ip1, j, k) - c2(i, j, kp1)
diag_c31 = c3(ip1, j, k) - c3(i, j, kp1)
# upper-right - lower-left.
diag_c12 = c1(ip1, j, kp1) - c1(i, j, k)
diag_c22 = c2(ip1, j, kp1) - c2(i, j, k)
diag_c32 = c3(ip1, j, kp1) - c3(i, j, k)
if cylindrical:
r1 = (c2(i, j, kp1) + c2(ip1, j, k )) / 2.
r2 = (c2(i, j, k ) + c2(ip1, j, kp1)) / 2.
else:
r1 = 1.
r2 = 1.
sc1 = 0.5 * ( r2 * diag_c21 * diag_c32 - r1 * diag_c22 * diag_c31)
sc2 = 0.5 * (-r2 * diag_c11 * diag_c32 + r1 * diag_c12 * diag_c31)
sc3 = 0.5 * ( diag_c11 * diag_c22 - diag_c12 * diag_c21)
return (sc1, sc2, sc3)
def _kface_normal(c1, c2, c3, i, j, k, cylindrical):
"""
Return non-dimensional vector normal to K face with magnitude equal to area.
"""
# FIXME: built-in ghosts
ip1 = i + 1
jp1 = j + 1
# upper-left - lower-right.
diag_c11 = c1(i, jp1, k) - c1(ip1, j, k)
diag_c21 = c2(i, jp1, k) - c2(ip1, j, k)
diag_c31 = c3(i, jp1, k) - c3(ip1, j, k)
# upper-right - lower-left.
diag_c12 = c1(ip1, jp1, k) - c1(i, j, k)
diag_c22 = c2(ip1, jp1, k) - c2(i, j, k)
diag_c32 = c3(ip1, jp1, k) - c3(i, j, k)
if cylindrical:
r1 = (c2(i, jp1, k) + c2(ip1, j, k)) / 2.
r2 = (c2(i, j, k) + c2(ip1, jp1, k)) / 2.
else:
r1 = 1.
r2 = 1.
sc1 = 0.5 * ( r2 * diag_c21 * diag_c32 - r1 * diag_c22 * diag_c31)
sc2 = 0.5 * (-r2 * diag_c11 * diag_c32 + r1 * diag_c12 * diag_c31)
sc3 = 0.5 * ( diag_c11 * diag_c22 - diag_c12 * diag_c21)
return (sc1, sc2, sc3)
def _cell_normal(c1, c2, c3, i, j, cylindrical):
"""
Return non-dimensional vector normal with magnitude equal to area.
If there is no 'z' coordinate, `c1` will be None in cylindrical
coordinates, otherwise `c3` will be None.
"""
# FIXME: built-in ghosts
ip1 = i + 1
jp1 = j + 1
# upper-left - lower-right.
diag_c11 = 0. if c1 is None else c1(i, jp1) - c1(ip1, j)
diag_c21 = c2(i, jp1) - c2(ip1, j)
diag_c31 = 0. if c3 is None else c3(i, jp1) - c3(ip1, j)
# upper-right - lower-left.
diag_c12 = 0. if c1 is None else c1(ip1, jp1) - c1(i, j)
diag_c22 = c2(ip1, jp1) - c2(i, j)
diag_c32 = 0. if c3 is None else c3(ip1, jp1) - c3(i, j)
if cylindrical:
r1 = (c2(i, jp1) + c2(ip1, j)) / 2.
r2 = (c2(i, j) + c2(ip1, jp1)) / 2.
else:
r1 = 1.
r2 = 1.
sc1 = 0.5 * ( r2 * diag_c21 * diag_c32 - r1 * diag_c22 * diag_c31)
sc2 = 0.5 * (-r2 * diag_c11 * diag_c32 + r1 * diag_c12 * diag_c31)
sc3 = 0.5 * ( diag_c11 * diag_c22 - diag_c12 * diag_c21)
return (sc1, sc2, sc3)
def _iedge_length(c1, c2, c3, loc, cylindrical):
""" Return length of edge along 'i'. """
if cylindrical:
if len(loc) > 2:
i, j, k = loc
ip1 = i + 1
theta = c3(ip1, j, k) - c3(i, j, k)
dx = c2(ip1, j, k) * cos(theta) - c2(i, j, k)
dy = c2(ip1, j, k) * sin(theta)
dz = c1(ip1, j, k) - c1(i, j, k)
elif len(loc) > 1:
i, j = loc
ip1 = i + 1
theta = c3(ip1, j) - c3(i, j)
dx = c2(ip1, j) * cos(theta) - c2(i, j)
dy = c2(ip1, j) * sin(theta)
dz = 0. if c1 is None else c1(ip1, j) - c1(i, j)
else:
i, = loc
ip1 = i + 1
theta = c3(ip1) - c3(i)
dx = c2(ip1) * cos(theta) - c2(i)
dy = c2(ip1) * sin(theta)
dz = 0. if c1 is None else c1(ip1) - c1(i)
else:
if len(loc) > 2:
i, j, k = loc
ip1 = i + 1
dx = c1(ip1, j, k) - c1(i, j, k)
dy = c2(ip1, j, k) - c2(i, j, k)
dz = c3(ip1, j, k) - c3(i, j, k)
elif len(loc) > 1:
i, j = loc
ip1 = i + 1
dx = c1(ip1, j) - c1(i, j)
dy = c2(ip1, j) - c2(i, j)
dz = 0. if c3 is None else c3(ip1, j) - c3(i, j)
else:
i, = loc
ip1 = i + 1
dx = c1(ip1) - c1(i)
dy = 0. if c2 is None else c2(ip1) - c2(i)
dz = 0. if c3 is None else c3(ip1) - c3(i)
return sqrt(dx*dx + dy*dy + dz*dz)
def _jedge_length(c1, c2, c3, loc, cylindrical):
""" Return length of edge along 'j'. """
if cylindrical:
if len(loc) > 2:
i, j, k = loc
jp1 = j + 1
theta = c3(i, jp1, k) - c3(i, j, k)
dx = c2(i, jp1, k) * cos(theta) - c2(i, j, k)
dy = c2(i, jp1, k) * sin(theta)
dz = c1(i, jp1, k) - c1(i, j, k)
else:
i, j = loc
jp1 = j + 1
theta = c3(i, jp1) - c3(i, j)
dx = c2(i, jp1) * cos(theta) - c2(i, j)
dy = c2(i, jp1) * sin(theta)
dz = 0. if c1 is None else c1(i, jp1) - c1(i, j)
else:
if len(loc) > 2:
i, j, k = loc
jp1 = j + 1
dx = c1(i, jp1, k) - c1(i, j, k)
dy = c2(i, jp1, k) - c2(i, j, k)
dz = c3(i, jp1, k) - c3(i, j, k)
else:
i, j = loc
jp1 = j + 1
dx = c1(i, jp1) - c1(i, j)
dy = c2(i, jp1) - c2(i, j)
dz = 0. if c3 is None else c3(i, jp1) - c3(i, j)
return sqrt(dx*dx + dy*dy + dz*dz)
def _kedge_length(c1, c2, c3, loc, cylindrical):
""" Return length of edge along 'k'. """
i, j, k = loc
kp1 = k + 1
if cylindrical:
theta = c3(i, j, kp1) - c3(i, j, k)
dx = c2(i, j, kp1) * cos(theta) - c2(i, j, k)
dy = c2(i, j, kp1) * sin(theta)
dz = c1(i, j, kp1) - c1(i, j, k)
else:
dx = c1(i, j, kp1) - c1(i, j, k)
dy = c2(i, j, kp1) - c2(i, j, k)
dz = c3(i, j, kp1) - c3(i, j, k)
return sqrt(dx*dx + dy*dy + dz*dz)
def _iface_cell_value(arr, loc):
""" Returns I face value for cell-centered data. """
i, j, k = loc
# FIXME: built-in ghosts
return 0.5 * (arr(i+1, j+1, k+1) + arr(i, j+1, k+1))
def _jface_cell_value(arr, loc):
""" Returns J face value for cell-centered data. """
i, j, k = loc
# FIXME: built-in ghosts
return 0.5 * (arr(i+1, j+1, k+1) + arr(i+1, j, k+1))
def _kface_cell_value(arr, loc):
""" Returns K face value for cell-centered data. """
i, j, k = loc
# FIXME: built-in ghosts
return 0.5 * (arr(i+1, j+1, k+1) + arr(i+1, j+1, k))
def _iface_node_value(arr, loc):
""" Returns I face value for vertex data. """
i, j, k = loc
return 0.25 * (arr(i, j+1, k+1) + arr(i, j, k+1) +
arr(i, j+1, k) + arr(i, j, k))
def _jface_node_value(arr, loc):
""" Returns J face value for vertex data. """
i, j, k = loc
return 0.25 * (arr(i+1, j, k+1) + arr(i, j, k+1) +
arr(i+1, j, k) + arr(i, j, k))
def _kface_node_value(arr, loc):
""" Returns K face value for vertex data. """
i, j, k = loc
return 0.25 * (arr(i+1, j+1, k) + arr(i, j+1, k) +
arr(i+1, j, k) + arr(i, j, k))
