Context Navigation

source: inundation-numpy-branch/pyvolution/quantity.py @ 3381

Last change on this file since 3381 was 2608, checked in by ole, 19 years ago
Played with custom exceptions for ANUGA
File size: 43.1 KB

Rev	Line
[2526]	1	"""Class Quantity - Implements values at each triangular element
	2
	3	To create:
	4
	5	Quantity(domain, vertex_values)
	6
	7	domain: Associated domain structure. Required.
	8
	9	vertex_values: N x 3 array of values at each vertex for each element.
	10	Default None
	11
	12	If vertex_values are None Create array of zeros compatible with domain.
	13	Otherwise check that it is compatible with dimenions of domain.
	14	Otherwise raise an exception
	15	"""
	16
[2608]	17	from numpy import concatenate, zeros, Float, Int, array, reshape, ArrayType, allclose, take, sum, equal
[2526]	18
[2608]	19
[2526]	20	class Quantity:
	21
	22	def __init__(self, domain, vertex_values=None):
	23
	24	from mesh import Mesh
[2608]	25	#from Numeric import array, zeros, Float
[2526]	26
	27	msg = 'First argument in Quantity.__init__ '
	28	msg += 'must be of class Mesh (or a subclass thereof)'
	29	assert isinstance(domain, Mesh), msg
	30
	31	if vertex_values is None:
	32	N = domain.number_of_elements
	33	self.vertex_values = zeros((N, 3), Float)
	34	else:
	35	self.vertex_values = array(vertex_values).astype(Float)
	36
	37	N, V = self.vertex_values.shape
	38	assert V == 3,\
	39	'Three vertex values per element must be specified'
	40
	41
	42	msg = 'Number of vertex values (%d) must be consistent with'\
	43	%N
	44	msg += 'number of elements in specified domain (%d).'\
	45	%domain.number_of_elements
	46
	47	assert N == domain.number_of_elements, msg
	48
	49	self.domain = domain
	50
	51	#Allocate space for other quantities
	52	self.centroid_values = zeros(N, Float)
	53	self.edge_values = zeros((N, 3), Float)
	54
	55	#Intialise centroid and edge_values
	56	self.interpolate()
	57
	58
	59
	60	#Methods for operator overloading
	61	def __len__(self):
	62	return self.centroid_values.shape[0]
	63
	64
	65	def __neg__(self):
	66	"""Negate all values in this quantity giving meaning to the
	67	expression -Q where Q is an instance of class Quantity
	68	"""
	69
	70	Q = Quantity(self.domain)
	71	Q.set_values(-self.vertex_values)
	72	return Q
	73
	74
	75	def __add__(self, other):
	76	"""Add to self anything that could populate a quantity
	77
	78	E.g other can be a constant, an array, a function, another quantity
	79	(except for a filename or points, attributes (for now))
	80	- see set_values for details
	81	"""
	82
	83	Q = Quantity(self.domain)
	84	Q.set_values(other)
	85
	86	result = Quantity(self.domain)
	87	result.set_values(self.vertex_values + Q.vertex_values)
	88	return result
	89
	90	def __radd__(self, other):
	91	"""Handle cases like 7+Q, where Q is an instance of class Quantity
	92	"""
	93	return self + other
	94
	95
	96	def __sub__(self, other):
	97	return self + -other #Invoke __neg__
	98
	99	def __mul__(self, other):
	100	"""Multiply self with anything that could populate a quantity
	101
	102	E.g other can be a constant, an array, a function, another quantity
	103	(except for a filename or points, attributes (for now))
	104	- see set_values for details
	105
	106	Note that if two quantitites q1 and q2 are multiplied,
	107	vertex values are multiplied entry by entry
	108	while centroid and edge values are re-interpolated.
	109	Hence they won't be the product of centroid or edge values
	110	from q1 and q2.
	111	"""
	112
	113	Q = Quantity(self.domain)
	114	Q.set_values(other)
	115
	116	result = Quantity(self.domain)
	117	result.set_values(self.vertex_values * Q.vertex_values)
	118	return result
	119
	120	def __rmul__(self, other):
	121	"""Handle cases like 3*Q, where Q is an instance of class Quantity
	122	"""
	123	return self * other
	124
	125	def __pow__(self, other):
	126	"""Raise quantity to (numerical) power
	127
	128	As with __mul__ vertex values are processed entry by entry
	129	while centroid and edge values are re-interpolated.
	130
	131	Example using __pow__:
	132	Q = (Q12 + Q22)**0.5
	133
	134	"""
	135
	136	result = Quantity(self.domain)
	137	result.set_values(self.vertex_values**other)
	138	return result
	139
	140
	141
	142	def interpolate(self):
	143	"""Compute interpolated values at edges and centroid
	144	Pre-condition: vertex_values have been set
	145	"""
	146
	147	N = self.vertex_values.shape[0]
	148	for i in range(N):
	149	v0 = self.vertex_values[i, 0]
	150	v1 = self.vertex_values[i, 1]
	151	v2 = self.vertex_values[i, 2]
	152
	153	self.centroid_values[i] = (v0 + v1 + v2)/3
	154
	155	self.interpolate_from_vertices_to_edges()
	156
	157
	158	def interpolate_from_vertices_to_edges(self):
	159	#Call correct module function
	160	#(either from this module or C-extension)
	161	interpolate_from_vertices_to_edges(self)
	162
	163
	164
	165
	166	#New leaner interface to setting values
	167	def set_values(self,
	168	numeric = None, # List, numeric array or constant
	169	quantity = None, # Another quantity
	170	function = None, # Callable object: f(x,y)
	171	geospatial_data = None, #Arbitrary dataset
	172	points = None, values = None, data_georef = None, #Input for least squares (obsoleted by use of geo_spatial object)
	173	filename = None, attribute_name = None, #Input from file
	174	alpha = None,
	175	location = 'vertices',
	176	indices = None,
	177	verbose = False,
	178	use_cache = False):
	179
	180	"""Set values for quantity based on different sources.
	181
	182	numeric:
	183	Compatible list, Numeric array (see below) or constant.
	184	If callable it will treated as a function (see below)
	185	If instance of another Quantity it will be treated as such.
	186	If geo_spatial object it will be treated as such
	187
	188	quantity:
	189	Another quantity (compatible quantity, e.g. obtained as a
	190	linear combination of quantities)
	191
	192	function:
	193	Any callable object that takes two 1d arrays x and y
	194	each of length N and returns an array also of length N.
	195	The function will be evaluated at points determined by
	196	location and indices in the underlying mesh.
	197
	198	geospatial_data:
	199	Arbitrary geo spatial dataset in the form of the class
	200	Geospatial_data. Mesh points are populated using least squares
	201	fitting
	202
	203	points:
	204	Nx2 array of data points for use with least squares fit
	205	If points are present, an N array of attribute
	206	values corresponding to
	207	each data point must be present.
	208
	209	values:
	210	If points is specified, values is an array of length N containing
	211	attribute values for each point.
	212
	213	data_georef:
	214	If points is specified, geo_reference applies to each point.
	215
	216	filename:
	217	Name of a .pts file containing data points and attributes for
	218	use with least squares.
	219
	220	attribute_name:
	221	If specified, any array matching that name
	222	will be used. from file or geospatial_data.
	223	Otherwise a default will be used.
	224
	225	alpha:
	226	Smoothing parameter to be used with least squares fits.
	227	See module least_squares for further details about alpha.
	228	Alpha will only be used with points, values or filename.
	229	Otherwise it will be ignored.
	230
	231
	232	location: Where values are to be stored.
	233	Permissible options are: vertices, edges, centroids
	234	Default is 'vertices'
	235
	236	In case of location == 'centroids' the dimension values must
	237	be a list of a Numerical array of length N,
	238	N being the number of elements.
	239	Otherwise it must be of dimension Nx3
	240
	241
	242	The values will be stored in elements following their
	243	internal ordering.
	244
	245	If location is not 'unique vertices' Indices is the
	246	set of element ids that the operation applies to.
	247	If location is 'unique vertices' Indices is the set
	248	of vertex ids that the operation applies to.
	249
	250	If selected location is vertices, values for
	251	centroid and edges will be assigned interpolated
	252	values. In any other case, only values for the
	253	specified locations will be assigned and the others
	254	will be left undefined.
	255
	256	verbose: True means that output to stdout is generated
	257
	258	use_cache: True means that caching of intermediate results is
	259	attempted for least squares fit.
	260
	261
	262
	263
	264	Exactly one of the arguments
	265	numeric, quantity, function, points, filename
	266	must be present.
	267	"""
	268
	269	from geospatial_data.geospatial_data import Geospatial_data
	270	from types import FloatType, IntType, LongType, ListType, NoneType
	271
	272	#General input checks
	273	L = [numeric, quantity, function, geospatial_data, points, filename]
	274	msg = 'Exactly one of the arguments '+\
	275	'numeric, quantity, function, geospatial_data, points, or filename '+\
	276	'must be present.'
	277	assert L.count(None) == len(L)-1, msg
	278
	279
	280	if location not in ['vertices', 'centroids', 'edges',
	281	'unique vertices']:
	282	msg = 'Invalid location: %s' %location
	283	raise msg
	284
	285
	286	msg = 'Indices must be a list or None'
	287	assert type(indices) in [ListType, NoneType, ArrayType], msg
	288
	289
	290	if not(points is None and values is None and data_georef is None):
	291	from warnings import warn
	292
	293	msg = 'Using points, values or data_georef with set_quantity '
	294	msg += 'is obsolete. Please use a Geospatial_data object instead.'
	295	warn(msg, DeprecationWarning)
	296
	297
	298
	299	#Determine which 'set_values_from_...' to use
	300
	301	if numeric is not None:
	302	if type(numeric) in [FloatType, IntType, LongType]:
	303	self.set_values_from_constant(numeric,
	304	location, indices, verbose)
	305	elif type(numeric) in [ArrayType, ListType]:
	306	self.set_values_from_array(numeric,
	307	location, indices, verbose)
	308	elif callable(numeric):
	309	self.set_values_from_function(numeric,
	310	location, indices, verbose)
	311	elif isinstance(numeric, Quantity):
	312	self.set_values_from_quantity(numeric,
	313	location, indices, verbose)
	314	elif isinstance(numeric, Geospatial_data):
	315	self.set_values_from_geospatial_data(numeric,
	316	alpha,
	317	location, indices,
	318	verbose = verbose,
	319	use_cache = use_cache)
	320	else:
	321	msg = 'Illegal type for argument numeric: %s' %str(numeric)
	322	raise msg
	323
	324	elif quantity is not None:
	325	self.set_values_from_quantity(quantity,
	326	location, indices, verbose)
	327	elif function is not None:
	328	msg = 'Argument function must be callable'
	329	assert callable(function), msg
	330	self.set_values_from_function(function,
	331	location, indices, verbose)
	332	elif geospatial_data is not None:
	333	self.set_values_from_geospatial_data(geospatial_data,
	334	alpha,
	335	location, indices,
	336	verbose = verbose,
	337	use_cache = use_cache)
	338	elif points is not None:
	339	print 'The usage of points in set_values will be deprecated.' +\
	340	'Please use the geospatial_data object.'
	341
	342	msg = 'When points are specified, associated values must also be.'
	343	assert values is not None, msg
	344	self.set_values_from_points(points, values, alpha,
	345	location, indices,
	346	data_georef = data_georef,
	347	verbose = verbose,
	348	use_cache = use_cache)
	349	elif filename is not None:
	350	self.set_values_from_file(filename, attribute_name, alpha,
	351	location, indices,
	352	verbose = verbose,
	353	use_cache = use_cache)
	354	else:
	355	raise 'This can\'t happen :-)'
	356
	357
	358	#Update all locations in triangles
	359	if location == 'vertices' or location == 'unique vertices':
	360	#Intialise centroid and edge_values
	361	self.interpolate()
	362
	363	if location == 'centroids':
	364	#Extrapolate 1st order - to capture notion of area being specified
	365	self.extrapolate_first_order()
	366
	367
	368
	369	#Specific functions for setting values
	370	def set_values_from_constant(self, X,
	371	location, indices, verbose):
	372	"""Set quantity values from specified constant X
	373	"""
	374
	375
	376	if location == 'centroids':
	377	if (indices == None):
	378	self.centroid_values[:] = X
	379	else:
	380	#Brute force
	381	for i in indices:
	382	self.centroid_values[i,:] = X
	383
	384	elif location == 'edges':
	385	if (indices == None):
	386	self.edge_values[:] = X
	387	else:
	388	#Brute force
	389	for i in indices:
	390	self.edge_values[i,:] = X
	391
	392	elif location == 'unique vertices':
	393	if (indices == None):
	394	self.edge_values[:] = X
	395	else:
	396
	397	#Go through list of unique vertices
	398	for unique_vert_id in indices:
	399	triangles = self.domain.vertexlist[unique_vert_id]
	400
	401	#In case there are unused points
	402	if triangles is None: continue
	403
	404	#Go through all triangle, vertex pairs
	405	#and set corresponding vertex value
	406	for triangle_id, vertex_id in triangles:
	407	self.vertex_values[triangle_id, vertex_id] = X
	408
	409	#Intialise centroid and edge_values
	410	self.interpolate()
	411	else:
	412	if (indices == None):
	413	self.vertex_values[:] = X
	414	else:
	415	#Brute force
	416	for i_vertex in indices:
	417	self.vertex_values[i_vertex,:] = X
	418
	419
	420
	421
	422
	423
	424	def set_values_from_array(self, values,
	425	location, indices, verbose):
	426	"""Set values for quantity
	427
	428	values: Numeric array
	429	location: Where values are to be stored.
	430	Permissible options are: vertices, edges, centroid, unique vertices
	431	Default is 'vertices'
	432
	433	indices - if this action is carried out on a subset of
	434	elements or unique vertices
	435	The element/unique vertex indices are specified here.
	436
	437	In case of location == 'centroid' the dimension values must
	438	be a list of a Numerical array of length N, N being the number
	439	of elements.
	440
	441	Otherwise it must be of dimension Nx3
	442
	443	The values will be stored in elements following their
	444	internal ordering.
	445
	446	If selected location is vertices, values for centroid and edges
	447	will be assigned interpolated values.
	448	In any other case, only values for the specified locations
	449	will be assigned and the others will be left undefined.
	450	"""
	451
	452	values = array(values).astype(Float)
	453
	454	if indices is not None:
	455	indices = array(indices).astype(Int)
	456	msg = 'Number of values must match number of indices'
	457	assert values.shape[0] == indices.shape[0], msg
	458
	459	N = self.centroid_values.shape[0]
	460
	461	if location == 'centroids':
	462	assert len(values.shape) == 1, 'Values array must be 1d'
	463
	464	if indices is None:
	465	msg = 'Number of values must match number of elements'
	466	assert values.shape[0] == N, msg
	467
	468	self.centroid_values = values
	469	else:
	470	msg = 'Number of values must match number of indices'
	471	assert values.shape[0] == indices.shape[0], msg
	472
	473	#Brute force
	474	for i in range(len(indices)):
	475	self.centroid_values[indices[i]] = values[i]
	476
	477	elif location == 'edges':
	478	assert len(values.shape) == 2, 'Values array must be 2d'
	479
	480	msg = 'Number of values must match number of elements'
	481	assert values.shape[0] == N, msg
	482
	483	msg = 'Array must be N x 3'
	484	assert values.shape[1] == 3, msg
	485
	486	self.edge_values = values
	487
	488	elif location == 'unique vertices':
	489	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
	490	'Values array must be 1d'
	491
	492	self.set_vertex_values(values.flat, indices = indices)
	493	else:
	494	if len(values.shape) == 1:
	495	self.set_vertex_values(values, indices = indices)
	496	#if indices == None:
	497	#Values are being specified once for each unique vertex
	498	# msg = 'Number of values must match number of vertices'
	499	# assert values.shape[0] == self.domain.coordinates.shape[0], msg
	500	# self.set_vertex_values(values)
	501	#else:
	502	# for element_index, value in map(None, indices, values):
	503	# self.vertex_values[element_index, :] = value
	504
	505	elif len(values.shape) == 2:
	506	#Vertex values are given as a triplet for each triangle
	507
	508	msg = 'Array must be N x 3'
	509	assert values.shape[1] == 3, msg
	510
	511	if indices == None:
	512	self.vertex_values = values
	513	else:
	514	for element_index, value in map(None, indices, values):
	515	self.vertex_values[element_index] = value
	516	else:
	517	msg = 'Values array must be 1d or 2d'
	518	raise msg
	519
	520	def set_values_from_quantity(self, q,
	521	location, indices, verbose):
	522	"""Set quantity values from specified quantity instance q
	523
	524	Location is ignored
	525	"""
	526
	527
	528	A = q.vertex_values
	529
	530	msg = 'Quantities are defined on different meshes. '+\
	531	'This might be a case for implementing interpolation '+\
	532	'between different meshes.'
	533	assert allclose(A.shape, self.vertex_values.shape), msg
	534
	535	self.set_values(A, location='vertices',
	536	indices=indices,
	537	verbose=verbose)
	538
	539
	540	def set_values_from_function(self, f,
	541	location, indices, verbose):
	542	"""Set values for quantity using specified function
	543
	544	f: x, y -> z Function where x, y and z are arrays
	545	location: Where values are to be stored.
	546	Permissible options are: vertices, centroid, edges,
	547	unique vertices
	548	Default is "vertices"
	549	"""
	550
	551	#FIXME: Should check that function returns something sensible and
	552	#raise a meaningfull exception if it returns None for example
	553
	554	#FIXME: Should supply absolute coordinates
	555
	556	if (indices is None):
	557	indices = range(len(self))
	558	is_subset = False
	559	else:
	560	is_subset = True
	561
	562	if location == 'centroids':
	563	P = take(self.domain.centroid_coordinates, indices)
	564	if is_subset:
	565	self.set_values(f(P[:,0], P[:,1]),
	566	location = location,
	567	indices = indices)
	568	else:
	569	self.set_values(f(P[:,0], P[:,1]), location = location)
	570	elif location == 'vertices':
	571	P = self.domain.vertex_coordinates
	572	if is_subset:
	573	#Brute force
	574	for e in indices:
	575	for i in range(3):
	576	self.vertex_values[e,i] = f(P[e,2i], P[e,2i+1])
	577	else:
	578	for i in range(3):
	579	self.vertex_values[:,i] = f(P[:,2i], P[:,2i+1])
	580	else:
	581	raise 'Not implemented: %s' %location
	582
	583
	584
	585	def set_values_from_geospatial_data(self, geospatial_data, alpha,
	586	location, indices,
	587	verbose = False,
	588	use_cache = False):
	589
	590	#FIXME: Use this function for the time being. Later move code in here
	591
	592	points = geospatial_data.get_data_points(absolute = False)
	593	values = geospatial_data.get_attributes()
	594	data_georef = geospatial_data.get_geo_reference()
	595
	596
	597
	598	self.set_values_from_points(points, values, alpha,
	599	location, indices,
	600	data_georef = data_georef,
	601	verbose = verbose,
	602	use_cache = use_cache)
	603
	604
	605
	606	def set_values_from_points(self, points, values, alpha,
	607	location, indices,
	608	data_georef = None,
	609	verbose = False,
	610	use_cache = False):
	611	"""Set quantity values from arbitray data points using least squares
	612	"""
	613
	614
	615	from utilities.numerical_tools import ensure_numeric
	616	from least_squares import fit_to_mesh
	617	from coordinate_transforms.geo_reference import Geo_reference
	618
	619
	620	points = ensure_numeric(points, Float)
	621	values = ensure_numeric(values, Float)
	622
	623	if location != 'vertices':
	624	msg = 'set_values_from_points is only defined for '+\
	625	'location=\'vertices\''
	626	raise msg
	627
	628	coordinates = self.domain.coordinates
	629	triangles = self.domain.triangles
	630
	631
	632	#Take care of georeferencing
	633	if data_georef is None:
	634	data_georef = Geo_reference()
	635
	636
	637	mesh_georef = self.domain.geo_reference
	638
	639	#print mesh_georef
	640	#print data_georef
	641	#print points
	642
	643
	644	#Call least squares method
	645	args = (coordinates, triangles, points, values)
	646	kwargs = {'data_origin': data_georef.get_origin(),
	647	'mesh_origin': mesh_georef.get_origin(),
	648	'alpha': alpha,
	649	'precrop': True,
	650	'verbose': verbose}
	651
	652	#print kwargs
	653
	654	if use_cache is True:
	655	try:
	656	from caching import cache
	657	except:
	658	msg = 'Caching was requested, but caching module'+\
	659	'could not be imported'
	660	raise msg
	661
	662	vertex_attributes = cache(fit_to_mesh,
	663	args, kwargs,
	664	verbose = verbose)
	665	else:
	666
	667	vertex_attributes = apply(fit_to_mesh,
	668	args, kwargs)
	669
	670	#Call underlying method using array values
	671	self.set_values_from_array(vertex_attributes,
	672	location, indices, verbose)
	673
	674
	675
	676	def set_values_from_file(self, filename, attribute_name, alpha,
	677	location, indices,
	678	verbose = False,
	679	use_cache = False):
	680	"""Set quantity based on arbitrary points in .pts file
	681	using least_squares attribute_name selects name of attribute
	682	present in file.
	683	If not specified try to use whatever is available in file.
	684	"""
	685
	686	from load_mesh.loadASCII import import_points_file
	687	from geospatial_data.geospatial_data import points_dictionary2geospatial_data
	688
	689	from types import StringType
	690	msg = 'Filename must be a text string'
	691	assert type(filename) == StringType, msg
	692
	693
	694	#Read from (NetCDF) file
	695	#FIXME (Ole): This function should really return a Geospatial_data object.
	696	points_dict = import_points_file(filename)
	697	points = points_dict['pointlist']
	698	attributes = points_dict['attributelist']
	699
	700	if attribute_name is None:
	701	names = attributes.keys()
	702	attribute_name = names[0]
	703
	704	msg = 'Attribute_name must be a text string'
	705	assert type(attribute_name) == StringType, msg
	706
	707
	708	if verbose:
	709	print 'Using attribute %s from file %s' %(attribute_name, filename)
	710	print 'Available attributes: %s' %(names)
	711
	712	#try:
	713	# z = attributes[attribute_name]
	714	#except:
	715	# msg = 'Could not extract attribute %s from file %s'\
	716	# %(attribute_name, filename)
	717	# raise msg
	718
	719
	720	#Take care of georeferencing
	721	if points_dict.has_key('geo_reference') and \
	722	points_dict['geo_reference'] is not None:
	723	data_georef = points_dict['geo_reference']
	724	else:
	725	data_georef = None
	726
	727
	728
	729	#Call underlying method for geospatial data
	730	geospatial_data = points_dictionary2geospatial_data(points_dict)
	731	geospatial_data.set_default_attribute_name(attribute_name)
	732
	733	self.set_values_from_geospatial_data(geospatial_data,
	734	alpha,
	735	location, indices,
	736	verbose = verbose,
	737	use_cache = use_cache)
	738
	739	#Call underlying method for points
	740	#self.set_values_from_points(points, z, alpha,
	741	# location, indices,
	742	# data_georef = data_georef,
	743	# verbose = verbose,
	744	# use_cache = use_cache)
	745
	746
	747
	748	def get_values(self, location='vertices', indices = None):
	749	"""get values for quantity
	750
	751	return X, Compatible list, Numeric array (see below)
	752	location: Where values are to be stored.
	753	Permissible options are: vertices, edges, centroid
	754	and unique vertices. Default is 'vertices'
	755
	756	In case of location == 'centroids' the dimension values must
	757	be a list of a Numerical array of length N, N being the number
	758	of elements. Otherwise it must be of dimension Nx3
	759
	760	The returned values with be a list the length of indices
	761	(N if indices = None). Each value will be a list of the three
	762	vertex values for this quantity.
	763
	764	Indices is the set of element ids that the operation applies to.
	765
	766	The values will be stored in elements following their
	767	internal ordering.
	768
	769	"""
	770
	771	if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
	772	msg = 'Invalid location: %s' %location
	773	raise msg
	774
[2608]	775	import types
[2526]	776	assert type(indices) in [types.ListType, types.NoneType,
[2608]	777	ArrayType],\
[2526]	778	'Indices must be a list or None'
	779
	780	if location == 'centroids':
	781	if (indices == None):
	782	indices = range(len(self))
	783	return take(self.centroid_values,indices)
	784	elif location == 'edges':
	785	if (indices == None):
	786	indices = range(len(self))
	787	return take(self.edge_values,indices)
	788	elif location == 'unique vertices':
	789	if (indices == None):
	790	indices=range(self.domain.coordinates.shape[0])
	791	vert_values = []
	792	#Go through list of unique vertices
	793	for unique_vert_id in indices:
	794	triangles = self.domain.vertexlist[unique_vert_id]
	795
	796	#In case there are unused points
	797	if triangles is None:
	798	msg = 'Unique vertex not associated with triangles'
	799	raise msg
	800
	801	# Go through all triangle, vertex pairs
	802	# Average the values
	803	sum = 0
	804	for triangle_id, vertex_id in triangles:
	805	sum += self.vertex_values[triangle_id, vertex_id]
	806	vert_values.append(sum/len(triangles))
[2608]	807	return array(vert_values)
[2526]	808	else:
	809	if (indices == None):
	810	indices = range(len(self))
	811	return take(self.vertex_values,indices)
	812
	813
	814
	815	def set_vertex_values(self, A, indices = None):
	816	"""Set vertex values for all unique vertices based on input array A
	817	which has one entry per unique vertex, i.e.
	818	one value for each row in array self.domain.coordinates or
	819	one value for each row in vertexlist.
	820
	821	indices is the list of vertex_id's that will be set.
	822
	823	This function is used by set_values_from_array
	824	"""
	825
	826	#Assert that A can be converted to a Numeric array of appropriate dim
	827	A = array(A, Float)
	828
	829	#print 'SHAPE A', A.shape
	830	assert len(A.shape) == 1
	831
	832	if indices == None:
	833	assert A.shape[0] == self.domain.coordinates.shape[0]
	834	vertex_list = range(A.shape[0])
	835	else:
	836	assert A.shape[0] == len(indices)
	837	vertex_list = indices
	838
	839	#Go through list of unique vertices
	840	for i_index, unique_vert_id in enumerate(vertex_list):
	841	triangles = self.domain.vertexlist[unique_vert_id]
	842
	843	if triangles is None: continue #In case there are unused points
	844
	845	#Go through all triangle, vertex pairs
	846	#touching vertex unique_vert_id and set corresponding vertex value
	847	for triangle_id, vertex_id in triangles:
	848	self.vertex_values[triangle_id, vertex_id] = A[i_index]
	849
	850	#Intialise centroid and edge_values
	851	self.interpolate()
	852
	853
	854	def smooth_vertex_values(self, value_array='field_values',
	855	precision = None):
	856	""" Smooths field_values or conserved_quantities data.
	857	TODO: be able to smooth individual fields
	858	NOTE: This function does not have a test.
	859	FIXME: NOT DONE - do we need it?
	860	FIXME: this function isn't called by anything.
	861	Maybe it should be removed..-DSG
	862	"""
	863
	864
	865	A,V = self.get_vertex_values(xy=False,
	866	value_array=value_array,
	867	smooth = True,
	868	precision = precision)
	869
	870	#Set some field values
	871	for volume in self:
	872	for i,v in enumerate(volume.vertices):
	873	if value_array == 'field_values':
	874	volume.set_field_values('vertex', i, A[v,:])
	875	elif value_array == 'conserved_quantities':
	876	volume.set_conserved_quantities('vertex', i, A[v,:])
	877
	878	if value_array == 'field_values':
	879	self.precompute()
	880	elif value_array == 'conserved_quantities':
	881	Volume.interpolate_conserved_quantities()
	882
	883
	884	#Method for outputting model results
	885	#FIXME: Split up into geometric and numeric stuff.
	886	#FIXME: Geometric (X,Y,V) should live in mesh.py
	887	#FIXME: STill remember to move XY to mesh
	888	def get_vertex_values(self,
	889	xy=True,
	890	smooth = None,
	891	precision = None,
	892	reduction = None):
	893	"""Return vertex values like an OBJ format
	894
	895	The vertex values are returned as one sequence in the 1D float array A.
	896	If requested the coordinates will be returned in 1D arrays X and Y.
	897
	898	The connectivity is represented as an integer array, V, of dimension
	899	M x 3, where M is the number of volumes. Each row has three indices
	900	into the X, Y, A arrays defining the triangle.
	901
	902	if smooth is True, vertex values corresponding to one common
	903	coordinate set will be smoothed according to the given
	904	reduction operator. In this case vertex coordinates will be
	905	de-duplicated.
	906
	907	If no smoothings is required, vertex coordinates and values will
	908	be aggregated as a concatenation of values at
	909	vertices 0, vertices 1 and vertices 2
	910
	911
	912	Calling convention
	913	if xy is True:
	914	X,Y,A,V = get_vertex_values
	915	else:
	916	A,V = get_vertex_values
	917
	918	"""
	919
	920	if smooth is None:
	921	smooth = self.domain.smooth
	922
	923	if precision is None:
	924	precision = Float
	925
	926	if reduction is None:
	927	reduction = self.domain.reduction
	928
	929	#Create connectivity
	930
	931	if smooth == True:
	932
	933	V = self.domain.get_vertices()
	934	N = len(self.domain.vertexlist)
	935	A = zeros(N, precision)
	936
	937	#Smoothing loop
	938	for k in range(N):
	939	L = self.domain.vertexlist[k]
	940
	941	#Go through all triangle, vertex pairs
	942	#contributing to vertex k and register vertex value
	943
	944	if L is None: continue #In case there are unused points
	945
	946	contributions = []
	947	for volume_id, vertex_id in L:
	948	v = self.vertex_values[volume_id, vertex_id]
	949	contributions.append(v)
	950
	951	A[k] = reduction(contributions)
	952
	953
	954	if xy is True:
	955	X = self.domain.coordinates[:,0].astype(precision)
	956	Y = self.domain.coordinates[:,1].astype(precision)
	957
	958	return X, Y, A, V
	959	else:
	960	return A, V
	961	else:
	962	#Don't smooth
	963	#obj machinery moved to general_mesh
	964
	965	# Create a V like [[0 1 2], [3 4 5]....[3m-2 3m-1 3*m]]
	966	# These vert_id's will relate to the verts created below
	967	#m = len(self.domain) #Number of volumes
	968	#M = 3*m #Total number of unique vertices
	969	#V = reshape(array(range(M)).astype(Int), (m,3))
	970
	971	V = self.domain.get_triangles(obj=True)
	972	#FIXME use get_vertices, when ready
	973
	974	A = self.vertex_values.flat
	975
	976	#Do vertex coordinates
	977	if xy is True:
	978	C = self.domain.get_vertex_coordinates()
	979
	980	X = C[:,0:6:2].copy()
	981	Y = C[:,1:6:2].copy()
	982
	983	return X.flat, Y.flat, A, V
	984	else:
	985	return A, V
	986
	987
	988	def extrapolate_first_order(self):
	989	"""Extrapolate conserved quantities from centroid to
	990	vertices for each volume using
	991	first order scheme.
	992	"""
	993
	994	qc = self.centroid_values
	995	qv = self.vertex_values
	996
	997	for i in range(3):
	998	qv[:,i] = qc
	999
	1000
	1001	def get_integral(self):
	1002	"""Compute the integral of quantity across entire domain
	1003	"""
	1004	integral = 0
	1005	for k in range(self.domain.number_of_elements):
	1006	area = self.domain.areas[k]
	1007	qc = self.centroid_values[k]
	1008	integral += qc*area
	1009
	1010	return integral
	1011
	1012
	1013
	1014
	1015	class Conserved_quantity(Quantity):
	1016	"""Class conserved quantity adds to Quantity:
	1017
	1018	boundary values, storage and method for updating, and
	1019	methods for (second order) extrapolation from centroid to vertices inluding
	1020	gradients and limiters
	1021	"""
	1022
	1023	def __init__(self, domain, vertex_values=None):
	1024	Quantity.__init__(self, domain, vertex_values)
	1025
	1026
	1027	#Allocate space for boundary values
	1028	L = len(domain.boundary)
	1029	self.boundary_values = zeros(L, Float)
	1030
	1031	#Allocate space for updates of conserved quantities by
	1032	#flux calculations and forcing functions
	1033
	1034	N = domain.number_of_elements
	1035	self.explicit_update = zeros(N, Float )
	1036	self.semi_implicit_update = zeros(N, Float )
	1037
	1038
	1039	def update(self, timestep):
	1040	#Call correct module function
	1041	#(either from this module or C-extension)
	1042	return update(self, timestep)
	1043
	1044
	1045	def compute_gradients(self):
	1046	#Call correct module function
	1047	#(either from this module or C-extension)
	1048	return compute_gradients(self)
	1049
	1050
	1051	def limit(self):
	1052	#Call correct module function
	1053	#(either from this module or C-extension)
	1054	limit(self)
	1055
	1056
	1057	def extrapolate_second_order(self):
	1058	#Call correct module function
	1059	#(either from this module or C-extension)
	1060	extrapolate_second_order(self)
	1061
	1062
	1063	def update(quantity, timestep):
	1064	"""Update centroid values based on values stored in
	1065	explicit_update and semi_implicit_update as well as given timestep
	1066
	1067	Function implementing forcing terms must take on argument
	1068	which is the domain and they must update either explicit
	1069	or implicit updates, e,g,:
	1070
	1071	def gravity(domain):
	1072	....
	1073	domain.quantities['xmomentum'].explicit_update = ...
	1074	domain.quantities['ymomentum'].explicit_update = ...
	1075
	1076
	1077
	1078	Explicit terms must have the form
	1079
	1080	G(q, t)
	1081
	1082	and explicit scheme is
	1083
	1084	q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
	1085
	1086
	1087	Semi implicit forcing terms are assumed to have the form
	1088
	1089	G(q, t) = H(q, t) q
	1090
	1091	and the semi implicit scheme will then be
	1092
	1093	q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
	1094
	1095
	1096	"""
	1097
[2608]	1098	#from Numeric import sum, equal, ones, Float
[2526]	1099
	1100	N = quantity.centroid_values.shape[0]
	1101
	1102
	1103	#Divide H by conserved quantity to obtain G (see docstring above)
	1104
	1105
	1106	for k in range(N):
	1107	x = quantity.centroid_values[k]
	1108	if x == 0.0:
	1109	#FIXME: Is this right
	1110	quantity.semi_implicit_update[k] = 0.0
	1111	else:
	1112	quantity.semi_implicit_update[k] /= x
	1113
	1114	#Explicit updates
	1115	quantity.centroid_values += timestep*quantity.explicit_update
	1116
	1117	#Semi implicit updates
	1118	denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
	1119
	1120	if sum(equal(denominator, 0.0)) > 0.0:
	1121	msg = 'Zero division in semi implicit update. Call Stephen :-)'
	1122	raise msg
	1123	else:
	1124	#Update conserved_quantities from semi implicit updates
	1125	quantity.centroid_values /= denominator
	1126
	1127
	1128	def interpolate_from_vertices_to_edges(quantity):
	1129	"""Compute edge values from vertex values using linear interpolation
	1130	"""
	1131
	1132	for k in range(quantity.vertex_values.shape[0]):
	1133	q0 = quantity.vertex_values[k, 0]
	1134	q1 = quantity.vertex_values[k, 1]
	1135	q2 = quantity.vertex_values[k, 2]
	1136
	1137	quantity.edge_values[k, 0] = 0.5*(q1+q2)
	1138	quantity.edge_values[k, 1] = 0.5*(q0+q2)
	1139	quantity.edge_values[k, 2] = 0.5*(q0+q1)
	1140
	1141
	1142
	1143	def extrapolate_second_order(quantity):
	1144	"""Extrapolate conserved quantities from centroid to
	1145	vertices for each volume using
	1146	second order scheme.
	1147	"""
	1148
	1149	a, b = quantity.compute_gradients()
	1150
	1151	X = quantity.domain.get_vertex_coordinates()
	1152	qc = quantity.centroid_values
	1153	qv = quantity.vertex_values
	1154
	1155	#Check each triangle
	1156	for k in range(quantity.domain.number_of_elements):
	1157	#Centroid coordinates
	1158	x, y = quantity.domain.centroid_coordinates[k]
	1159
	1160	#vertex coordinates
	1161	x0, y0, x1, y1, x2, y2 = X[k,:]
	1162
	1163	#Extrapolate
	1164	qv[k,0] = qc[k] + a[k](x0-x) + b[k](y0-y)
	1165	qv[k,1] = qc[k] + a[k](x1-x) + b[k](y1-y)
	1166	qv[k,2] = qc[k] + a[k](x2-x) + b[k](y2-y)
	1167
	1168
	1169	def compute_gradients(quantity):
	1170	"""Compute gradients of triangle surfaces defined by centroids of
	1171	neighbouring volumes.
	1172	If one edge is on the boundary, use own centroid as neighbour centroid.
	1173	If two or more are on the boundary, fall back to first order scheme.
	1174	"""
	1175
	1176	from utilitites.numerical_tools import gradient
	1177
	1178	centroid_coordinates = quantity.domain.centroid_coordinates
	1179	surrogate_neighbours = quantity.domain.surrogate_neighbours
	1180	centroid_values = quantity.centroid_values
	1181	number_of_boundaries = quantity.domain.number_of_boundaries
	1182
	1183	N = centroid_values.shape[0]
	1184
	1185	a = zeros(N, Float)
	1186	b = zeros(N, Float)
	1187
	1188	for k in range(N):
	1189	if number_of_boundaries[k] < 2:
	1190	#Two or three true neighbours
	1191
	1192	#Get indices of neighbours (or self when used as surrogate)
	1193	k0, k1, k2 = surrogate_neighbours[k,:]
	1194
	1195	#Get data
	1196	q0 = centroid_values[k0]
	1197	q1 = centroid_values[k1]
	1198	q2 = centroid_values[k2]
	1199
	1200	x0, y0 = centroid_coordinates[k0] #V0 centroid
	1201	x1, y1 = centroid_coordinates[k1] #V1 centroid
	1202	x2, y2 = centroid_coordinates[k2] #V2 centroid
	1203
	1204	#Gradient
	1205	a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
	1206
	1207	elif number_of_boundaries[k] == 2:
	1208	#One true neighbour
	1209
	1210	#Get index of the one neighbour
	1211	for k0 in surrogate_neighbours[k,:]:
	1212	if k0 != k: break
	1213	assert k0 != k
	1214
	1215	k1 = k #self
	1216
	1217	#Get data
	1218	q0 = centroid_values[k0]
	1219	q1 = centroid_values[k1]
	1220
	1221	x0, y0 = centroid_coordinates[k0] #V0 centroid
	1222	x1, y1 = centroid_coordinates[k1] #V1 centroid
	1223
	1224	#Gradient
	1225	a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
	1226	else:
	1227	#No true neighbours -
	1228	#Fall back to first order scheme
	1229	pass
	1230
	1231
	1232	return a, b
	1233
	1234
	1235
	1236	def limit(quantity):
	1237	"""Limit slopes for each volume to eliminate artificial variance
	1238	introduced by e.g. second order extrapolator
	1239
	1240	This is an unsophisticated limiter as it does not take into
	1241	account dependencies among quantities.
	1242
	1243	precondition:
	1244	vertex values are estimated from gradient
	1245	postcondition:
	1246	vertex values are updated
	1247	"""
	1248
	1249	N = quantity.domain.number_of_elements
	1250
	1251	beta_w = quantity.domain.beta_w
	1252
	1253	qc = quantity.centroid_values
	1254	qv = quantity.vertex_values
	1255
	1256	#Find min and max of this and neighbour's centroid values
	1257	qmax = zeros(qc.shape, Float)
	1258	qmin = zeros(qc.shape, Float)
	1259
	1260	for k in range(N):
	1261	qmax[k] = qmin[k] = qc[k]
	1262	for i in range(3):
	1263	n = quantity.domain.neighbours[k,i]
	1264	if n >= 0:
	1265	qn = qc[n] #Neighbour's centroid value
	1266
	1267	qmin[k] = min(qmin[k], qn)
	1268	qmax[k] = max(qmax[k], qn)
	1269
	1270
	1271	#Diffences between centroids and maxima/minima
	1272	dqmax = qmax - qc
	1273	dqmin = qmin - qc
	1274
	1275	#Deltas between vertex and centroid values
	1276	dq = zeros(qv.shape, Float)
	1277	for i in range(3):
	1278	dq[:,i] = qv[:,i] - qc
	1279
	1280	#Phi limiter
	1281	for k in range(N):
	1282
	1283	#Find the gradient limiter (phi) across vertices
	1284	phi = 1.0
	1285	for i in range(3):
	1286	r = 1.0
	1287	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
	1288	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
	1289
	1290	phi = min( min(r*beta_w, 1), phi )
	1291
	1292	#Then update using phi limiter
	1293	for i in range(3):
	1294	qv[k,i] = qc[k] + phi*dq[k,i]
	1295
	1296
	1297
	1298	from utilities import compile
	1299	if compile.can_use_C_extension('quantity_ext.c'):
	1300	#Replace python version with c implementations
	1301
	1302	from quantity_ext import limit, compute_gradients,\
	1303	extrapolate_second_order, interpolate_from_vertices_to_edges, update

Note: See TracBrowser for help on using the repository browser.

Download in other formats: