Source code for upsp.processing.grids
import math
import numpy as np
###############################################################################
[docs]class StructGrid:
"""Manage plot3d-style structured grid and write formatted to file
Attributes
----------
sz : array_like
sizes of each grid [3xN]
x : array_like
x-position [N]
y : array_like,
y-position [N]
z : array_like,
z-position [N]
zones : np.ndarray
vertices ordered by zones [N]
"""
def __init__(self):
self.sz = []
self.x = []
self.y = []
self.z = []
self.zones = np.array([], dtype=np.float32)
[docs] def load_grid(self, grid_file):
"""Read a formatted p3d file
Parameters
----------
grid_file : str
formatted plot3d file
"""
with open(grid_file, "r") as f:
n_zones = int(f.readline())
for i in range(0, n_zones):
zone_sz = [int(x) for x in f.readline().split()]
self.sz.append(zone_sz)
zone_list = []
for z in range(0, n_zones):
total_sz = self.sz[z][0] * self.sz[z][1] * self.sz[z][2]
vals = []
curr_sz = 0
for line in f:
# Note the vertex-zone mapping
zone_list.append(z)
# Load all XYZ values
new_vals = [float(x) for x in line.split()]
vals.extend(new_vals)
curr_sz += len(new_vals)
if curr_sz >= 3 * total_sz:
break
self.x.extend(vals[:total_sz])
self.y.extend(vals[total_sz : 2 * total_sz])
self.z.extend(vals[2 * total_sz :])
self.zones = np.array(zone_list, dtype=np.float32)
[docs] def num_zones(self):
"""Return the number of grids (or zones)"""
return len(self.sz)
[docs] def size(self):
"""Return the number of grid nodes"""
total_size = 0
for i in range(0, self.num_zones()):
total_size += np.product(self.sz[i])
return total_size
[docs] def num_faces(self, zone=None):
"""Return the number of faces in a zone
Parameters
----------
zone : int, optional
zone number (0-based)
Returns
-------
n_faces : int
Number of faces in the given `zone` (if provided), otherwise the total
number of faces in the grid.
Raises
------
RuntimeError: invalid zone number
"""
n_faces = 0
if zone is None:
for z in range(0, self.num_zones()):
n_faces += self.num_faces(zone=z)
else:
# Make sure it is a real zone
if zone < 0 or zone >= len(self.sz):
raise RuntimeError("Invalid zone number")
# Get the size (handle surface grids in j,k,or l)
n_faces = (
max(1, self.sz[zone][0] - 1)
* max(1, self.sz[zone][1] - 1)
* max(1, self.sz[zone][2] - 1)
)
return n_faces
[docs] def write_p3d(self, fileout):
"""Write formatted p3d file
Parameters
----------
fileout : str
output file
"""
with open(fileout, "w") as f:
f.write("{}\n".format(len(self.sz)))
for z in range(0, len(self.sz)):
f.write("{} {} {}\n".format(*self.sz[z]))
idx = 0
for z in range(0, len(self.sz)):
num_pts = self.sz[z][0] * self.sz[z][1] * self.sz[z][2]
for i in range(0, num_pts):
f.write("{:7.4e} ".format(self.x[idx + i]))
for i in range(0, num_pts):
f.write("{:7.4e} ".format(self.y[idx + i]))
for i in range(0, num_pts):
f.write("{:7.4e} ".format(self.z[idx + i]))
f.write("\n")
idx += num_pts
[docs] def write_zones_mapping(self, fileout):
"""Write out the binary vertex zones mapping from a plot3d grid
Parameters
----------
fileout : str
output file
"""
self.zones.tofile(fileout)
[docs]class UnstructGrid:
"""Manages triangulated unstructured grid
Attributes
----------
n_comps : int
number of components
tris :array_like
node ids in each triangle [3,T]
comps : array_like
component id for each node [N]
x : array_like
x-position of each node [N]
y : array_like
y-position of each node [N]
z : array_like
z-position of each node [N]
"""
def __init__(self):
self.x = []
self.y = []
self.z = []
self.tris = [] # 2d array of sets of indices
self.comps = [] # len(comps) == len(tri)
self.n_comps = 0
[docs] def num_comps(self):
"""Return the number of components"""
return self.n_comps
[docs] def num_nodes(self):
"""Return the number of nodes"""
return len(self.x)
[docs] def num_faces(self, comp=None):
"""Return the number of faces in a component
Parameters
----------
comp : int, optional
component number (id)
Returns
-------
n_faces : int
Number of faces in the given `comp` (if provided), otherwise the total
number of faces in the grid.
"""
n_faces = 0
if comp is None:
n_faces = len(self.tris)
else:
for i in range(0, len(self.comps)):
if self.comps[i] == comp:
n_faces += 1
return n_faces
[docs] def get_area(self, t):
"""Return the area of a triangle
Parameters
----------
t : int
triangle index
Returns
-------
float
area of the triangle
Raises
------
RuntimeError : The triangle index is invalid
"""
if t < 0 or t >= self.num_faces():
raise RuntimeError("Invalid triangle, cannot compute area")
# Get points
p1 = np.array(
[self.x[self.tris[t][0]], self.y[self.tris[t][0]], self.z[self.tris[t][0]]]
)
p2 = np.array(
[self.x[self.tris[t][1]], self.y[self.tris[t][1]], self.z[self.tris[t][1]]]
)
p3 = np.array(
[self.x[self.tris[t][2]], self.y[self.tris[t][2]], self.z[self.tris[t][2]]]
)
# Use numerically stable Heron's formula
a = np.linalg.norm(p2 - p1)
b = np.linalg.norm(p3 - p2)
c = np.linalg.norm(p3 - p1)
if b > a:
tmp = a
a = b
b = tmp
if c > a:
tmp = a
a = c
c = b
elif c > b:
tmp = b
b = c
c = tmp
# could be negative due to small numerical error, so shortcut
pos_neg = abs(c - (a - b))
return 0.25 * math.sqrt((a + (b + c)) * pos_neg * (c + (a - b)) * (a + (b - c)))