Context Navigation

source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py @ 4768

Last change on this file since 4768 was 4768, checked in by ole, 17 years ago
Backlog of comments and errormessages
File size: 52.4 KB

Line
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
17	from Numeric import array, zeros, Float, less, concatenate, NewAxis,\
18	argmax, argmin, allclose, take, reshape
19
20	from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
21	from anuga.utilities.polygon import inside_polygon
22
23	from anuga.geospatial_data.geospatial_data import Geospatial_data
24	from anuga.fit_interpolate.fit import fit_to_mesh
25	from anuga.config import points_file_block_line_size as default_block_line_size
26
27	class Quantity:
28
29	def __init__(self, domain, vertex_values=None):
30
31	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
32
33	msg = 'First argument in Quantity.__init__ '
34	msg += 'must be of class Mesh (or a subclass thereof)'
35	assert isinstance(domain, Mesh), msg
36
37	if vertex_values is None:
38	N = len(domain) # number_of_elements
39	self.vertex_values = zeros((N, 3), Float)
40	else:
41	self.vertex_values = array(vertex_values).astype(Float)
42
43	N, V = self.vertex_values.shape
44	assert V == 3,\
45	'Three vertex values per element must be specified'
46
47
48	msg = 'Number of vertex values (%d) must be consistent with'\
49	%N
50	msg += 'number of elements in specified domain (%d).'\
51	%len(domain)
52
53	assert N == len(domain), msg
54
55	self.domain = domain
56
57	#Allocate space for other quantities
58	self.centroid_values = zeros(N, Float)
59	self.edge_values = zeros((N, 3), Float)
60
61	#Intialise centroid and edge_values
62	self.interpolate()
63
64
65
66	#Methods for operator overloading
67	def __len__(self):
68	return self.centroid_values.shape[0]
69
70
71	def __neg__(self):
72	"""Negate all values in this quantity giving meaning to the
73	expression -Q where Q is an instance of class Quantity
74	"""
75
76	Q = Quantity(self.domain)
77	Q.set_values(-self.vertex_values)
78	return Q
79
80
81	def __add__(self, other):
82	"""Add to self anything that could populate a quantity
83
84	E.g other can be a constant, an array, a function, another quantity
85	(except for a filename or points, attributes (for now))
86	- see set_values for details
87	"""
88
89	Q = Quantity(self.domain)
90	Q.set_values(other)
91
92	result = Quantity(self.domain)
93	result.set_values(self.vertex_values + Q.vertex_values)
94	return result
95
96	def __radd__(self, other):
97	"""Handle cases like 7+Q, where Q is an instance of class Quantity
98	"""
99	return self + other
100
101
102	def __sub__(self, other):
103	return self + -other #Invoke __neg__
104
105	def __mul__(self, other):
106	"""Multiply self with anything that could populate a quantity
107
108	E.g other can be a constant, an array, a function, another quantity
109	(except for a filename or points, attributes (for now))
110	- see set_values for details
111
112	Note that if two quantitites q1 and q2 are multiplied,
113	vertex values are multiplied entry by entry
114	while centroid and edge values are re-interpolated.
115	Hence they won't be the product of centroid or edge values
116	from q1 and q2.
117	"""
118
119	Q = Quantity(self.domain)
120	Q.set_values(other)
121
122	result = Quantity(self.domain)
123	result.set_values(self.vertex_values * Q.vertex_values)
124	return result
125
126	def __rmul__(self, other):
127	"""Handle cases like 3*Q, where Q is an instance of class Quantity
128	"""
129	return self * other
130
131	def __pow__(self, other):
132	"""Raise quantity to (numerical) power
133
134	As with __mul__ vertex values are processed entry by entry
135	while centroid and edge values are re-interpolated.
136
137	Example using __pow__:
138	Q = (Q12 + Q22)**0.5
139
140	"""
141
142	result = Quantity(self.domain)
143	result.set_values(self.vertex_values**other)
144	return result
145
146
147
148	def interpolate(self):
149	"""Compute interpolated values at edges and centroid
150	Pre-condition: vertex_values have been set
151	"""
152	# FIXME (Ole): This function is not called
153	# in real model runs. Consider removing it.
154
155	N = self.vertex_values.shape[0]
156	for i in range(N):
157	v0 = self.vertex_values[i, 0]
158	v1 = self.vertex_values[i, 1]
159	v2 = self.vertex_values[i, 2]
160
161	self.centroid_values[i] = (v0 + v1 + v2)/3
162
163	self.interpolate_from_vertices_to_edges()
164
165
166	def interpolate_from_vertices_to_edges(self):
167	#Call correct module function
168	#(either from this module or C-extension)
169	interpolate_from_vertices_to_edges(self)
170
171
172
173
174	# New leaner interface to setting values
175	def set_values(self,
176	numeric = None, # List, numeric array or constant
177	quantity = None, # Another quantity
178	function = None, # Callable object: f(x,y)
179	geospatial_data = None, # Arbitrary dataset
180	points = None, values = None, data_georef = None, # Obsoleted by use of geo_spatial object
181	filename = None, attribute_name = None, #Input from file
182	alpha = None,
183	location = 'vertices',
184	polygon = None,
185	indices = None,
186	verbose = False,
187	use_cache = False):
188
189	"""Set values for quantity based on different sources.
190
191	numeric:
192	Compatible list, Numeric array (see below) or constant.
193	If callable it will treated as a function (see below)
194	If instance of another Quantity it will be treated as such.
195	If geo_spatial object it will be treated as such
196
197	quantity:
198	Another quantity (compatible quantity, e.g. obtained as a
199	linear combination of quantities)
200
201	function:
202	Any callable object that takes two 1d arrays x and y
203	each of length N and returns an array also of length N.
204	The function will be evaluated at points determined by
205	location and indices in the underlying mesh.
206
207	geospatial_data:
208	Arbitrary geo spatial dataset in the form of the class
209	Geospatial_data. Mesh points are populated using
210	fit_interpolate.fit fitting
211
212	points:
213	Nx2 array of data points for use with fit_interpolate.fit
214	If points are present, an N array of attribute
215	values corresponding to
216	each data point must be present.
217	(Obsoleted by geospatial_data)
218
219	values:
220	If points is specified, values is an array of length N containing
221	attribute values for each point.
222	(Obsoleted by geospatial_data)
223
224	data_georef:
225	If points is specified, geo_reference applies to each point.
226	(Obsoleted by geospatial_data)
227
228	filename:
229	Name of a points file containing data points and attributes for
230	use with fit_interpolate.fit.
231
232	attribute_name:
233	If specified, any array matching that name
234	will be used. from file or geospatial_data.
235	Otherwise a default will be used.
236
237	alpha:
238	Smoothing parameter to be used with fit_interpolate.fit.
239	See module fit_interpolate.fit for further details about alpha.
240	Alpha will only be used with points, values or filename.
241	Otherwise it will be ignored.
242
243
244	location: Where values are to be stored.
245	Permissible options are: vertices, edges, centroids
246	Default is 'vertices'
247
248	In case of location == 'centroids' the dimension values must
249	be a list of a Numerical array of length N,
250	N being the number of elements.
251	Otherwise it must be of dimension Nx3
252
253
254	The values will be stored in elements following their
255	internal ordering.
256
257	If location is not 'unique vertices' Indices is the
258	set of element ids that the operation applies to.
259	If location is 'unique vertices' Indices is the set
260	of vertex ids that the operation applies to.
261
262	If selected location is vertices, values for
263	centroid and edges will be assigned interpolated
264	values. In any other case, only values for the
265	specified locations will be assigned and the others
266	will be left undefined.
267
268
269	polygon: Restrict update of quantity to locations that fall
270	inside polygon. Polygon works by selecting indices
271	and calling set_values recursively.
272	Polygon mode has only been implemented for
273	constant values so far.
274
275	indices: Restrict update of quantity to locations that are
276	identified by indices (e.g. node ids if location
277	is 'vertices')
278
279	verbose: True means that output to stdout is generated
280
281	use_cache: True means that caching of intermediate results is
282	attempted for fit_interpolate.fit.
283
284
285
286
287	Exactly one of the arguments
288	numeric, quantity, function, points, filename
289	must be present.
290	"""
291
292	from anuga.geospatial_data.geospatial_data import Geospatial_data
293	from types import FloatType, IntType, LongType, ListType, NoneType
294	from Numeric import ArrayType
295
296
297	# Treat special case: Polygon situation
298	# Location will be ignored and set to 'centroids'
299	# FIXME (Ole): This needs to be generalised and
300	# perhaps the notion of location and indices simplified
301
302	if polygon is not None:
303	if indices is not None:
304	msg = 'Only one of polygon and indices can be specified'
305	raise Exception, msg
306
307	msg = 'With polygon selected, set_quantity must provide '
308	msg += 'the keyword numeric and it must (currently) be '
309	msg += 'a constant.'
310	if numeric is None:
311	raise Exception, msg
312	else:
313	# Check that numeric is as constant
314	assert type(numeric) in [FloatType, IntType, LongType], msg
315
316
317	location = 'centroids'
318
319
320	points = self.domain.get_centroid_coordinates(absolute=True)
321	indices = inside_polygon(points, polygon)
322
323	self.set_values_from_constant(numeric,
324	location, indices, verbose)
325
326	self.extrapolate_first_order()
327	return
328
329
330
331
332
333
334	#General input checks
335	L = [numeric, quantity, function, geospatial_data, points, filename]
336	msg = 'Exactly one of the arguments '+\
337	'numeric, quantity, function, geospatial_data, points, '+\
338	'or filename must be present.'
339	assert L.count(None) == len(L)-1, msg
340
341
342	if location not in ['vertices', 'centroids', 'edges',
343	'unique vertices']:
344	msg = 'Invalid location: %s' %location
345	raise Exception, msg
346
347
348	msg = 'Indices must be a list or None'
349	assert type(indices) in [ListType, NoneType, ArrayType], msg
350
351
352	if not(points is None and values is None and data_georef is None):
353	from warnings import warn
354
355	msg = 'Using points, values or data_georef with set_quantity '
356	msg += 'is obsolete. Please use a Geospatial_data object instead.'
357	#warn(msg, DeprecationWarning)
358	raise Exception, msg
359
360
361
362	#Determine which 'set_values_from_...' to use
363
364	if numeric is not None:
365	if type(numeric) in [FloatType, IntType, LongType]:
366	self.set_values_from_constant(numeric,
367	location, indices, verbose)
368	elif type(numeric) in [ArrayType, ListType]:
369	self.set_values_from_array(numeric,
370	location, indices, verbose)
371	elif callable(numeric):
372	self.set_values_from_function(numeric,
373	location, indices, verbose)
374	elif isinstance(numeric, Quantity):
375	self.set_values_from_quantity(numeric,
376	location, indices, verbose)
377	elif isinstance(numeric, Geospatial_data):
378	self.set_values_from_geospatial_data(numeric,
379	alpha,
380	location, indices,
381	verbose = verbose,
382	use_cache = use_cache)
383	else:
384	msg = 'Illegal type for argument numeric: %s' %str(numeric)
385	raise msg
386
387	elif quantity is not None:
388	self.set_values_from_quantity(quantity,
389	location, indices, verbose)
390	elif function is not None:
391	msg = 'Argument function must be callable'
392	assert callable(function), msg
393	self.set_values_from_function(function,
394	location, indices, verbose)
395	elif geospatial_data is not None:
396	self.set_values_from_geospatial_data(geospatial_data,
397	alpha,
398	location, indices,
399	verbose = verbose,
400	use_cache = use_cache)
401	elif points is not None:
402	msg = 'The usage of points in set_values has been deprecated.' +\
403	'Please use the geospatial_data object instead.'
404	raise Exception, msg
405
406
407
408	elif filename is not None:
409	if hasattr(self.domain, 'points_file_block_line_size'):
410	max_read_lines = self.domain.points_file_block_line_size
411	else:
412	max_read_lines = default_block_line_size
413	self.set_values_from_file(filename, attribute_name, alpha,
414	location, indices,
415	verbose = verbose,
416	max_read_lines=max_read_lines,
417	use_cache = use_cache)
418	else:
419	raise Exception, 'This can\'t happen :-)'
420
421
422
423	# Update all locations in triangles
424	if location == 'vertices' or location == 'unique vertices':
425	# Intialise centroid and edge_values
426	self.interpolate()
427
428	if location == 'centroids':
429	# Extrapolate 1st order - to capture notion of area being specified
430	self.extrapolate_first_order()
431
432
433
434	#Specific functions for setting values
435	def set_values_from_constant(self, X,
436	location, indices, verbose):
437	"""Set quantity values from specified constant X
438	"""
439
440	# FIXME (Ole): Somehow indices refer to centroids
441	# rather than vertices as default. See unit test
442	# test_set_vertex_values_using_general_interface_with_subset(self):
443
444
445	if location == 'centroids':
446	if indices is None:
447	self.centroid_values[:] = X
448	else:
449	#Brute force
450	for i in indices:
451	self.centroid_values[i] = X
452
453	elif location == 'edges':
454	if indices is None:
455	self.edge_values[:] = X
456	else:
457	#Brute force
458	for i in indices:
459	self.edge_values[i] = X
460
461	elif location == 'unique vertices':
462	if indices is None:
463	self.edge_values[:] = X #FIXME (Ole): Shouldn't this be vertex_values?
464	else:
465
466	#Go through list of unique vertices
467	for unique_vert_id in indices:
468
469	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
470
471	#In case there are unused points
472	if len(triangles) == 0:
473	continue
474
475	#Go through all triangle, vertex pairs
476	#and set corresponding vertex value
477	for triangle_id, vertex_id in triangles:
478	self.vertex_values[triangle_id, vertex_id] = X
479
480	#Intialise centroid and edge_values
481	self.interpolate()
482	else:
483	if indices is None:
484	self.vertex_values[:] = X
485	else:
486	#Brute force
487	for i_vertex in indices:
488	self.vertex_values[i_vertex] = X
489
490
491
492
493	def set_values_from_array(self, values,
494	location='vertices',
495	indices=None,
496	verbose=False):
497	"""Set values for quantity
498
499	values: Numeric array
500	location: Where values are to be stored.
501	Permissible options are: vertices, edges, centroid, unique vertices
502	Default is 'vertices'
503
504	indices - if this action is carried out on a subset of
505	elements or unique vertices
506	The element/unique vertex indices are specified here.
507
508	In case of location == 'centroid' the dimension values must
509	be a list of a Numerical array of length N, N being the number
510	of elements.
511
512	Otherwise it must be of dimension Nx3
513
514	The values will be stored in elements following their
515	internal ordering.
516
517	If selected location is vertices, values for centroid and edges
518	will be assigned interpolated values.
519	In any other case, only values for the specified locations
520	will be assigned and the others will be left undefined.
521	"""
522
523	from Numeric import array, Float, Int, allclose
524
525	values = array(values).astype(Float)
526
527	if indices is not None:
528	indices = array(indices).astype(Int)
529	msg = 'Number of values must match number of indices:'
530	msg += 'You specified %d values and %d indices'\
531	%(values.shape[0], indices.shape[0])
532	assert values.shape[0] == indices.shape[0], msg
533
534	N = self.centroid_values.shape[0]
535
536	if location == 'centroids':
537	assert len(values.shape) == 1, 'Values array must be 1d'
538
539	if indices is None:
540	msg = 'Number of values must match number of elements'
541	assert values.shape[0] == N, msg
542
543	self.centroid_values = values
544	else:
545	msg = 'Number of values must match number of indices'
546	assert values.shape[0] == indices.shape[0], msg
547
548	#Brute force
549	for i in range(len(indices)):
550	self.centroid_values[indices[i]] = values[i]
551
552	elif location == 'edges':
553	# FIXME (Ole): No mention of indices here. However, I don't
554	# think we ever need to set values at edges anyway
555	assert len(values.shape) == 2, 'Values array must be 2d'
556
557	msg = 'Number of values must match number of elements'
558	assert values.shape[0] == N, msg
559
560	msg = 'Array must be N x 3'
561	assert values.shape[1] == 3, msg
562
563	self.edge_values = values
564
565	elif location == 'unique vertices':
566	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
567	'Values array must be 1d'
568
569	self.set_vertex_values(values.flat, indices=indices)
570
571	else:
572	# Location vertices
573	if len(values.shape) == 1:
574	self.set_vertex_values(values, indices=indices)
575
576	elif len(values.shape) == 2:
577	#Vertex values are given as a triplet for each triangle
578
579	msg = 'Array must be N x 3'
580	assert values.shape[1] == 3, msg
581
582	if indices is None:
583	self.vertex_values = values
584	else:
585	for element_index, value in map(None, indices, values):
586	self.vertex_values[element_index] = value
587	else:
588	msg = 'Values array must be 1d or 2d'
589	raise msg
590
591
592	def set_values_from_quantity(self, q,
593	location, indices, verbose):
594	"""Set quantity values from specified quantity instance q
595
596	Location is ignored - vertices will always be used here.
597	"""
598
599
600	A = q.vertex_values
601
602	from Numeric import allclose
603	msg = 'Quantities are defined on different meshes. '+\
604	'This might be a case for implementing interpolation '+\
605	'between different meshes.'
606	assert allclose(A.shape, self.vertex_values.shape), msg
607
608	self.set_values(A, location='vertices',
609	indices=indices,
610	verbose=verbose)
611
612
613	def set_values_from_function(self, f,
614	location='vertices',
615	indices=None,
616	verbose=False):
617	"""Set values for quantity using specified function
618
619	Input
620
621	f: x, y -> z Function where x, y and z are arrays
622	location: Where values are to be stored.
623	Permissible options are: vertices, centroid, edges,
624	unique vertices
625	Default is "vertices"
626	indices:
627
628
629	"""
630
631	#FIXME: Should check that function returns something sensible and
632	#raise a meaningfull exception if it returns None for example
633
634	#FIXME: Should supply absolute coordinates
635
636
637	# Compute the function values and call set_values again
638	if location == 'centroids':
639	if indices is None:
640	indices = range(len(self))
641
642	V = take(self.domain.get_centroid_coordinates(), indices)
643	self.set_values(f(V[:,0], V[:,1]),
644	location=location,
645	indices=indices)
646
647	elif location == 'vertices':
648
649	M = self.domain.number_of_triangles
650	V = self.domain.get_vertex_coordinates()
651
652	x = V[:,0]; y = V[:,1];
653	values = f(x, y)
654
655
656	# FIXME (Ole): This code should replace all the
657	# rest of this function and it would work, except
658	# one unit test in test_region fails.
659	# If that could be resolved this one will be
660	# more robust and simple.
661
662	#values = reshape(values, (M,3))
663	#self.set_values(values,
664	# location='vertices',
665	# indices=indices)
666
667
668	# This should be removed
669	if is_scalar(values):
670	# Function returned a constant value
671	self.set_values_from_constant(values,
672	location, indices, verbose)
673	return
674
675	# This should be removed
676	if indices is None:
677	for j in range(3):
678	self.vertex_values[:,j] = values[j::3]
679	else:
680	#Brute force
681	for i in indices:
682	for j in range(3):
683	self.vertex_values[i,j] = values[3*i+j]
684
685
686	else:
687	raise 'Not implemented: %s' %location
688
689
690
691	def set_values_from_geospatial_data(self, geospatial_data, alpha,
692	location, indices,
693	verbose = False,
694	use_cache = False):
695
696	#FIXME: Use this function for the time being. Later move code in here
697
698	points = geospatial_data.get_data_points(absolute = False)
699	values = geospatial_data.get_attributes()
700	data_georef = geospatial_data.get_geo_reference()
701
702
703
704	from anuga.coordinate_transforms.geo_reference import Geo_reference
705
706
707	points = ensure_numeric(points, Float)
708	values = ensure_numeric(values, Float)
709
710	if location != 'vertices':
711	msg = 'set_values_from_points is only defined for '+\
712	'location=\'vertices\''
713	raise ms
714
715	coordinates = self.domain.get_nodes()
716	triangles = self.domain.triangles #FIXME
717
718
719	#Take care of georeferencing
720	if data_georef is None:
721	data_georef = Geo_reference()
722
723
724	mesh_georef = self.domain.geo_reference
725
726
727	# Call fit_interpolate.fit function
728	args = (coordinates, triangles, points, values)
729	kwargs = {'data_origin': data_georef.get_origin(),
730	'mesh_origin': mesh_georef.get_origin(),
731	'alpha': alpha,
732	'verbose': verbose}
733
734	vertex_attributes = apply(fit_to_mesh,
735	args, kwargs)
736
737	# Call underlying method using array values
738	self.set_values_from_array(vertex_attributes,
739	location, indices, verbose)
740
741
742
743	def set_values_from_points(self, points, values, alpha,
744	location, indices,
745	data_georef = None,
746	verbose = False,
747	use_cache = False):
748	"""
749	Set quantity values from arbitray data points using
750	fit_interpolate.fit
751	"""
752
753	raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
754
755
756	def set_values_from_file(self, filename, attribute_name, alpha,
757	location, indices,
758	verbose = False,
759	use_cache = False,
760	max_read_lines=None):
761	"""Set quantity based on arbitrary points in a points file
762	using attribute_name selects name of attribute
763	present in file.
764	If attribute_name is not specified, use first available attribute
765	as defined in geospatial_data.
766	"""
767
768	from types import StringType
769	msg = 'Filename must be a text string'
770	assert type(filename) == StringType, msg
771
772
773	if location != 'vertices':
774	msg = 'set_values_from_file is only defined for '+\
775	'location=\'vertices\''
776	raise msg
777
778	coordinates = self.domain.get_nodes(absolute=True)
779	triangles = self.domain.triangles #FIXME
780
781	vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
782	alpha=alpha,
783	attribute_name=attribute_name,
784	use_cache=use_cache,
785	verbose=verbose,
786	max_read_lines=max_read_lines)
787
788	# Call underlying method using array values
789	self.set_values_from_array(vertex_attributes,
790	location, indices, verbose)
791
792
793	def get_extremum_index(self, mode=None, indices=None):
794	"""Return index for maximum or minimum value of quantity (on centroids)
795
796	Optional arguments:
797	mode is either 'max'(default) or 'min'.
798	indices is the set of element ids that the operation applies to.
799
800	Usage:
801	i = get_extreme_index()
802
803	Notes:
804	We do not seek the extremum at vertices as each vertex can
805	have multiple values - one for each triangle sharing it.
806
807	If there are multiple cells with same maximum value, the
808	first cell encountered in the triangle array is returned.
809	"""
810
811	V = self.get_values(location='centroids', indices=indices)
812
813	# Always return absolute indices
814	if mode is None or mode == 'max':
815	i = argmax(V)
816	elif mode == 'min':
817	i = argmin(V)
818
819
820	if indices is None:
821	return i
822	else:
823	return indices[i]
824
825
826	def get_maximum_index(self, indices=None):
827	"""See get extreme index for details
828	"""
829
830	return self.get_extremum_index(mode='max',
831	indices=indices)
832
833
834
835	def get_maximum_value(self, indices=None):
836	"""Return maximum value of quantity (on centroids)
837
838	Optional argument:
839	indices is the set of element ids that the operation applies to.
840
841	Usage:
842	v = get_maximum_value()
843
844	Note, we do not seek the maximum at vertices as each vertex can
845	have multiple values - one for each triangle sharing it
846	"""
847
848
849	i = self.get_maximum_index(indices)
850	V = self.get_values(location='centroids') #, indices=indices)
851
852	return V[i]
853
854
855	def get_maximum_location(self, indices=None):
856	"""Return location of maximum value of quantity (on centroids)
857
858	Optional argument:
859	indices is the set of element ids that the operation applies to.
860
861	Usage:
862	x, y = get_maximum_location()
863
864
865	Notes:
866	We do not seek the maximum at vertices as each vertex can
867	have multiple values - one for each triangle sharing it.
868
869	If there are multiple cells with same maximum value, the
870	first cell encountered in the triangle array is returned.
871	"""
872
873	i = self.get_maximum_index(indices)
874	x, y = self.domain.get_centroid_coordinates()[i]
875
876	return x, y
877
878
879	def get_minimum_index(self, indices=None):
880	"""See get extreme index for details
881	"""
882
883	return self.get_extremum_index(mode='min',
884	indices=indices)
885
886
887	def get_minimum_value(self, indices=None):
888	"""Return minimum value of quantity (on centroids)
889
890	Optional argument:
891	indices is the set of element ids that the operation applies to.
892
893	Usage:
894	v = get_minimum_value()
895
896	See get_maximum_value for more details.
897	"""
898
899
900	i = self.get_minimum_index(indices)
901	V = self.get_values(location='centroids')
902
903	return V[i]
904
905
906	def get_minimum_location(self, indices=None):
907	"""Return location of minimum value of quantity (on centroids)
908
909	Optional argument:
910	indices is the set of element ids that the operation applies to.
911
912	Usage:
913	x, y = get_minimum_location()
914
915
916	Notes:
917	We do not seek the maximum at vertices as each vertex can
918	have multiple values - one for each triangle sharing it.
919
920	If there are multiple cells with same maximum value, the
921	first cell encountered in the triangle array is returned.
922	"""
923
924	i = self.get_minimum_index(indices)
925	x, y = self.domain.get_centroid_coordinates()[i]
926
927	return x, y
928
929
930
931
932	def get_interpolated_values(self, interpolation_points):
933
934	# Interpolation object based on internal (discontinuous triangles)
935	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
936	smooth=False)
937	# FIXME: This concat should roll into get_vertex_values
938	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
939	axis=1)
940
941	can_reuse = False
942	if hasattr(self, 'interpolation_object'):
943	# Reuse to save time
944	I = self.interpolation_object
945
946	if allclose(interpolation_points, I._point_coordinates):
947	can_reuse = True
948
949
950	if can_reuse is True:
951	# Use absence of points to indicate reuse in I.interpolate
952	result = I.interpolate(vertex_values)
953	else:
954	from anuga.fit_interpolate.interpolate import Interpolate
955
956	# Create interpolation object with matrix
957	I = Interpolate(vertex_coordinates, triangles)
958	self.interpolation_object = I
959
960	# Call interpolate with points the first time
961	interpolation_points = ensure_numeric(interpolation_points, Float)
962	result = I.interpolate(vertex_values, interpolation_points)
963
964	return result
965
966
967	def get_values(self, interpolation_points=None,
968	location='vertices',
969	indices = None):
970	"""get values for quantity
971
972	return X, Compatible list, Numeric array (see below)
973	interpolation_points: List of x, y coordinates where value is
974	sought (using interpolation). If points are given, values of
975	location and indices are ignored
976
977	location: Where values are to be stored.
978	Permissible options are: vertices, edges, centroids
979	and unique vertices. Default is 'vertices'
980
981
982	The returned values with be a list the length of indices
983	(N if indices = None).
984
985	In case of location == 'centroids' the dimension of returned
986	values will be a list or a Numerical array of length N, N being
987	the number of elements.
988
989	In case of location == 'vertices' or 'edges' the dimension of
990	returned values will be of dimension Nx3
991
992	In case of location == 'unique vertices' the average value at
993	each vertex will be returned and the dimension of returned values
994	will be a 1d array of length "number of vertices"
995
996	Indices is the set of element ids that the operation applies to.
997
998	The values will be stored in elements following their
999	internal ordering.
1000	"""
1001	from Numeric import take
1002
1003	# FIXME (Ole): I reckon we should have the option of passing a
1004	# polygon into get_values. The question becomes how
1005	# resulting values should be ordered.
1006
1007	if interpolation_points is not None:
1008	return self.get_interpolated_values(interpolation_points)
1009
1010
1011
1012	if location not in ['vertices', 'centroids', 'edges',
1013	'unique vertices']:
1014	msg = 'Invalid location: %s' %location
1015	raise msg
1016
1017	import types, Numeric
1018	assert type(indices) in [types.ListType, types.NoneType,
1019	Numeric.ArrayType],\
1020	'Indices must be a list or None'
1021
1022	if location == 'centroids':
1023	if (indices == None):
1024	indices = range(len(self))
1025	return take(self.centroid_values,indices)
1026	elif location == 'edges':
1027	if (indices == None):
1028	indices = range(len(self))
1029	return take(self.edge_values,indices)
1030	elif location == 'unique vertices':
1031	if (indices == None):
1032	indices=range(self.domain.number_of_nodes)
1033	vert_values = []
1034	#Go through list of unique vertices
1035	for unique_vert_id in indices:
1036	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1037
1038	#In case there are unused points
1039	if len(triangles) == 0:
1040	msg = 'Unique vertex not associated with triangles'
1041	raise msg
1042
1043	# Go through all triangle, vertex pairs
1044	# Average the values
1045
1046	# FIXME (Ole): Should we merge this with get_vertex_values
1047	sum = 0
1048	for triangle_id, vertex_id in triangles:
1049	sum += self.vertex_values[triangle_id, vertex_id]
1050	vert_values.append(sum/len(triangles))
1051	return Numeric.array(vert_values)
1052	else:
1053	if (indices is None):
1054	indices = range(len(self))
1055	return take(self.vertex_values, indices)
1056
1057
1058
1059	def set_vertex_values(self, A, indices = None):
1060	"""Set vertex values for all unique vertices based on input array A
1061	which has one entry per unique vertex, i.e.
1062	one value for each row in array self.domain.nodes.
1063
1064	indices is the list of vertex_id's that will be set.
1065
1066	This function is used by set_values_from_array
1067	"""
1068
1069	from Numeric import array, Float
1070
1071	#Assert that A can be converted to a Numeric array of appropriate dim
1072	A = array(A, Float)
1073
1074	#print 'SHAPE A', A.shape
1075	assert len(A.shape) == 1
1076
1077	if indices is None:
1078	assert A.shape[0] == self.domain.get_nodes().shape[0]
1079	vertex_list = range(A.shape[0])
1080	else:
1081	assert A.shape[0] == len(indices)
1082	vertex_list = indices
1083
1084	#Go through list of unique vertices
1085
1086	for i_index, unique_vert_id in enumerate(vertex_list):
1087
1088
1089	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1090
1091	#In case there are unused points
1092	if len(triangles) == 0: continue
1093
1094	#Go through all triangle, vertex pairs
1095	#touching vertex unique_vert_id and set corresponding vertex value
1096	for triangle_id, vertex_id in triangles:
1097	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1098
1099	#Intialise centroid and edge_values
1100	self.interpolate()
1101
1102
1103	def smooth_vertex_values(self, value_array='field_values',
1104	precision = None):
1105	""" Smooths field_values or conserved_quantities data.
1106	TODO: be able to smooth individual fields
1107	NOTE: This function does not have a test.
1108	FIXME: NOT DONE - do we need it?
1109	FIXME: this function isn't called by anything.
1110	Maybe it should be removed..-DSG
1111	"""
1112
1113	from Numeric import concatenate, zeros, Float, Int, array, reshape
1114
1115
1116	A,V = self.get_vertex_values(xy=False,
1117	value_array=value_array,
1118	smooth = True,
1119	precision = precision)
1120
1121	#Set some field values
1122	for volume in self:
1123	for i,v in enumerate(volume.vertices):
1124	if value_array == 'field_values':
1125	volume.set_field_values('vertex', i, A[v,:])
1126	elif value_array == 'conserved_quantities':
1127	volume.set_conserved_quantities('vertex', i, A[v,:])
1128
1129	if value_array == 'field_values':
1130	self.precompute()
1131	elif value_array == 'conserved_quantities':
1132	Volume.interpolate_conserved_quantities()
1133
1134
1135	# Methods for outputting model results
1136	def get_vertex_values(self,
1137	xy=True,
1138	smooth=None,
1139	precision=None):
1140	"""Return vertex values like an OBJ format i.e. one value per node.
1141
1142	The vertex values are returned as one sequence in the 1D float array A.
1143	If requested the coordinates will be returned in 1D arrays X and Y.
1144
1145	The connectivity is represented as an integer array, V, of dimension
1146	Mx3, where M is the number of triangles. Each row has three indices
1147	defining the triangle and they correspond to elements in the arrays
1148	X, Y and A.
1149
1150	if smooth is True, vertex values corresponding to one common
1151	coordinate set will be smoothed by taking the average of vertex values for each node.
1152	In this case vertex coordinates will be
1153	de-duplicated corresponding to the original nodes as obtained from
1154	the method general_mesh.get_nodes()
1155
1156	If no smoothings is required, vertex coordinates and values will
1157	be aggregated as a concatenation of values at
1158	vertices 0, vertices 1 and vertices 2. This corresponds to
1159	the node coordinates obtained from the method
1160	general_mesh.get_vertex_coordinates()
1161
1162
1163	Calling convention
1164	if xy is True:
1165	X,Y,A,V = get_vertex_values
1166	else:
1167	A,V = get_vertex_values
1168
1169	"""
1170
1171	from Numeric import concatenate, zeros, Float, Int, array, reshape
1172
1173
1174	if smooth is None:
1175	# Take default from domain
1176	smooth = self.domain.smooth
1177
1178	if precision is None:
1179	precision = Float
1180
1181
1182	if smooth is True:
1183	# Ensure continuous vertex values by averaging
1184	# values at each node
1185
1186	V = self.domain.get_triangles()
1187	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1188	A = zeros(N, Float)
1189	points = self.domain.get_nodes()
1190
1191	if 1:
1192	# Fast C version
1193	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1194	ensure_numeric(self.domain.number_of_triangles_per_node),
1195	ensure_numeric(self.vertex_values),
1196	A)
1197	A = A.astype(precision)
1198	else:
1199
1200	# Slow Python version
1201
1202	current_node = 0
1203	k = 0 # Track triangles touching on node
1204	total = 0.0
1205	for index in self.domain.vertex_value_indices:
1206	if current_node == N:
1207	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1208	raise msg
1209
1210
1211
1212	k += 1
1213
1214	volume_id = index / 3
1215	vertex_id = index % 3
1216
1217	#assert V[volume_id, vertex_id] == current_node
1218
1219	v = self.vertex_values[volume_id, vertex_id]
1220	total += v
1221
1222	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1223	if self.domain.number_of_triangles_per_node[current_node] == k:
1224	A[current_node] = total/k
1225
1226
1227	# Move on to next node
1228	total = 0.0
1229	k = 0
1230	current_node += 1
1231
1232
1233
1234	else:
1235	# Allow discontinuous vertex values
1236	V = self.domain.get_disconnected_triangles()
1237	points = self.domain.get_vertex_coordinates()
1238	A = self.vertex_values.flat.astype(precision)
1239
1240
1241	# Return
1242	if xy is True:
1243	X = points[:,0].astype(precision)
1244	Y = points[:,1].astype(precision)
1245
1246	return X, Y, A, V
1247	else:
1248	return A, V
1249
1250
1251
1252	def extrapolate_first_order(self):
1253	"""Extrapolate conserved quantities from centroid to
1254	vertices for each volume using
1255	first order scheme.
1256	"""
1257
1258	qc = self.centroid_values
1259	qv = self.vertex_values
1260
1261	for i in range(3):
1262	qv[:,i] = qc
1263
1264
1265	def get_integral(self):
1266	"""Compute the integral of quantity across entire domain
1267	"""
1268	integral = 0
1269	for k in range(len(self.domain)):
1270	area = self.domain.areas[k]
1271	qc = self.centroid_values[k]
1272	integral += qc*area
1273
1274	return integral
1275
1276
1277
1278
1279	class Conserved_quantity(Quantity):
1280	"""Class conserved quantity adds to Quantity:
1281
1282	boundary values, storage and method for updating, and
1283	methods for (second order) extrapolation from centroid to vertices inluding
1284	gradients and limiters
1285	"""
1286
1287	def __init__(self, domain, vertex_values=None):
1288	Quantity.__init__(self, domain, vertex_values)
1289
1290	from Numeric import zeros, Float
1291
1292	#Allocate space for boundary values
1293	L = len(domain.boundary)
1294	self.boundary_values = zeros(L, Float)
1295
1296	#Allocate space for updates of conserved quantities by
1297	#flux calculations and forcing functions
1298
1299	N = len(domain) # number_of_triangles
1300	self.explicit_update = zeros(N, Float )
1301	self.semi_implicit_update = zeros(N, Float )
1302	self.centroid_backup_values = zeros(N, Float)
1303
1304
1305
1306
1307	def update(self, timestep):
1308	#Call correct module function
1309	#(either from this module or C-extension)
1310	return update(self, timestep)
1311
1312
1313	def compute_gradients(self):
1314	#Call correct module function
1315	#(either from this module or C-extension)
1316	return compute_gradients(self)
1317
1318	def limit(self):
1319	#Call correct module function
1320	#(either from this module or C-extension)
1321	limit(self)
1322
1323	def extrapolate_second_order(self):
1324	#Call correct module function
1325	#(either from this module or C-extension)
1326	extrapolate_second_order(self)
1327
1328	def backup_centroid_values(self):
1329	#Call correct module function
1330	#(either from this module or C-extension)
1331	backup_centroid_values(self)
1332
1333	def saxpy_centroid_values(self,a,b):
1334	#Call correct module function
1335	#(either from this module or C-extension)
1336	saxpy_centroid_values(self,a,b)
1337
1338
1339	def update(quantity, timestep):
1340	"""Update centroid values based on values stored in
1341	explicit_update and semi_implicit_update as well as given timestep
1342
1343	Function implementing forcing terms must take on argument
1344	which is the domain and they must update either explicit
1345	or implicit updates, e,g,:
1346
1347	def gravity(domain):
1348	....
1349	domain.quantities['xmomentum'].explicit_update = ...
1350	domain.quantities['ymomentum'].explicit_update = ...
1351
1352
1353
1354	Explicit terms must have the form
1355
1356	G(q, t)
1357
1358	and explicit scheme is
1359
1360	q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
1361
1362
1363	Semi implicit forcing terms are assumed to have the form
1364
1365	G(q, t) = H(q, t) q
1366
1367	and the semi implicit scheme will then be
1368
1369	q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
1370
1371
1372	"""
1373
1374	from Numeric import sum, equal, ones, exp, Float
1375
1376	N = quantity.centroid_values.shape[0]
1377
1378
1379	# Divide H by conserved quantity to obtain G (see docstring above)
1380
1381
1382	for k in range(N):
1383	x = quantity.centroid_values[k]
1384	if x == 0.0:
1385	#FIXME: Is this right
1386	quantity.semi_implicit_update[k] = 0.0
1387	else:
1388	quantity.semi_implicit_update[k] /= x
1389
1390
1391	# Semi implicit updates
1392	denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1393
1394	if sum(less(denominator, 1.0)) > 0.0:
1395	msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
1396	raise msg
1397
1398	if sum(equal(denominator, 0.0)) > 0.0:
1399	msg = 'Zero division in semi implicit update. Call Stephen :-)'
1400	raise msg
1401	else:
1402	#Update conserved_quantities from semi implicit updates
1403	quantity.centroid_values /= denominator
1404
1405	# quantity.centroid_values = exp(timestepquantity.semi_implicit_update)quantity.centroid_values
1406
1407	#Explicit updates
1408	quantity.centroid_values += timestep*quantity.explicit_update
1409
1410	def interpolate_from_vertices_to_edges(quantity):
1411	"""Compute edge values from vertex values using linear interpolation
1412	"""
1413
1414	for k in range(quantity.vertex_values.shape[0]):
1415	q0 = quantity.vertex_values[k, 0]
1416	q1 = quantity.vertex_values[k, 1]
1417	q2 = quantity.vertex_values[k, 2]
1418
1419	quantity.edge_values[k, 0] = 0.5*(q1+q2)
1420	quantity.edge_values[k, 1] = 0.5*(q0+q2)
1421	quantity.edge_values[k, 2] = 0.5*(q0+q1)
1422
1423
1424
1425	def backup_centroid_values(quantity):
1426	"""Copy centroid values to backup array"""
1427
1428	qc = quantity.centroid_values
1429	qb = quantity.centroid_backup_values
1430
1431	#Check each triangle
1432	for k in range(len(quantity.domain)):
1433	qb[k] = qc[k]
1434
1435
1436	def saxpy_centroid_values(quantity,a,b):
1437	"""saxpy operation between centroid value and backup"""
1438
1439	qc = quantity.centroid_values
1440	qb = quantity.centroid_backup_values
1441
1442	#Check each triangle
1443	for k in range(len(quantity.domain)):
1444	qc[k] = aqc[k]+bqb[k]
1445
1446
1447	def compute_gradients(quantity):
1448	"""Compute gradients of triangle surfaces defined by centroids of
1449	neighbouring volumes.
1450	If one edge is on the boundary, use own centroid as neighbour centroid.
1451	If two or more are on the boundary, fall back to first order scheme.
1452	"""
1453
1454	from Numeric import zeros, Float
1455	from utilitites.numerical_tools import gradient
1456
1457	centroid_coordinates = quantity.domain.centroid_coordinates
1458	surrogate_neighbours = quantity.domain.surrogate_neighbours
1459	centroid_values = quantity.centroid_values
1460	number_of_boundaries = quantity.domain.number_of_boundaries
1461
1462	N = centroid_values.shape[0]
1463
1464	a = zeros(N, Float)
1465	b = zeros(N, Float)
1466
1467	for k in range(N):
1468	if number_of_boundaries[k] < 2:
1469	#Two or three true neighbours
1470
1471	#Get indices of neighbours (or self when used as surrogate)
1472	k0, k1, k2 = surrogate_neighbours[k,:]
1473
1474	#Get data
1475	q0 = centroid_values[k0]
1476	q1 = centroid_values[k1]
1477	q2 = centroid_values[k2]
1478
1479	x0, y0 = centroid_coordinates[k0] #V0 centroid
1480	x1, y1 = centroid_coordinates[k1] #V1 centroid
1481	x2, y2 = centroid_coordinates[k2] #V2 centroid
1482
1483	#Gradient
1484	a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1485
1486	elif number_of_boundaries[k] == 2:
1487	#One true neighbour
1488
1489	#Get index of the one neighbour
1490	for k0 in surrogate_neighbours[k,:]:
1491	if k0 != k: break
1492	assert k0 != k
1493
1494	k1 = k #self
1495
1496	#Get data
1497	q0 = centroid_values[k0]
1498	q1 = centroid_values[k1]
1499
1500	x0, y0 = centroid_coordinates[k0] #V0 centroid
1501	x1, y1 = centroid_coordinates[k1] #V1 centroid
1502
1503	#Gradient
1504	a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1505	else:
1506	#No true neighbours -
1507	#Fall back to first order scheme
1508	pass
1509
1510
1511	return a, b
1512
1513
1514
1515	def limit(quantity):
1516	"""Limit slopes for each volume to eliminate artificial variance
1517	introduced by e.g. second order extrapolator
1518
1519	This is an unsophisticated limiter as it does not take into
1520	account dependencies among quantities.
1521
1522	precondition:
1523	vertex values are estimated from gradient
1524	postcondition:
1525	vertex values are updated
1526	"""
1527
1528	from Numeric import zeros, Float
1529
1530	N = quantity.domain.number_of_nodes
1531
1532	beta_w = quantity.domain.beta_w
1533
1534	qc = quantity.centroid_values
1535	qv = quantity.vertex_values
1536
1537	#Find min and max of this and neighbour's centroid values
1538	qmax = zeros(qc.shape, Float)
1539	qmin = zeros(qc.shape, Float)
1540
1541	for k in range(N):
1542	qmax[k] = qc[k]
1543	qmin[k] = qc[k]
1544	for i in range(3):
1545	n = quantity.domain.neighbours[k,i]
1546	if n >= 0:
1547	qn = qc[n] #Neighbour's centroid value
1548
1549	qmin[k] = min(qmin[k], qn)
1550	qmax[k] = max(qmax[k], qn)
1551	qmax[k] = min(qmax[k], 2.0*qc[k])
1552	qmin[k] = max(qmin[k], 0.5*qc[k])
1553
1554
1555	#Diffences between centroids and maxima/minima
1556	dqmax = qmax - qc
1557	dqmin = qmin - qc
1558
1559	#Deltas between vertex and centroid values
1560	dq = zeros(qv.shape, Float)
1561	for i in range(3):
1562	dq[:,i] = qv[:,i] - qc
1563
1564	#Phi limiter
1565	for k in range(N):
1566
1567	#Find the gradient limiter (phi) across vertices
1568	phi = 1.0
1569	for i in range(3):
1570	r = 1.0
1571	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1572	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1573
1574	phi = min( min(r*beta_w, 1), phi )
1575
1576	#Then update using phi limiter
1577	for i in range(3):
1578	qv[k,i] = qc[k] + phi*dq[k,i]
1579
1580
1581
1582	from anuga.utilities import compile
1583	if compile.can_use_C_extension('quantity_ext.c'):
1584	#Replace python version with c implementations
1585
1586	from quantity_ext import average_vertex_values, backup_centroid_values, \
1587	saxpy_centroid_values
1588
1589	from quantity_ext import compute_gradients, limit,\
1590	extrapolate_second_order, interpolate_from_vertices_to_edges, update

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

Download in other formats: