"""
Define region
"""
__author__="steve"
__date__ ="$09/03/2012 4:46:39 PM$"
from anuga.utilities.numerical_tools import ensure_numeric
from anuga import Domain
from anuga import Quantity
import numpy as num
from pprint import pprint
from anuga.geometry.polygon import inside_polygon, line_intersect
from anuga.utilities.function_utils import determine_function_type
#from anuga import indent
[docs]class Region(object):
""" Object which defines a region within the domain
"""
[docs] def __init__(self,
domain,
indices=None,
polygon=None,
center=None,
radius=None,
line=None,
poly=None,
expand_polygon=False,
verbose = False):
"""Create a Region object
:param domain: Region must be defined wrt a domain
:param indices: Define the region by triangle IDs
:param polygon: List of [x,y] points to define region
:param center: point [x,y] which defines the centre of a circle
:param radius: radius of a circle which defines a region
:param line: List of [x,y] points defining a polyline
:param poly: An old argument which was used to define a polyline or polygon
:param expand_polygon: If set true, then calculation of intersection of polygon with triangles based on vertices, otherwise based just on centroids
:param verbose: Set to True for more verbose output
Setup region (defined by indices, polygon or center/radius).
Useful in defining where to apply certain operations
"""
#------------------------------------------
# Local variables
#------------------------------------------
self.indices = indices
self.center = center
self.radius = radius
self.polygon = polygon
self.line = line
self.poly = poly
self.type = 'full_region'
self.expand_polygon = expand_polygon
self.verbose = verbose
#------------------------------------------
# Useful aliases
#------------------------------------------
self.domain = domain
self.coord_c = self.domain.centroid_coordinates
self.areas = self.domain.areas
#-------------------------------------------
# Work out indices associated with region:
#-------------------------------------------
if self.indices is not None:
# This overrides polygon, center and radius, line
assert self.radius is None
assert self.center is None
assert self.polygon is None
assert self.line is None
self.indices = num.asarray(self.indices)
if self.indices.size == 0:
self.indices = []
self.type = 'empty'
else:
self.type = 'indices_specified'
elif (self.center is not None) and (self.radius is not None):
assert self.indices is None
assert self.polygon is None
assert self.line is None
self._setup_indices_circle()
self.type = 'circle'
elif (self.polygon is not None):
assert self.indices is None
assert self.radius is None
assert self.center is None
assert self.line is None
self._setup_indices_polygon()
self.type = 'polygon'
elif (self.line is not None):
assert self.indices is None
assert self.radius is None
assert self.center is None
assert self.polygon is None
self._setup_indices_line()
self.type = 'line'
elif (self.poly is not None):
# could be either a line or a polygon
# This is essentially for backwards compatibility
assert self.indices is None
assert self.radius is None
assert self.center is None
assert self.polygon is None
assert self.line is None
self.poly = num.asarray(self.poly)
if len(self.poly) > 2:
self.polygon = self.poly
self._setup_indices_polygon()
self.type = 'polygon'
else:
self.line = self.poly
self._setup_indices_line()
self.type = 'line'
else:
assert self.indices is None or self.indices is []
if self.indices is None:
self.full_indices = num.where(self.domain.tri_full_flag ==1)[0]
elif len(self.indices) == 0:
self.full_indices = []
else:
self.full_indices = num.array(self.indices)[self.domain.tri_full_flag[self.indices]==1]
def __repr__(self):
return "%s(%r)" % (self.__class__, self.__dict__)
def get_type(self):
return self.type
def plot_region(self, filename=None):
try:
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.tri as tri
except:
msg ="Couldn't import module from matplotlib, probably you need to update matplotlib"
raise msg
vertices = self.domain.get_vertex_coordinates()
full_mask = num.repeat(self.domain.tri_full_flag == 1, 3)
region_mask = num.zeros(len(self.domain),int).astype(bool)
region_mask[self.indices] = True
region_mask = num.repeat(region_mask,3)
#num.repeat(self.indices, 3)
#pprint(region_mask)
# Gather full and region nodes
fx = vertices[full_mask,0]
fy = vertices[full_mask,1]
gx = vertices[region_mask,0]
gy = vertices[region_mask,1]
# Plot mesh
n = len(fx) // 3
triang = num.array(list(range(0,3*n)))
triang.shape = (n, 3)
plt.triplot(fx, fy, triang, 'b-')
# Plot region
n = len(gx) // 3
if n > 0:
triang = num.array(list(range(0,3*n)))
triang.shape = (n, 3)
plt.triplot(gx, gy, triang, 'r-')
# Save triangulation to location pointed by filename
if filename is not None: plt.savefig(filename)
plt.show()
def _setup_indices_circle(self):
# Determine indices in circular region
N = self.domain.get_number_of_triangles()
points = self.domain.get_centroid_coordinates(absolute=True)
indices = []
c = self.center
r = self.radius
intersect = False
for k in range(N):
x, y = points[k,:] # Centroid
if ((x-c[0])**2+(y-c[1])**2) < r**2:
intersect = True
indices.append(k)
if len(indices) != 0:
self.indices = indices
else:
self.indices = num.asarray(indices)
if not self.domain.parallel:
msg = 'No centroids intersect circle center'+str(c)+' radius '+str(r)
if not intersect: raise Exception(msg)
def _setup_indices_polygon(self):
# Determine indices for polygonal region
points = self.domain.get_centroid_coordinates(absolute=True)
vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
indices = inside_polygon(points, self.polygon)
if self.expand_polygon :
n = len(self.polygon)
for j in range(n):
tris_0 = line_intersect(vertex_coordinates,
[self.polygon[j],self.polygon[(j+1)%n]])
indices = num.union1d(tris_0, indices)
if len(indices) == 0:
self.indices = indices
else:
self.indices = num.asarray(indices)
if not self.domain.parallel:
# only warn if not parallel as we should get lots of subdomains without indices
if len(indices) == 0:
msg = 'No centroids found for polygon %s '% str(self.polygon)
import warnings
warnings.warn(msg)
def _setup_indices_line(self):
# Determine indices for triangles intersecting a line region
vertex_coordinates = self.domain.get_vertex_coordinates(absolute=True)
indices = line_intersect(vertex_coordinates, self.line)
if len(indices) == 0:
self.indices = indices
else:
self.indices = num.asarray(indices)
if not self.domain.parallel:
msg = 'No centroids intersecting line %s '% str(self.line)
if len(indices) == 0: raise Exception(msg)
def get_indices(self, full_only=True):
if full_only:
return self.full_indices
else:
return self.indices
def set_verbose(self, verbose=True):
self.verbose = verbose
class Centroid_field(object):
def __init__(self, region, value, verbose=None):
self.region = region
self.value = value
self.domain = self.region.domain
def set_value(self,value):
"""Set value
Can change value while running
Can be a scalar, or a function of t or x,y or x,y,t or a quantity
"""
# Test if rate is a quantity
if isinstance(value, Quantity):
self.value_type = 'quantity'
else:
# Possible types are 'scalar', 't', 'x,y' and 'x,y,t'
from anuga.utilities.function_utils import determine_function_type
self.value_type = determine_function_type(value)
self.value = value
if self.value_type == 'scalar':
self.value_callable = False
self.value_spatial = False
elif self.value_type == 'quantity':
self.value_callable = False
self.valuespatial = False
elif self.value_type == 't':
self.value_callable = True
self.value_spatial = False
else:
self.value_callable = True
self.value_spatial = True