Context Navigation

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

Last change on this file since 4808 was 4808, checked in by duncan, 16 years ago
bug fix
File size: 46.9 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
153	# FIXME (Ole): Maybe this function
154	# should move to the C-interface?
155	# However, it isn't called by validate_all.py, so it
156	# may not be that important to optimise it?
157
158	N = self.vertex_values.shape[0]
159	for i in range(N):
160	v0 = self.vertex_values[i, 0]
161	v1 = self.vertex_values[i, 1]
162	v2 = self.vertex_values[i, 2]
163
164	self.centroid_values[i] = (v0 + v1 + v2)/3
165
166	self.interpolate_from_vertices_to_edges()
167
168
169	def interpolate_from_vertices_to_edges(self):
170	#Call correct module function
171	#(either from this module or C-extension)
172	interpolate_from_vertices_to_edges(self)
173
174
175
176
177	# New leaner interface to setting values
178	def set_values(self,
179	numeric = None, # List, numeric array or constant
180	quantity = None, # Another quantity
181	function = None, # Callable object: f(x,y)
182	geospatial_data = None, # Arbitrary dataset
183	points = None, values = None, data_georef = None, # Obsoleted by use of geo_spatial object
184	filename = None, attribute_name = None, #Input from file
185	alpha = None,
186	location = 'vertices',
187	polygon = None,
188	indices = None,
189	verbose = False,
190	use_cache = False):
191
192	"""Set values for quantity based on different sources.
193
194	numeric:
195	Compatible list, Numeric array (see below) or constant.
196	If callable it will treated as a function (see below)
197	If instance of another Quantity it will be treated as such.
198	If geo_spatial object it will be treated as such
199
200	quantity:
201	Another quantity (compatible quantity, e.g. obtained as a
202	linear combination of quantities)
203
204	function:
205	Any callable object that takes two 1d arrays x and y
206	each of length N and returns an array also of length N.
207	The function will be evaluated at points determined by
208	location and indices in the underlying mesh.
209
210	geospatial_data:
211	Arbitrary geo spatial dataset in the form of the class
212	Geospatial_data. Mesh points are populated using
213	fit_interpolate.fit fitting
214
215	points:
216	Nx2 array of data points for use with fit_interpolate.fit
217	If points are present, an N array of attribute
218	values corresponding to
219	each data point must be present.
220	(Obsoleted by geospatial_data)
221
222	values:
223	If points is specified, values is an array of length N containing
224	attribute values for each point.
225	(Obsoleted by geospatial_data)
226
227	data_georef:
228	If points is specified, geo_reference applies to each point.
229	(Obsoleted by geospatial_data)
230
231	filename:
232	Name of a points file containing data points and attributes for
233	use with fit_interpolate.fit.
234
235	attribute_name:
236	If specified, any array matching that name
237	will be used. from file or geospatial_data.
238	Otherwise a default will be used.
239
240	alpha:
241	Smoothing parameter to be used with fit_interpolate.fit.
242	See module fit_interpolate.fit for further details about alpha.
243	Alpha will only be used with points, values or filename.
244	Otherwise it will be ignored.
245
246
247	location: Where values are to be stored.
248	Permissible options are: vertices, edges, centroids
249	Default is 'vertices'
250
251	In case of location == 'centroids' the dimension values must
252	be a list of a Numerical array of length N,
253	N being the number of elements.
254	Otherwise it must be of dimension Nx3
255
256
257	The values will be stored in elements following their
258	internal ordering.
259
260	If location is not 'unique vertices' Indices is the
261	set of element ids that the operation applies to.
262	If location is 'unique vertices' Indices is the set
263	of vertex ids that the operation applies to.
264
265	If selected location is vertices, values for
266	centroid and edges will be assigned interpolated
267	values. In any other case, only values for the
268	specified locations will be assigned and the others
269	will be left undefined.
270
271
272	polygon: Restrict update of quantity to locations that fall
273	inside polygon. Polygon works by selecting indices
274	and calling set_values recursively.
275	Polygon mode has only been implemented for
276	constant values so far.
277
278	indices: Restrict update of quantity to locations that are
279	identified by indices (e.g. node ids if location
280	is 'vertices')
281
282	verbose: True means that output to stdout is generated
283
284	use_cache: True means that caching of intermediate results is
285	attempted for fit_interpolate.fit.
286
287
288
289
290	Exactly one of the arguments
291	numeric, quantity, function, points, filename
292	must be present.
293	"""
294
295	from anuga.geospatial_data.geospatial_data import Geospatial_data
296	from types import FloatType, IntType, LongType, ListType, NoneType
297	from Numeric import ArrayType
298
299
300	# Treat special case: Polygon situation
301	# Location will be ignored and set to 'centroids'
302	# FIXME (Ole): This needs to be generalised and
303	# perhaps the notion of location and indices simplified
304
305	if polygon is not None:
306	if indices is not None:
307	msg = 'Only one of polygon and indices can be specified'
308	raise Exception, msg
309
310	msg = 'With polygon selected, set_quantity must provide '
311	msg += 'the keyword numeric and it must (currently) be '
312	msg += 'a constant.'
313	if numeric is None:
314	raise Exception, msg
315	else:
316	# Check that numeric is as constant
317	assert type(numeric) in [FloatType, IntType, LongType], msg
318
319
320	location = 'centroids'
321
322
323	points = self.domain.get_centroid_coordinates(absolute=True)
324	indices = inside_polygon(points, polygon)
325
326	self.set_values_from_constant(numeric,
327	location, indices, verbose)
328
329	self.extrapolate_first_order()
330	return
331
332
333
334
335
336
337	# General input checks
338	L = [numeric, quantity, function, geospatial_data, points, filename]
339	msg = 'Exactly one of the arguments '+\
340	'numeric, quantity, function, geospatial_data, points, '+\
341	'or filename must be present.'
342	assert L.count(None) == len(L)-1, msg
343
344
345	if location not in ['vertices', 'centroids', 'edges',
346	'unique vertices']:
347	msg = 'Invalid location: %s' %location
348	raise Exception, msg
349
350
351	msg = 'Indices must be a list or None'
352	assert type(indices) in [ListType, NoneType, ArrayType], msg
353
354
355	if not(points is None and values is None and data_georef is None):
356	from warnings import warn
357
358	msg = 'Using points, values or data_georef with set_quantity '
359	msg += 'is obsolete. Please use a Geospatial_data object instead.'
360	#warn(msg, DeprecationWarning)
361	raise Exception, msg
362
363
364
365	# Determine which 'set_values_from_...' to use
366
367	if numeric is not None:
368	if type(numeric) in [FloatType, IntType, LongType]:
369	self.set_values_from_constant(numeric,
370	location, indices, verbose)
371	elif type(numeric) in [ArrayType, ListType]:
372	self.set_values_from_array(numeric,
373	location, indices, verbose)
374	elif callable(numeric):
375	self.set_values_from_function(numeric,
376	location, indices, verbose)
377	elif isinstance(numeric, Quantity):
378	self.set_values_from_quantity(numeric,
379	location, indices, verbose)
380	elif isinstance(numeric, Geospatial_data):
381	self.set_values_from_geospatial_data(numeric,
382	alpha,
383	location, indices,
384	verbose = verbose,
385	use_cache = use_cache)
386	else:
387	msg = 'Illegal type for argument numeric: %s' %str(numeric)
388	raise msg
389
390	elif quantity is not None:
391	self.set_values_from_quantity(quantity,
392	location, indices, verbose)
393	elif function is not None:
394	msg = 'Argument function must be callable'
395	assert callable(function), msg
396	self.set_values_from_function(function,
397	location, indices, verbose)
398	elif geospatial_data is not None:
399	self.set_values_from_geospatial_data(geospatial_data,
400	alpha,
401	location, indices,
402	verbose = verbose,
403	use_cache = use_cache)
404	elif points is not None:
405	msg = 'The usage of points in set_values has been deprecated.' +\
406	'Please use the geospatial_data object instead.'
407	raise Exception, msg
408
409
410
411	elif filename is not None:
412	if hasattr(self.domain, 'points_file_block_line_size'):
413	max_read_lines = self.domain.points_file_block_line_size
414	else:
415	max_read_lines = default_block_line_size
416	self.set_values_from_file(filename, attribute_name, alpha,
417	location, indices,
418	verbose = verbose,
419	max_read_lines=max_read_lines,
420	use_cache = use_cache)
421	else:
422	raise Exception, 'This can\'t happen :-)'
423
424
425
426	# Update all locations in triangles
427	if location == 'vertices' or location == 'unique vertices':
428	# Intialise centroid and edge_values
429	self.interpolate()
430
431	if location == 'centroids':
432	# Extrapolate 1st order - to capture notion of area being specified
433	self.extrapolate_first_order()
434
435
436
437	#Specific functions for setting values
438	def set_values_from_constant(self, X,
439	location, indices, verbose):
440	"""Set quantity values from specified constant X
441	"""
442
443	# FIXME (Ole): Somehow indices refer to centroids
444	# rather than vertices as default. See unit test
445	# test_set_vertex_values_using_general_interface_with_subset(self):
446
447
448	if location == 'centroids':
449	if indices is None:
450	self.centroid_values[:] = X
451	else:
452	#Brute force
453	for i in indices:
454	self.centroid_values[i] = X
455
456	elif location == 'edges':
457	if indices is None:
458	self.edge_values[:] = X
459	else:
460	#Brute force
461	for i in indices:
462	self.edge_values[i] = X
463
464	elif location == 'unique vertices':
465	if indices is None:
466	self.edge_values[:] = X #FIXME (Ole): Shouldn't this be vertex_values?
467	else:
468
469	#Go through list of unique vertices
470	for unique_vert_id in indices:
471
472	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
473
474	#In case there are unused points
475	if len(triangles) == 0:
476	continue
477
478	#Go through all triangle, vertex pairs
479	#and set corresponding vertex value
480	for triangle_id, vertex_id in triangles:
481	self.vertex_values[triangle_id, vertex_id] = X
482
483	#Intialise centroid and edge_values
484	self.interpolate()
485	else:
486	if indices is None:
487	self.vertex_values[:] = X
488	else:
489	#Brute force
490	for i_vertex in indices:
491	self.vertex_values[i_vertex] = X
492
493
494
495
496	def set_values_from_array(self, values,
497	location='vertices',
498	indices=None,
499	verbose=False):
500	"""Set values for quantity
501
502	values: Numeric array
503	location: Where values are to be stored.
504	Permissible options are: vertices, edges, centroid, unique vertices
505	Default is 'vertices'
506
507	indices - if this action is carried out on a subset of
508	elements or unique vertices
509	The element/unique vertex indices are specified here.
510
511	In case of location == 'centroid' the dimension values must
512	be a list of a Numerical array of length N, N being the number
513	of elements.
514
515	Otherwise it must be of dimension Nx3
516
517	The values will be stored in elements following their
518	internal ordering.
519
520	If selected location is vertices, values for centroid and edges
521	will be assigned interpolated values.
522	In any other case, only values for the specified locations
523	will be assigned and the others will be left undefined.
524	"""
525
526	from Numeric import array, Float, Int, allclose
527
528	values = array(values).astype(Float)
529
530	if indices is not None:
531	indices = array(indices).astype(Int)
532	msg = 'Number of values must match number of indices:'
533	msg += 'You specified %d values and %d indices'\
534	%(values.shape[0], indices.shape[0])
535	assert values.shape[0] == indices.shape[0], msg
536
537	N = self.centroid_values.shape[0]
538
539	if location == 'centroids':
540	assert len(values.shape) == 1, 'Values array must be 1d'
541
542	if indices is None:
543	msg = 'Number of values must match number of elements'
544	assert values.shape[0] == N, msg
545
546	self.centroid_values = values
547	else:
548	msg = 'Number of values must match number of indices'
549	assert values.shape[0] == indices.shape[0], msg
550
551	#Brute force
552	for i in range(len(indices)):
553	self.centroid_values[indices[i]] = values[i]
554
555	elif location == 'edges':
556	# FIXME (Ole): No mention of indices here. However, I don't
557	# think we ever need to set values at edges anyway
558	assert len(values.shape) == 2, 'Values array must be 2d'
559
560	msg = 'Number of values must match number of elements'
561	assert values.shape[0] == N, msg
562
563	msg = 'Array must be N x 3'
564	assert values.shape[1] == 3, msg
565
566	self.edge_values = values
567
568	elif location == 'unique vertices':
569	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
570	'Values array must be 1d'
571
572	self.set_vertex_values(values.flat, indices=indices)
573
574	else:
575	# Location vertices
576	if len(values.shape) == 1:
577	self.set_vertex_values(values, indices=indices)
578
579	elif len(values.shape) == 2:
580	#Vertex values are given as a triplet for each triangle
581
582	msg = 'Array must be N x 3'
583	assert values.shape[1] == 3, msg
584
585	if indices is None:
586	self.vertex_values = values
587	else:
588	for element_index, value in map(None, indices, values):
589	self.vertex_values[element_index] = value
590	else:
591	msg = 'Values array must be 1d or 2d'
592	raise msg
593
594
595	def set_values_from_quantity(self, q,
596	location, indices, verbose):
597	"""Set quantity values from specified quantity instance q
598
599	Location is ignored - vertices will always be used here.
600	"""
601
602
603	A = q.vertex_values
604
605	from Numeric import allclose
606	msg = 'Quantities are defined on different meshes. '+\
607	'This might be a case for implementing interpolation '+\
608	'between different meshes.'
609	assert allclose(A.shape, self.vertex_values.shape), msg
610
611	self.set_values(A, location='vertices',
612	indices=indices,
613	verbose=verbose)
614
615
616	def set_values_from_function(self, f,
617	location='vertices',
618	indices=None,
619	verbose=False):
620	"""Set values for quantity using specified function
621
622	Input
623
624	f: x, y -> z Function where x, y and z are arrays
625	location: Where values are to be stored.
626	Permissible options are: vertices, centroid, edges,
627	unique vertices
628	Default is "vertices"
629	indices:
630
631
632	"""
633
634	#FIXME: Should check that function returns something sensible and
635	#raise a meaningfull exception if it returns None for example
636
637	#FIXME: Should supply absolute coordinates
638
639
640	# Compute the function values and call set_values again
641	if location == 'centroids':
642	if indices is None:
643	indices = range(len(self))
644
645	V = take(self.domain.get_centroid_coordinates(), indices)
646	self.set_values(f(V[:,0], V[:,1]),
647	location=location,
648	indices=indices)
649
650	elif location == 'vertices':
651
652	M = self.domain.number_of_triangles
653	V = self.domain.get_vertex_coordinates()
654
655	x = V[:,0]; y = V[:,1];
656	values = f(x, y)
657
658
659	# FIXME (Ole): This code should replace all the
660	# rest of this function and it would work, except
661	# one unit test in test_region fails.
662	# If that could be resolved this one will be
663	# more robust and simple.
664
665	#values = reshape(values, (M,3))
666	#self.set_values(values,
667	# location='vertices',
668	# indices=indices)
669
670
671	# This should be removed
672	if is_scalar(values):
673	# Function returned a constant value
674	self.set_values_from_constant(values,
675	location, indices, verbose)
676	return
677
678	# This should be removed
679	if indices is None:
680	for j in range(3):
681	self.vertex_values[:,j] = values[j::3]
682	else:
683	#Brute force
684	for i in indices:
685	for j in range(3):
686	self.vertex_values[i,j] = values[3*i+j]
687
688
689	else:
690	raise 'Not implemented: %s' %location
691
692
693
694	def set_values_from_geospatial_data(self, geospatial_data, alpha,
695	location, indices,
696	verbose = False,
697	use_cache = False):
698	#FIXME: Use this function for the time being. Later move code in here
699
700	points = geospatial_data.get_data_points(absolute = False)
701	values = geospatial_data.get_attributes()
702	data_georef = geospatial_data.get_geo_reference()
703
704
705
706	from anuga.coordinate_transforms.geo_reference import Geo_reference
707
708
709	points = ensure_numeric(points, Float)
710	values = ensure_numeric(values, Float)
711
712	if location != 'vertices':
713	msg = 'set_values_from_points is only defined for '+\
714	'location=\'vertices\''
715	raise ms
716
717	coordinates = self.domain.get_nodes()
718	triangles = self.domain.triangles #FIXME
719
720
721	#Take care of georeferencing
722	if data_georef is None:
723	data_georef = Geo_reference()
724
725
726	mesh_georef = self.domain.geo_reference
727
728
729	# Call fit_interpolate.fit function
730	#args = (coordinates, triangles, points, values)
731	args = (points, )
732	kwargs = {'vertex_coordinates': coordinates,
733	'triangles': triangles,
734	'mesh': None,
735	'point_attributes': values,
736	'data_origin': data_georef.get_origin(),
737	'mesh_origin': mesh_georef.get_origin(),
738	'alpha': alpha,
739	'verbose': verbose}
740
741	vertex_attributes = apply(fit_to_mesh,
742	args, kwargs)
743
744	# Call underlying method using array values
745	self.set_values_from_array(vertex_attributes,
746	location, indices, verbose)
747
748
749
750	def set_values_from_points(self, points, values, alpha,
751	location, indices,
752	data_georef = None,
753	verbose = False,
754	use_cache = False):
755	"""
756	Set quantity values from arbitray data points using
757	fit_interpolate.fit
758	"""
759
760	raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
761
762
763	def set_values_from_file(self, filename, attribute_name, alpha,
764	location, indices,
765	verbose = False,
766	use_cache = False,
767	max_read_lines=None):
768	"""Set quantity based on arbitrary points in a points file
769	using attribute_name selects name of attribute
770	present in file.
771	If attribute_name is not specified, use first available attribute
772	as defined in geospatial_data.
773	"""
774
775	from types import StringType
776	msg = 'Filename must be a text string'
777	assert type(filename) == StringType, msg
778
779
780	if location != 'vertices':
781	msg = 'set_values_from_file is only defined for '+\
782	'location=\'vertices\''
783	raise msg
784
785	if True:
786	vertex_attributes = fit_to_mesh(filename,
787	mesh = self.domain,
788	alpha=alpha,
789	attribute_name=attribute_name,
790	use_cache=use_cache,
791	verbose=verbose,
792	max_read_lines=max_read_lines)
793	else:
794
795	coordinates = self.domain.get_nodes(absolute=True)
796	triangles = self.domain.triangles #FIXME
797	vertex_attributes = fit_to_mesh(filename,
798	coordinates, triangles,
799	alpha=alpha,
800	attribute_name=attribute_name,
801	use_cache=use_cache,
802	verbose=verbose,
803	max_read_lines=max_read_lines)
804
805	# Call underlying method using array values
806	self.set_values_from_array(vertex_attributes,
807	location, indices, verbose)
808
809
810	def get_extremum_index(self, mode=None, indices=None):
811	"""Return index for maximum or minimum value of quantity (on centroids)
812
813	Optional arguments:
814	mode is either 'max'(default) or 'min'.
815	indices is the set of element ids that the operation applies to.
816
817	Usage:
818	i = get_extreme_index()
819
820	Notes:
821	We do not seek the extremum at vertices as each vertex can
822	have multiple values - one for each triangle sharing it.
823
824	If there are multiple cells with same maximum value, the
825	first cell encountered in the triangle array is returned.
826	"""
827
828	V = self.get_values(location='centroids', indices=indices)
829
830	# Always return absolute indices
831	if mode is None or mode == 'max':
832	i = argmax(V)
833	elif mode == 'min':
834	i = argmin(V)
835
836
837	if indices is None:
838	return i
839	else:
840	return indices[i]
841
842
843	def get_maximum_index(self, indices=None):
844	"""See get extreme index for details
845	"""
846
847	return self.get_extremum_index(mode='max',
848	indices=indices)
849
850
851
852	def get_maximum_value(self, indices=None):
853	"""Return maximum value of quantity (on centroids)
854
855	Optional argument:
856	indices is the set of element ids that the operation applies to.
857
858	Usage:
859	v = get_maximum_value()
860
861	Note, we do not seek the maximum at vertices as each vertex can
862	have multiple values - one for each triangle sharing it
863	"""
864
865
866	i = self.get_maximum_index(indices)
867	V = self.get_values(location='centroids') #, indices=indices)
868
869	return V[i]
870
871
872	def get_maximum_location(self, indices=None):
873	"""Return location of maximum value of quantity (on centroids)
874
875	Optional argument:
876	indices is the set of element ids that the operation applies to.
877
878	Usage:
879	x, y = get_maximum_location()
880
881
882	Notes:
883	We do not seek the maximum at vertices as each vertex can
884	have multiple values - one for each triangle sharing it.
885
886	If there are multiple cells with same maximum value, the
887	first cell encountered in the triangle array is returned.
888	"""
889
890	i = self.get_maximum_index(indices)
891	x, y = self.domain.get_centroid_coordinates()[i]
892
893	return x, y
894
895
896	def get_minimum_index(self, indices=None):
897	"""See get extreme index for details
898	"""
899
900	return self.get_extremum_index(mode='min',
901	indices=indices)
902
903
904	def get_minimum_value(self, indices=None):
905	"""Return minimum value of quantity (on centroids)
906
907	Optional argument:
908	indices is the set of element ids that the operation applies to.
909
910	Usage:
911	v = get_minimum_value()
912
913	See get_maximum_value for more details.
914	"""
915
916
917	i = self.get_minimum_index(indices)
918	V = self.get_values(location='centroids')
919
920	return V[i]
921
922
923	def get_minimum_location(self, indices=None):
924	"""Return location of minimum value of quantity (on centroids)
925
926	Optional argument:
927	indices is the set of element ids that the operation applies to.
928
929	Usage:
930	x, y = get_minimum_location()
931
932
933	Notes:
934	We do not seek the maximum at vertices as each vertex can
935	have multiple values - one for each triangle sharing it.
936
937	If there are multiple cells with same maximum value, the
938	first cell encountered in the triangle array is returned.
939	"""
940
941	i = self.get_minimum_index(indices)
942	x, y = self.domain.get_centroid_coordinates()[i]
943
944	return x, y
945
946
947
948
949	def get_interpolated_values(self, interpolation_points):
950
951	# Interpolation object based on internal (discontinuous triangles)
952	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
953	smooth=False)
954	# FIXME: This concat should roll into get_vertex_values
955	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
956	axis=1)
957
958	can_reuse = False
959	if hasattr(self, 'interpolation_object'):
960	# Reuse to save time
961	I = self.interpolation_object
962
963	if allclose(interpolation_points, I._point_coordinates):
964	can_reuse = True
965
966
967	if can_reuse is True:
968	# Use absence of points to indicate reuse in I.interpolate
969	result = I.interpolate(vertex_values)
970	else:
971	from anuga.fit_interpolate.interpolate import Interpolate
972
973	# Create interpolation object with matrix
974	I = Interpolate(vertex_coordinates, triangles)
975	self.interpolation_object = I
976
977	# Call interpolate with points the first time
978	interpolation_points = ensure_numeric(interpolation_points, Float)
979	result = I.interpolate(vertex_values, interpolation_points)
980
981	return result
982
983
984	def get_values(self, interpolation_points=None,
985	location='vertices',
986	indices = None):
987	"""get values for quantity
988
989	return X, Compatible list, Numeric array (see below)
990	interpolation_points: List of x, y coordinates where value is
991	sought (using interpolation). If points are given, values of
992	location and indices are ignored
993
994	location: Where values are to be stored.
995	Permissible options are: vertices, edges, centroids
996	and unique vertices. Default is 'vertices'
997
998
999	The returned values with be a list the length of indices
1000	(N if indices = None).
1001
1002	In case of location == 'centroids' the dimension of returned
1003	values will be a list or a Numerical array of length N, N being
1004	the number of elements.
1005
1006	In case of location == 'vertices' or 'edges' the dimension of
1007	returned values will be of dimension Nx3
1008
1009	In case of location == 'unique vertices' the average value at
1010	each vertex will be returned and the dimension of returned values
1011	will be a 1d array of length "number of vertices"
1012
1013	Indices is the set of element ids that the operation applies to.
1014
1015	The values will be stored in elements following their
1016	internal ordering.
1017	"""
1018	from Numeric import take
1019
1020	# FIXME (Ole): I reckon we should have the option of passing a
1021	# polygon into get_values. The question becomes how
1022	# resulting values should be ordered.
1023
1024	if interpolation_points is not None:
1025	return self.get_interpolated_values(interpolation_points)
1026
1027
1028
1029	if location not in ['vertices', 'centroids', 'edges',
1030	'unique vertices']:
1031	msg = 'Invalid location: %s' %location
1032	raise msg
1033
1034	import types, Numeric
1035	assert type(indices) in [types.ListType, types.NoneType,
1036	Numeric.ArrayType],\
1037	'Indices must be a list or None'
1038
1039	if location == 'centroids':
1040	if (indices == None):
1041	indices = range(len(self))
1042	return take(self.centroid_values,indices)
1043	elif location == 'edges':
1044	if (indices == None):
1045	indices = range(len(self))
1046	return take(self.edge_values,indices)
1047	elif location == 'unique vertices':
1048	if (indices == None):
1049	indices=range(self.domain.number_of_nodes)
1050	vert_values = []
1051	#Go through list of unique vertices
1052	for unique_vert_id in indices:
1053	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1054
1055	#In case there are unused points
1056	if len(triangles) == 0:
1057	msg = 'Unique vertex not associated with triangles'
1058	raise msg
1059
1060	# Go through all triangle, vertex pairs
1061	# Average the values
1062
1063	# FIXME (Ole): Should we merge this with get_vertex_values
1064	sum = 0
1065	for triangle_id, vertex_id in triangles:
1066	sum += self.vertex_values[triangle_id, vertex_id]
1067	vert_values.append(sum/len(triangles))
1068	return Numeric.array(vert_values)
1069	else:
1070	if (indices is None):
1071	indices = range(len(self))
1072	return take(self.vertex_values, indices)
1073
1074
1075
1076	def set_vertex_values(self, A, indices = None):
1077	"""Set vertex values for all unique vertices based on input array A
1078	which has one entry per unique vertex, i.e.
1079	one value for each row in array self.domain.nodes.
1080
1081	indices is the list of vertex_id's that will be set.
1082
1083	This function is used by set_values_from_array
1084	"""
1085
1086	from Numeric import array, Float
1087
1088	#Assert that A can be converted to a Numeric array of appropriate dim
1089	A = array(A, Float)
1090
1091	#print 'SHAPE A', A.shape
1092	assert len(A.shape) == 1
1093
1094	if indices is None:
1095	assert A.shape[0] == self.domain.get_nodes().shape[0]
1096	vertex_list = range(A.shape[0])
1097	else:
1098	assert A.shape[0] == len(indices)
1099	vertex_list = indices
1100
1101	#Go through list of unique vertices
1102
1103	for i_index, unique_vert_id in enumerate(vertex_list):
1104
1105
1106	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1107
1108	#In case there are unused points
1109	if len(triangles) == 0: continue
1110
1111	#Go through all triangle, vertex pairs
1112	#touching vertex unique_vert_id and set corresponding vertex value
1113	for triangle_id, vertex_id in triangles:
1114	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1115
1116	#Intialise centroid and edge_values
1117	self.interpolate()
1118
1119
1120	def smooth_vertex_values(self, value_array='field_values',
1121	precision = None):
1122	""" Smooths field_values or conserved_quantities data.
1123	TODO: be able to smooth individual fields
1124	NOTE: This function does not have a test.
1125	FIXME: NOT DONE - do we need it?
1126	FIXME: this function isn't called by anything.
1127	Maybe it should be removed..-DSG
1128	"""
1129
1130	from Numeric import concatenate, zeros, Float, Int, array, reshape
1131
1132
1133	A,V = self.get_vertex_values(xy=False,
1134	value_array=value_array,
1135	smooth = True,
1136	precision = precision)
1137
1138	#Set some field values
1139	for volume in self:
1140	for i,v in enumerate(volume.vertices):
1141	if value_array == 'field_values':
1142	volume.set_field_values('vertex', i, A[v,:])
1143	elif value_array == 'conserved_quantities':
1144	volume.set_conserved_quantities('vertex', i, A[v,:])
1145
1146	if value_array == 'field_values':
1147	self.precompute()
1148	elif value_array == 'conserved_quantities':
1149	Volume.interpolate_conserved_quantities()
1150
1151
1152	# Methods for outputting model results
1153	def get_vertex_values(self,
1154	xy=True,
1155	smooth=None,
1156	precision=None):
1157	"""Return vertex values like an OBJ format i.e. one value per node.
1158
1159	The vertex values are returned as one sequence in the 1D float array A.
1160	If requested the coordinates will be returned in 1D arrays X and Y.
1161
1162	The connectivity is represented as an integer array, V, of dimension
1163	Mx3, where M is the number of triangles. Each row has three indices
1164	defining the triangle and they correspond to elements in the arrays
1165	X, Y and A.
1166
1167	if smooth is True, vertex values corresponding to one common
1168	coordinate set will be smoothed by taking the average of vertex values for each node.
1169	In this case vertex coordinates will be
1170	de-duplicated corresponding to the original nodes as obtained from
1171	the method general_mesh.get_nodes()
1172
1173	If no smoothings is required, vertex coordinates and values will
1174	be aggregated as a concatenation of values at
1175	vertices 0, vertices 1 and vertices 2. This corresponds to
1176	the node coordinates obtained from the method
1177	general_mesh.get_vertex_coordinates()
1178
1179
1180	Calling convention
1181	if xy is True:
1182	X,Y,A,V = get_vertex_values
1183	else:
1184	A,V = get_vertex_values
1185
1186	"""
1187
1188	from Numeric import concatenate, zeros, Float, Int, array, reshape
1189
1190
1191	if smooth is None:
1192	# Take default from domain
1193	smooth = self.domain.smooth
1194
1195	if precision is None:
1196	precision = Float
1197
1198
1199	if smooth is True:
1200	# Ensure continuous vertex values by averaging
1201	# values at each node
1202
1203	V = self.domain.get_triangles()
1204	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1205	A = zeros(N, Float)
1206	points = self.domain.get_nodes()
1207
1208	if 1:
1209	# Fast C version
1210	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1211	ensure_numeric(self.domain.number_of_triangles_per_node),
1212	ensure_numeric(self.vertex_values),
1213	A)
1214	A = A.astype(precision)
1215	else:
1216
1217	# Slow Python version
1218
1219	current_node = 0
1220	k = 0 # Track triangles touching on node
1221	total = 0.0
1222	for index in self.domain.vertex_value_indices:
1223	if current_node == N:
1224	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1225	raise msg
1226
1227
1228
1229	k += 1
1230
1231	volume_id = index / 3
1232	vertex_id = index % 3
1233
1234	#assert V[volume_id, vertex_id] == current_node
1235
1236	v = self.vertex_values[volume_id, vertex_id]
1237	total += v
1238
1239	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1240	if self.domain.number_of_triangles_per_node[current_node] == k:
1241	A[current_node] = total/k
1242
1243
1244	# Move on to next node
1245	total = 0.0
1246	k = 0
1247	current_node += 1
1248
1249
1250
1251	else:
1252	# Allow discontinuous vertex values
1253	V = self.domain.get_disconnected_triangles()
1254	points = self.domain.get_vertex_coordinates()
1255	A = self.vertex_values.flat.astype(precision)
1256
1257
1258	# Return
1259	if xy is True:
1260	X = points[:,0].astype(precision)
1261	Y = points[:,1].astype(precision)
1262
1263	return X, Y, A, V
1264	else:
1265	return A, V
1266
1267
1268
1269	def extrapolate_first_order(self):
1270	"""Extrapolate conserved quantities from centroid to
1271	vertices for each volume using
1272	first order scheme.
1273	"""
1274
1275	qc = self.centroid_values
1276	qv = self.vertex_values
1277
1278	for i in range(3):
1279	qv[:,i] = qc
1280
1281
1282	def get_integral(self):
1283	"""Compute the integral of quantity across entire domain
1284	"""
1285	integral = 0
1286	for k in range(len(self.domain)):
1287	area = self.domain.areas[k]
1288	qc = self.centroid_values[k]
1289	integral += qc*area
1290
1291	return integral
1292
1293
1294
1295
1296	class Conserved_quantity(Quantity):
1297	"""Class conserved quantity adds to Quantity:
1298
1299	boundary values, storage and method for updating, and
1300	methods for (second order) extrapolation from centroid to vertices inluding
1301	gradients and limiters
1302	"""
1303
1304	def __init__(self, domain, vertex_values=None):
1305	Quantity.__init__(self, domain, vertex_values)
1306
1307	from Numeric import zeros, Float
1308
1309	#Allocate space for boundary values
1310	L = len(domain.boundary)
1311	self.boundary_values = zeros(L, Float)
1312
1313	#Allocate space for updates of conserved quantities by
1314	#flux calculations and forcing functions
1315
1316	N = len(domain) # number_of_triangles
1317	self.explicit_update = zeros(N, Float )
1318	self.semi_implicit_update = zeros(N, Float )
1319	self.centroid_backup_values = zeros(N, Float)
1320
1321
1322
1323
1324	def update(self, timestep):
1325	#Call correct module function
1326	#(either from this module or C-extension)
1327	return update(self, timestep)
1328
1329
1330	def compute_gradients(self):
1331	#Call correct module function
1332	#(either from this module or C-extension)
1333	return compute_gradients(self)
1334
1335	def limit(self):
1336	#Call correct module function
1337	#(either from this module or C-extension)
1338	limit(self)
1339
1340	def extrapolate_second_order(self):
1341	#Call correct module function
1342	#(either from this module or C-extension)
1343	extrapolate_second_order(self)
1344
1345	def backup_centroid_values(self):
1346	#Call correct module function
1347	#(either from this module or C-extension)
1348	backup_centroid_values(self)
1349
1350	def saxpy_centroid_values(self,a,b):
1351	#Call correct module function
1352	#(either from this module or C-extension)
1353	saxpy_centroid_values(self,a,b)
1354
1355
1356
1357	from anuga.utilities import compile
1358	if compile.can_use_C_extension('quantity_ext.c'):
1359	# Underlying C implementations can be accessed
1360
1361	from quantity_ext import average_vertex_values, backup_centroid_values, \
1362	saxpy_centroid_values
1363
1364	from quantity_ext import compute_gradients, limit,\
1365	extrapolate_second_order, interpolate_from_vertices_to_edges, update
1366	else:
1367	msg = 'C implementations could not be accessed by %s.\n ' %__file__
1368	msg += 'Make sure compile_all.py has been run as described in '
1369	msg += 'the ANUGA installation guide.'
1370	raise Exception, msg
1371
1372
1373

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

Download in other formats: