Context Navigation

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

Last change on this file since 4779 was 4779, checked in by duncan, 16 years ago
reduce memory use in quantity.set_value. fit_to_mesh can now use an existing mesh instance, which it does in quantity.set_value.
File size: 46.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	args = (points, )
730	kwargs = {'vertex_coordinates': coordinates,
731	'triangles': triangles,
732	'mesh': None,
733	'point_attributes': values,
734	'data_origin': data_georef.get_origin(),
735	'mesh_origin': mesh_georef.get_origin(),
736	'alpha': alpha,
737	'verbose': verbose}
738
739	vertex_attributes = apply(fit_to_mesh,
740	args, kwargs)
741
742	# Call underlying method using array values
743	self.set_values_from_array(vertex_attributes,
744	location, indices, verbose)
745
746
747
748	def set_values_from_points(self, points, values, alpha,
749	location, indices,
750	data_georef = None,
751	verbose = False,
752	use_cache = False):
753	"""
754	Set quantity values from arbitray data points using
755	fit_interpolate.fit
756	"""
757
758	raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
759
760
761	def set_values_from_file(self, filename, attribute_name, alpha,
762	location, indices,
763	verbose = False,
764	use_cache = False,
765	max_read_lines=None):
766	"""Set quantity based on arbitrary points in a points file
767	using attribute_name selects name of attribute
768	present in file.
769	If attribute_name is not specified, use first available attribute
770	as defined in geospatial_data.
771	"""
772
773	from types import StringType
774	msg = 'Filename must be a text string'
775	assert type(filename) == StringType, msg
776
777
778	if location != 'vertices':
779	msg = 'set_values_from_file is only defined for '+\
780	'location=\'vertices\''
781	raise msg
782
783	#coordinates = self.domain.get_nodes(absolute=True)
784	#triangles = self.domain.triangles #FIXME
785	#vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
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
794	# Call underlying method using array values
795	self.set_values_from_array(vertex_attributes,
796	location, indices, verbose)
797
798
799	def get_extremum_index(self, mode=None, indices=None):
800	"""Return index for maximum or minimum value of quantity (on centroids)
801
802	Optional arguments:
803	mode is either 'max'(default) or 'min'.
804	indices is the set of element ids that the operation applies to.
805
806	Usage:
807	i = get_extreme_index()
808
809	Notes:
810	We do not seek the extremum at vertices as each vertex can
811	have multiple values - one for each triangle sharing it.
812
813	If there are multiple cells with same maximum value, the
814	first cell encountered in the triangle array is returned.
815	"""
816
817	V = self.get_values(location='centroids', indices=indices)
818
819	# Always return absolute indices
820	if mode is None or mode == 'max':
821	i = argmax(V)
822	elif mode == 'min':
823	i = argmin(V)
824
825
826	if indices is None:
827	return i
828	else:
829	return indices[i]
830
831
832	def get_maximum_index(self, indices=None):
833	"""See get extreme index for details
834	"""
835
836	return self.get_extremum_index(mode='max',
837	indices=indices)
838
839
840
841	def get_maximum_value(self, indices=None):
842	"""Return maximum value of quantity (on centroids)
843
844	Optional argument:
845	indices is the set of element ids that the operation applies to.
846
847	Usage:
848	v = get_maximum_value()
849
850	Note, we do not seek the maximum at vertices as each vertex can
851	have multiple values - one for each triangle sharing it
852	"""
853
854
855	i = self.get_maximum_index(indices)
856	V = self.get_values(location='centroids') #, indices=indices)
857
858	return V[i]
859
860
861	def get_maximum_location(self, indices=None):
862	"""Return location of maximum value of quantity (on centroids)
863
864	Optional argument:
865	indices is the set of element ids that the operation applies to.
866
867	Usage:
868	x, y = get_maximum_location()
869
870
871	Notes:
872	We do not seek the maximum at vertices as each vertex can
873	have multiple values - one for each triangle sharing it.
874
875	If there are multiple cells with same maximum value, the
876	first cell encountered in the triangle array is returned.
877	"""
878
879	i = self.get_maximum_index(indices)
880	x, y = self.domain.get_centroid_coordinates()[i]
881
882	return x, y
883
884
885	def get_minimum_index(self, indices=None):
886	"""See get extreme index for details
887	"""
888
889	return self.get_extremum_index(mode='min',
890	indices=indices)
891
892
893	def get_minimum_value(self, indices=None):
894	"""Return minimum value of quantity (on centroids)
895
896	Optional argument:
897	indices is the set of element ids that the operation applies to.
898
899	Usage:
900	v = get_minimum_value()
901
902	See get_maximum_value for more details.
903	"""
904
905
906	i = self.get_minimum_index(indices)
907	V = self.get_values(location='centroids')
908
909	return V[i]
910
911
912	def get_minimum_location(self, indices=None):
913	"""Return location of minimum value of quantity (on centroids)
914
915	Optional argument:
916	indices is the set of element ids that the operation applies to.
917
918	Usage:
919	x, y = get_minimum_location()
920
921
922	Notes:
923	We do not seek the maximum at vertices as each vertex can
924	have multiple values - one for each triangle sharing it.
925
926	If there are multiple cells with same maximum value, the
927	first cell encountered in the triangle array is returned.
928	"""
929
930	i = self.get_minimum_index(indices)
931	x, y = self.domain.get_centroid_coordinates()[i]
932
933	return x, y
934
935
936
937
938	def get_interpolated_values(self, interpolation_points):
939
940	# Interpolation object based on internal (discontinuous triangles)
941	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
942	smooth=False)
943	# FIXME: This concat should roll into get_vertex_values
944	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
945	axis=1)
946
947	can_reuse = False
948	if hasattr(self, 'interpolation_object'):
949	# Reuse to save time
950	I = self.interpolation_object
951
952	if allclose(interpolation_points, I._point_coordinates):
953	can_reuse = True
954
955
956	if can_reuse is True:
957	# Use absence of points to indicate reuse in I.interpolate
958	result = I.interpolate(vertex_values)
959	else:
960	from anuga.fit_interpolate.interpolate import Interpolate
961
962	# Create interpolation object with matrix
963	I = Interpolate(vertex_coordinates, triangles)
964	self.interpolation_object = I
965
966	# Call interpolate with points the first time
967	interpolation_points = ensure_numeric(interpolation_points, Float)
968	result = I.interpolate(vertex_values, interpolation_points)
969
970	return result
971
972
973	def get_values(self, interpolation_points=None,
974	location='vertices',
975	indices = None):
976	"""get values for quantity
977
978	return X, Compatible list, Numeric array (see below)
979	interpolation_points: List of x, y coordinates where value is
980	sought (using interpolation). If points are given, values of
981	location and indices are ignored
982
983	location: Where values are to be stored.
984	Permissible options are: vertices, edges, centroids
985	and unique vertices. Default is 'vertices'
986
987
988	The returned values with be a list the length of indices
989	(N if indices = None).
990
991	In case of location == 'centroids' the dimension of returned
992	values will be a list or a Numerical array of length N, N being
993	the number of elements.
994
995	In case of location == 'vertices' or 'edges' the dimension of
996	returned values will be of dimension Nx3
997
998	In case of location == 'unique vertices' the average value at
999	each vertex will be returned and the dimension of returned values
1000	will be a 1d array of length "number of vertices"
1001
1002	Indices is the set of element ids that the operation applies to.
1003
1004	The values will be stored in elements following their
1005	internal ordering.
1006	"""
1007	from Numeric import take
1008
1009	# FIXME (Ole): I reckon we should have the option of passing a
1010	# polygon into get_values. The question becomes how
1011	# resulting values should be ordered.
1012
1013	if interpolation_points is not None:
1014	return self.get_interpolated_values(interpolation_points)
1015
1016
1017
1018	if location not in ['vertices', 'centroids', 'edges',
1019	'unique vertices']:
1020	msg = 'Invalid location: %s' %location
1021	raise msg
1022
1023	import types, Numeric
1024	assert type(indices) in [types.ListType, types.NoneType,
1025	Numeric.ArrayType],\
1026	'Indices must be a list or None'
1027
1028	if location == 'centroids':
1029	if (indices == None):
1030	indices = range(len(self))
1031	return take(self.centroid_values,indices)
1032	elif location == 'edges':
1033	if (indices == None):
1034	indices = range(len(self))
1035	return take(self.edge_values,indices)
1036	elif location == 'unique vertices':
1037	if (indices == None):
1038	indices=range(self.domain.number_of_nodes)
1039	vert_values = []
1040	#Go through list of unique vertices
1041	for unique_vert_id in indices:
1042	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1043
1044	#In case there are unused points
1045	if len(triangles) == 0:
1046	msg = 'Unique vertex not associated with triangles'
1047	raise msg
1048
1049	# Go through all triangle, vertex pairs
1050	# Average the values
1051
1052	# FIXME (Ole): Should we merge this with get_vertex_values
1053	sum = 0
1054	for triangle_id, vertex_id in triangles:
1055	sum += self.vertex_values[triangle_id, vertex_id]
1056	vert_values.append(sum/len(triangles))
1057	return Numeric.array(vert_values)
1058	else:
1059	if (indices is None):
1060	indices = range(len(self))
1061	return take(self.vertex_values, indices)
1062
1063
1064
1065	def set_vertex_values(self, A, indices = None):
1066	"""Set vertex values for all unique vertices based on input array A
1067	which has one entry per unique vertex, i.e.
1068	one value for each row in array self.domain.nodes.
1069
1070	indices is the list of vertex_id's that will be set.
1071
1072	This function is used by set_values_from_array
1073	"""
1074
1075	from Numeric import array, Float
1076
1077	#Assert that A can be converted to a Numeric array of appropriate dim
1078	A = array(A, Float)
1079
1080	#print 'SHAPE A', A.shape
1081	assert len(A.shape) == 1
1082
1083	if indices is None:
1084	assert A.shape[0] == self.domain.get_nodes().shape[0]
1085	vertex_list = range(A.shape[0])
1086	else:
1087	assert A.shape[0] == len(indices)
1088	vertex_list = indices
1089
1090	#Go through list of unique vertices
1091
1092	for i_index, unique_vert_id in enumerate(vertex_list):
1093
1094
1095	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1096
1097	#In case there are unused points
1098	if len(triangles) == 0: continue
1099
1100	#Go through all triangle, vertex pairs
1101	#touching vertex unique_vert_id and set corresponding vertex value
1102	for triangle_id, vertex_id in triangles:
1103	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1104
1105	#Intialise centroid and edge_values
1106	self.interpolate()
1107
1108
1109	def smooth_vertex_values(self, value_array='field_values',
1110	precision = None):
1111	""" Smooths field_values or conserved_quantities data.
1112	TODO: be able to smooth individual fields
1113	NOTE: This function does not have a test.
1114	FIXME: NOT DONE - do we need it?
1115	FIXME: this function isn't called by anything.
1116	Maybe it should be removed..-DSG
1117	"""
1118
1119	from Numeric import concatenate, zeros, Float, Int, array, reshape
1120
1121
1122	A,V = self.get_vertex_values(xy=False,
1123	value_array=value_array,
1124	smooth = True,
1125	precision = precision)
1126
1127	#Set some field values
1128	for volume in self:
1129	for i,v in enumerate(volume.vertices):
1130	if value_array == 'field_values':
1131	volume.set_field_values('vertex', i, A[v,:])
1132	elif value_array == 'conserved_quantities':
1133	volume.set_conserved_quantities('vertex', i, A[v,:])
1134
1135	if value_array == 'field_values':
1136	self.precompute()
1137	elif value_array == 'conserved_quantities':
1138	Volume.interpolate_conserved_quantities()
1139
1140
1141	# Methods for outputting model results
1142	def get_vertex_values(self,
1143	xy=True,
1144	smooth=None,
1145	precision=None):
1146	"""Return vertex values like an OBJ format i.e. one value per node.
1147
1148	The vertex values are returned as one sequence in the 1D float array A.
1149	If requested the coordinates will be returned in 1D arrays X and Y.
1150
1151	The connectivity is represented as an integer array, V, of dimension
1152	Mx3, where M is the number of triangles. Each row has three indices
1153	defining the triangle and they correspond to elements in the arrays
1154	X, Y and A.
1155
1156	if smooth is True, vertex values corresponding to one common
1157	coordinate set will be smoothed by taking the average of vertex values for each node.
1158	In this case vertex coordinates will be
1159	de-duplicated corresponding to the original nodes as obtained from
1160	the method general_mesh.get_nodes()
1161
1162	If no smoothings is required, vertex coordinates and values will
1163	be aggregated as a concatenation of values at
1164	vertices 0, vertices 1 and vertices 2. This corresponds to
1165	the node coordinates obtained from the method
1166	general_mesh.get_vertex_coordinates()
1167
1168
1169	Calling convention
1170	if xy is True:
1171	X,Y,A,V = get_vertex_values
1172	else:
1173	A,V = get_vertex_values
1174
1175	"""
1176
1177	from Numeric import concatenate, zeros, Float, Int, array, reshape
1178
1179
1180	if smooth is None:
1181	# Take default from domain
1182	smooth = self.domain.smooth
1183
1184	if precision is None:
1185	precision = Float
1186
1187
1188	if smooth is True:
1189	# Ensure continuous vertex values by averaging
1190	# values at each node
1191
1192	V = self.domain.get_triangles()
1193	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1194	A = zeros(N, Float)
1195	points = self.domain.get_nodes()
1196
1197	if 1:
1198	# Fast C version
1199	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1200	ensure_numeric(self.domain.number_of_triangles_per_node),
1201	ensure_numeric(self.vertex_values),
1202	A)
1203	A = A.astype(precision)
1204	else:
1205
1206	# Slow Python version
1207
1208	current_node = 0
1209	k = 0 # Track triangles touching on node
1210	total = 0.0
1211	for index in self.domain.vertex_value_indices:
1212	if current_node == N:
1213	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1214	raise msg
1215
1216
1217
1218	k += 1
1219
1220	volume_id = index / 3
1221	vertex_id = index % 3
1222
1223	#assert V[volume_id, vertex_id] == current_node
1224
1225	v = self.vertex_values[volume_id, vertex_id]
1226	total += v
1227
1228	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1229	if self.domain.number_of_triangles_per_node[current_node] == k:
1230	A[current_node] = total/k
1231
1232
1233	# Move on to next node
1234	total = 0.0
1235	k = 0
1236	current_node += 1
1237
1238
1239
1240	else:
1241	# Allow discontinuous vertex values
1242	V = self.domain.get_disconnected_triangles()
1243	points = self.domain.get_vertex_coordinates()
1244	A = self.vertex_values.flat.astype(precision)
1245
1246
1247	# Return
1248	if xy is True:
1249	X = points[:,0].astype(precision)
1250	Y = points[:,1].astype(precision)
1251
1252	return X, Y, A, V
1253	else:
1254	return A, V
1255
1256
1257
1258	def extrapolate_first_order(self):
1259	"""Extrapolate conserved quantities from centroid to
1260	vertices for each volume using
1261	first order scheme.
1262	"""
1263
1264	qc = self.centroid_values
1265	qv = self.vertex_values
1266
1267	for i in range(3):
1268	qv[:,i] = qc
1269
1270
1271	def get_integral(self):
1272	"""Compute the integral of quantity across entire domain
1273	"""
1274	integral = 0
1275	for k in range(len(self.domain)):
1276	area = self.domain.areas[k]
1277	qc = self.centroid_values[k]
1278	integral += qc*area
1279
1280	return integral
1281
1282
1283
1284
1285	class Conserved_quantity(Quantity):
1286	"""Class conserved quantity adds to Quantity:
1287
1288	boundary values, storage and method for updating, and
1289	methods for (second order) extrapolation from centroid to vertices inluding
1290	gradients and limiters
1291	"""
1292
1293	def __init__(self, domain, vertex_values=None):
1294	Quantity.__init__(self, domain, vertex_values)
1295
1296	from Numeric import zeros, Float
1297
1298	#Allocate space for boundary values
1299	L = len(domain.boundary)
1300	self.boundary_values = zeros(L, Float)
1301
1302	#Allocate space for updates of conserved quantities by
1303	#flux calculations and forcing functions
1304
1305	N = len(domain) # number_of_triangles
1306	self.explicit_update = zeros(N, Float )
1307	self.semi_implicit_update = zeros(N, Float )
1308	self.centroid_backup_values = zeros(N, Float)
1309
1310
1311
1312
1313	def update(self, timestep):
1314	#Call correct module function
1315	#(either from this module or C-extension)
1316	return update(self, timestep)
1317
1318
1319	def compute_gradients(self):
1320	#Call correct module function
1321	#(either from this module or C-extension)
1322	return compute_gradients(self)
1323
1324	def limit(self):
1325	#Call correct module function
1326	#(either from this module or C-extension)
1327	limit(self)
1328
1329	def extrapolate_second_order(self):
1330	#Call correct module function
1331	#(either from this module or C-extension)
1332	extrapolate_second_order(self)
1333
1334	def backup_centroid_values(self):
1335	#Call correct module function
1336	#(either from this module or C-extension)
1337	backup_centroid_values(self)
1338
1339	def saxpy_centroid_values(self,a,b):
1340	#Call correct module function
1341	#(either from this module or C-extension)
1342	saxpy_centroid_values(self,a,b)
1343
1344
1345
1346	from anuga.utilities import compile
1347	if compile.can_use_C_extension('quantity_ext.c'):
1348	# Underlying C implementations can be accessed
1349
1350	from quantity_ext import average_vertex_values, backup_centroid_values, \
1351	saxpy_centroid_values
1352
1353	from quantity_ext import compute_gradients, limit,\
1354	extrapolate_second_order, interpolate_from_vertices_to_edges, update
1355	else:
1356	msg = 'C implementations could not be accessed by %s.\n ' %__file__
1357	msg += 'Make sure compile_all.py has been run as described in '
1358	msg += 'the ANUGA installation guide.'
1359	raise Exception, msg
1360
1361
1362

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

Download in other formats: