Context Navigation

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

Last change on this file since 4735 was 4735, checked in by ole, 16 years ago
Cleaned up some obsolete code
File size: 52.1 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	assert values.shape[0] == indices.shape[0], msg
531
532	N = self.centroid_values.shape[0]
533
534	if location == 'centroids':
535	assert len(values.shape) == 1, 'Values array must be 1d'
536
537	if indices is None:
538	msg = 'Number of values must match number of elements'
539	assert values.shape[0] == N, msg
540
541	self.centroid_values = values
542	else:
543	msg = 'Number of values must match number of indices'
544	assert values.shape[0] == indices.shape[0], msg
545
546	#Brute force
547	for i in range(len(indices)):
548	self.centroid_values[indices[i]] = values[i]
549
550	elif location == 'edges':
551	# FIXME (Ole): No mention of indices here. However, I don't
552	# think we ever need to set values at edges anyway
553	assert len(values.shape) == 2, 'Values array must be 2d'
554
555	msg = 'Number of values must match number of elements'
556	assert values.shape[0] == N, msg
557
558	msg = 'Array must be N x 3'
559	assert values.shape[1] == 3, msg
560
561	self.edge_values = values
562
563	elif location == 'unique vertices':
564	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
565	'Values array must be 1d'
566
567	self.set_vertex_values(values.flat, indices=indices)
568
569	else:
570	# Location vertices
571	if len(values.shape) == 1:
572	self.set_vertex_values(values, indices=indices)
573
574	elif len(values.shape) == 2:
575	#Vertex values are given as a triplet for each triangle
576
577	msg = 'Array must be N x 3'
578	assert values.shape[1] == 3, msg
579
580	if indices is None:
581	self.vertex_values = values
582	else:
583	for element_index, value in map(None, indices, values):
584	self.vertex_values[element_index] = value
585	else:
586	msg = 'Values array must be 1d or 2d'
587	raise msg
588
589
590	def set_values_from_quantity(self, q,
591	location, indices, verbose):
592	"""Set quantity values from specified quantity instance q
593
594	Location is ignored - vertices will always be used here.
595	"""
596
597
598	A = q.vertex_values
599
600	from Numeric import allclose
601	msg = 'Quantities are defined on different meshes. '+\
602	'This might be a case for implementing interpolation '+\
603	'between different meshes.'
604	assert allclose(A.shape, self.vertex_values.shape), msg
605
606	self.set_values(A, location='vertices',
607	indices=indices,
608	verbose=verbose)
609
610
611	def set_values_from_function(self, f,
612	location='vertices',
613	indices=None,
614	verbose=False):
615	"""Set values for quantity using specified function
616
617	Input
618
619	f: x, y -> z Function where x, y and z are arrays
620	location: Where values are to be stored.
621	Permissible options are: vertices, centroid, edges,
622	unique vertices
623	Default is "vertices"
624	indices:
625
626
627	"""
628
629	#FIXME: Should check that function returns something sensible and
630	#raise a meaningfull exception if it returns None for example
631
632	#FIXME: Should supply absolute coordinates
633
634
635	# Compute the function values and call set_values again
636	if location == 'centroids':
637	if indices is None:
638	indices = range(len(self))
639
640	V = take(self.domain.get_centroid_coordinates(), indices)
641	self.set_values(f(V[:,0], V[:,1]),
642	location=location,
643	indices=indices)
644
645	elif location == 'vertices':
646
647	M = self.domain.number_of_triangles
648	V = self.domain.get_vertex_coordinates()
649
650	x = V[:,0]; y = V[:,1];
651	values = f(x, y)
652
653
654	# FIXME (Ole): This code should replace all the
655	# rest of this function and it would work, except
656	# one unit test in test_region fails.
657	# If that could be resolved this one will be
658	# more robust and simple.
659
660	#values = reshape(values, (M,3))
661	#self.set_values(values,
662	# location='vertices',
663	# indices=indices)
664
665
666	# This should be removed
667	if is_scalar(values):
668	# Function returned a constant value
669	self.set_values_from_constant(values,
670	location, indices, verbose)
671	return
672
673	# This should be removed
674	if indices is None:
675	for j in range(3):
676	self.vertex_values[:,j] = values[j::3]
677	else:
678	#Brute force
679	for i in indices:
680	for j in range(3):
681	self.vertex_values[i,j] = values[3*i+j]
682
683
684	else:
685	raise 'Not implemented: %s' %location
686
687
688
689	def set_values_from_geospatial_data(self, geospatial_data, alpha,
690	location, indices,
691	verbose = False,
692	use_cache = False):
693
694	#FIXME: Use this function for the time being. Later move code in here
695
696	points = geospatial_data.get_data_points(absolute = False)
697	values = geospatial_data.get_attributes()
698	data_georef = geospatial_data.get_geo_reference()
699
700
701
702	from anuga.coordinate_transforms.geo_reference import Geo_reference
703
704
705	points = ensure_numeric(points, Float)
706	values = ensure_numeric(values, Float)
707
708	if location != 'vertices':
709	msg = 'set_values_from_points is only defined for '+\
710	'location=\'vertices\''
711	raise ms
712
713	coordinates = self.domain.get_nodes()
714	triangles = self.domain.triangles #FIXME
715
716
717	#Take care of georeferencing
718	if data_georef is None:
719	data_georef = Geo_reference()
720
721
722	mesh_georef = self.domain.geo_reference
723
724
725	# Call fit_interpolate.fit function
726	args = (coordinates, triangles, points, values)
727	kwargs = {'data_origin': data_georef.get_origin(),
728	'mesh_origin': mesh_georef.get_origin(),
729	'alpha': alpha,
730	'verbose': verbose}
731
732	vertex_attributes = apply(fit_to_mesh,
733	args, kwargs)
734
735	# Call underlying method using array values
736	self.set_values_from_array(vertex_attributes,
737	location, indices, verbose)
738
739
740
741	def set_values_from_points(self, points, values, alpha,
742	location, indices,
743	data_georef = None,
744	verbose = False,
745	use_cache = False):
746	"""
747	Set quantity values from arbitray data points using
748	fit_interpolate.fit
749	"""
750
751	raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
752
753
754	def set_values_from_file(self, filename, attribute_name, alpha,
755	location, indices,
756	verbose = False,
757	use_cache = False,
758	max_read_lines=None):
759	"""Set quantity based on arbitrary points in a points file
760	using attribute_name selects name of attribute
761	present in file.
762	If attribute_name is not specified, use first available attribute
763	as defined in geospatial_data.
764	"""
765
766	from types import StringType
767	msg = 'Filename must be a text string'
768	assert type(filename) == StringType, msg
769
770
771	if location != 'vertices':
772	msg = 'set_values_from_file is only defined for '+\
773	'location=\'vertices\''
774	raise msg
775
776	coordinates = self.domain.get_nodes(absolute=True)
777	triangles = self.domain.triangles #FIXME
778
779	vertex_attributes = fit_to_mesh(coordinates, triangles, filename,
780	alpha=alpha,
781	attribute_name=attribute_name,
782	use_cache=use_cache,
783	verbose=verbose,
784	max_read_lines=max_read_lines)
785
786	# Call underlying method using array values
787	self.set_values_from_array(vertex_attributes,
788	location, indices, verbose)
789
790
791	def get_extremum_index(self, mode=None, indices=None):
792	"""Return index for maximum or minimum value of quantity (on centroids)
793
794	Optional arguments:
795	mode is either 'max'(default) or 'min'.
796	indices is the set of element ids that the operation applies to.
797
798	Usage:
799	i = get_extreme_index()
800
801	Notes:
802	We do not seek the extremum at vertices as each vertex can
803	have multiple values - one for each triangle sharing it.
804
805	If there are multiple cells with same maximum value, the
806	first cell encountered in the triangle array is returned.
807	"""
808
809	V = self.get_values(location='centroids', indices=indices)
810
811	# Always return absolute indices
812	if mode is None or mode == 'max':
813	i = argmax(V)
814	elif mode == 'min':
815	i = argmin(V)
816
817
818	if indices is None:
819	return i
820	else:
821	return indices[i]
822
823
824	def get_maximum_index(self, indices=None):
825	"""See get extreme index for details
826	"""
827
828	return self.get_extremum_index(mode='max',
829	indices=indices)
830
831
832
833	def get_maximum_value(self, indices=None):
834	"""Return maximum value of quantity (on centroids)
835
836	Optional argument:
837	indices is the set of element ids that the operation applies to.
838
839	Usage:
840	v = get_maximum_value()
841
842	Note, we do not seek the maximum at vertices as each vertex can
843	have multiple values - one for each triangle sharing it
844	"""
845
846
847	i = self.get_maximum_index(indices)
848	V = self.get_values(location='centroids') #, indices=indices)
849
850	return V[i]
851
852
853	def get_maximum_location(self, indices=None):
854	"""Return location of maximum value of quantity (on centroids)
855
856	Optional argument:
857	indices is the set of element ids that the operation applies to.
858
859	Usage:
860	x, y = get_maximum_location()
861
862
863	Notes:
864	We do not seek the maximum at vertices as each vertex can
865	have multiple values - one for each triangle sharing it.
866
867	If there are multiple cells with same maximum value, the
868	first cell encountered in the triangle array is returned.
869	"""
870
871	i = self.get_maximum_index(indices)
872	x, y = self.domain.get_centroid_coordinates()[i]
873
874	return x, y
875
876
877	def get_minimum_index(self, indices=None):
878	"""See get extreme index for details
879	"""
880
881	return self.get_extremum_index(mode='min',
882	indices=indices)
883
884
885	def get_minimum_value(self, indices=None):
886	"""Return minimum value of quantity (on centroids)
887
888	Optional argument:
889	indices is the set of element ids that the operation applies to.
890
891	Usage:
892	v = get_minimum_value()
893
894	See get_maximum_value for more details.
895	"""
896
897
898	i = self.get_minimum_index(indices)
899	V = self.get_values(location='centroids')
900
901	return V[i]
902
903
904	def get_minimum_location(self, indices=None):
905	"""Return location of 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	x, y = get_minimum_location()
912
913
914	Notes:
915	We do not seek the maximum at vertices as each vertex can
916	have multiple values - one for each triangle sharing it.
917
918	If there are multiple cells with same maximum value, the
919	first cell encountered in the triangle array is returned.
920	"""
921
922	i = self.get_minimum_index(indices)
923	x, y = self.domain.get_centroid_coordinates()[i]
924
925	return x, y
926
927
928
929
930	def get_interpolated_values(self, interpolation_points):
931
932	# Interpolation object based on internal (discontinuous triangles)
933	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
934	smooth=False)
935	# FIXME: This concat should roll into get_vertex_values
936	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
937	axis=1)
938
939	can_reuse = False
940	if hasattr(self, 'interpolation_object'):
941	# Reuse to save time
942	I = self.interpolation_object
943
944	if allclose(interpolation_points, I._point_coordinates):
945	can_reuse = True
946
947
948	if can_reuse is True:
949	# Use absence of points to indicate reuse in I.interpolate
950	result = I.interpolate(vertex_values)
951	else:
952	from anuga.fit_interpolate.interpolate import Interpolate
953
954	# Create interpolation object with matrix
955	I = Interpolate(vertex_coordinates, triangles)
956	self.interpolation_object = I
957
958	# Call interpolate with points the first time
959	interpolation_points = ensure_numeric(interpolation_points, Float)
960	result = I.interpolate(vertex_values, interpolation_points)
961
962	return result
963
964
965	def get_values(self, interpolation_points=None,
966	location='vertices',
967	indices = None):
968	"""get values for quantity
969
970	return X, Compatible list, Numeric array (see below)
971	interpolation_points: List of x, y coordinates where value is
972	sought (using interpolation). If points are given, values of
973	location and indices are ignored
974
975	location: Where values are to be stored.
976	Permissible options are: vertices, edges, centroids
977	and unique vertices. Default is 'vertices'
978
979
980	The returned values with be a list the length of indices
981	(N if indices = None).
982
983	In case of location == 'centroids' the dimension of returned
984	values will be a list or a Numerical array of length N, N being
985	the number of elements.
986
987	In case of location == 'vertices' or 'edges' the dimension of
988	returned values will be of dimension Nx3
989
990	In case of location == 'unique vertices' the average value at
991	each vertex will be returned and the dimension of returned values
992	will be a 1d array of length "number of vertices"
993
994	Indices is the set of element ids that the operation applies to.
995
996	The values will be stored in elements following their
997	internal ordering.
998
999	"""
1000	from Numeric import take
1001
1002	if interpolation_points is not None:
1003	return self.get_interpolated_values(interpolation_points)
1004
1005
1006
1007	if location not in ['vertices', 'centroids', 'edges',
1008	'unique vertices']:
1009	msg = 'Invalid location: %s' %location
1010	raise msg
1011
1012	import types, Numeric
1013	assert type(indices) in [types.ListType, types.NoneType,
1014	Numeric.ArrayType],\
1015	'Indices must be a list or None'
1016
1017	if location == 'centroids':
1018	if (indices == None):
1019	indices = range(len(self))
1020	return take(self.centroid_values,indices)
1021	elif location == 'edges':
1022	if (indices == None):
1023	indices = range(len(self))
1024	return take(self.edge_values,indices)
1025	elif location == 'unique vertices':
1026	if (indices == None):
1027	indices=range(self.domain.number_of_nodes)
1028	vert_values = []
1029	#Go through list of unique vertices
1030	for unique_vert_id in indices:
1031	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1032
1033	#In case there are unused points
1034	if len(triangles) == 0:
1035	msg = 'Unique vertex not associated with triangles'
1036	raise msg
1037
1038	# Go through all triangle, vertex pairs
1039	# Average the values
1040
1041	# FIXME (Ole): Should we merge this with get_vertex_values
1042	sum = 0
1043	for triangle_id, vertex_id in triangles:
1044	sum += self.vertex_values[triangle_id, vertex_id]
1045	vert_values.append(sum/len(triangles))
1046	return Numeric.array(vert_values)
1047	else:
1048	if (indices is None):
1049	indices = range(len(self))
1050	return take(self.vertex_values, indices)
1051
1052
1053
1054	def set_vertex_values(self, A, indices = None):
1055	"""Set vertex values for all unique vertices based on input array A
1056	which has one entry per unique vertex, i.e.
1057	one value for each row in array self.domain.nodes.
1058
1059	indices is the list of vertex_id's that will be set.
1060
1061	This function is used by set_values_from_array
1062	"""
1063
1064	from Numeric import array, Float
1065
1066	#Assert that A can be converted to a Numeric array of appropriate dim
1067	A = array(A, Float)
1068
1069	#print 'SHAPE A', A.shape
1070	assert len(A.shape) == 1
1071
1072	if indices is None:
1073	assert A.shape[0] == self.domain.get_nodes().shape[0]
1074	vertex_list = range(A.shape[0])
1075	else:
1076	assert A.shape[0] == len(indices)
1077	vertex_list = indices
1078
1079	#Go through list of unique vertices
1080
1081	for i_index, unique_vert_id in enumerate(vertex_list):
1082
1083
1084	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1085
1086	#In case there are unused points
1087	if len(triangles) == 0: continue
1088
1089	#Go through all triangle, vertex pairs
1090	#touching vertex unique_vert_id and set corresponding vertex value
1091	for triangle_id, vertex_id in triangles:
1092	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1093
1094	#Intialise centroid and edge_values
1095	self.interpolate()
1096
1097
1098	def smooth_vertex_values(self, value_array='field_values',
1099	precision = None):
1100	""" Smooths field_values or conserved_quantities data.
1101	TODO: be able to smooth individual fields
1102	NOTE: This function does not have a test.
1103	FIXME: NOT DONE - do we need it?
1104	FIXME: this function isn't called by anything.
1105	Maybe it should be removed..-DSG
1106	"""
1107
1108	from Numeric import concatenate, zeros, Float, Int, array, reshape
1109
1110
1111	A,V = self.get_vertex_values(xy=False,
1112	value_array=value_array,
1113	smooth = True,
1114	precision = precision)
1115
1116	#Set some field values
1117	for volume in self:
1118	for i,v in enumerate(volume.vertices):
1119	if value_array == 'field_values':
1120	volume.set_field_values('vertex', i, A[v,:])
1121	elif value_array == 'conserved_quantities':
1122	volume.set_conserved_quantities('vertex', i, A[v,:])
1123
1124	if value_array == 'field_values':
1125	self.precompute()
1126	elif value_array == 'conserved_quantities':
1127	Volume.interpolate_conserved_quantities()
1128
1129
1130	# Methods for outputting model results
1131	def get_vertex_values(self,
1132	xy=True,
1133	smooth=None,
1134	precision=None):
1135	"""Return vertex values like an OBJ format i.e. one value per node.
1136
1137	The vertex values are returned as one sequence in the 1D float array A.
1138	If requested the coordinates will be returned in 1D arrays X and Y.
1139
1140	The connectivity is represented as an integer array, V, of dimension
1141	Mx3, where M is the number of triangles. Each row has three indices
1142	defining the triangle and they correspond to elements in the arrays
1143	X, Y and A.
1144
1145	if smooth is True, vertex values corresponding to one common
1146	coordinate set will be smoothed by taking the average of vertex values for each node.
1147	In this case vertex coordinates will be
1148	de-duplicated corresponding to the original nodes as obtained from
1149	the method general_mesh.get_nodes()
1150
1151	If no smoothings is required, vertex coordinates and values will
1152	be aggregated as a concatenation of values at
1153	vertices 0, vertices 1 and vertices 2. This corresponds to
1154	the node coordinates obtained from the method
1155	general_mesh.get_vertex_coordinates()
1156
1157
1158	Calling convention
1159	if xy is True:
1160	X,Y,A,V = get_vertex_values
1161	else:
1162	A,V = get_vertex_values
1163
1164	"""
1165
1166	from Numeric import concatenate, zeros, Float, Int, array, reshape
1167
1168
1169	if smooth is None:
1170	# Take default from domain
1171	smooth = self.domain.smooth
1172
1173	if precision is None:
1174	precision = Float
1175
1176
1177	if smooth is True:
1178	# Ensure continuous vertex values by averaging
1179	# values at each node
1180
1181	V = self.domain.get_triangles()
1182	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1183	A = zeros(N, Float)
1184	points = self.domain.get_nodes()
1185
1186	if 1:
1187	# Fast C version
1188	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1189	ensure_numeric(self.domain.number_of_triangles_per_node),
1190	ensure_numeric(self.vertex_values),
1191	A)
1192	A = A.astype(precision)
1193	else:
1194
1195	# Slow Python version
1196
1197	current_node = 0
1198	k = 0 # Track triangles touching on node
1199	total = 0.0
1200	for index in self.domain.vertex_value_indices:
1201	if current_node == N:
1202	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1203	raise msg
1204
1205
1206
1207	k += 1
1208
1209	volume_id = index / 3
1210	vertex_id = index % 3
1211
1212	#assert V[volume_id, vertex_id] == current_node
1213
1214	v = self.vertex_values[volume_id, vertex_id]
1215	total += v
1216
1217	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1218	if self.domain.number_of_triangles_per_node[current_node] == k:
1219	A[current_node] = total/k
1220
1221
1222	# Move on to next node
1223	total = 0.0
1224	k = 0
1225	current_node += 1
1226
1227
1228
1229	else:
1230	# Allow discontinuous vertex values
1231	V = self.domain.get_disconnected_triangles()
1232	points = self.domain.get_vertex_coordinates()
1233	A = self.vertex_values.flat.astype(precision)
1234
1235
1236	# Return
1237	if xy is True:
1238	X = points[:,0].astype(precision)
1239	Y = points[:,1].astype(precision)
1240
1241	return X, Y, A, V
1242	else:
1243	return A, V
1244
1245
1246
1247	def extrapolate_first_order(self):
1248	"""Extrapolate conserved quantities from centroid to
1249	vertices for each volume using
1250	first order scheme.
1251	"""
1252
1253	qc = self.centroid_values
1254	qv = self.vertex_values
1255
1256	for i in range(3):
1257	qv[:,i] = qc
1258
1259
1260	def get_integral(self):
1261	"""Compute the integral of quantity across entire domain
1262	"""
1263	integral = 0
1264	for k in range(len(self.domain)):
1265	area = self.domain.areas[k]
1266	qc = self.centroid_values[k]
1267	integral += qc*area
1268
1269	return integral
1270
1271
1272
1273
1274	class Conserved_quantity(Quantity):
1275	"""Class conserved quantity adds to Quantity:
1276
1277	boundary values, storage and method for updating, and
1278	methods for (second order) extrapolation from centroid to vertices inluding
1279	gradients and limiters
1280	"""
1281
1282	def __init__(self, domain, vertex_values=None):
1283	Quantity.__init__(self, domain, vertex_values)
1284
1285	from Numeric import zeros, Float
1286
1287	#Allocate space for boundary values
1288	L = len(domain.boundary)
1289	self.boundary_values = zeros(L, Float)
1290
1291	#Allocate space for updates of conserved quantities by
1292	#flux calculations and forcing functions
1293
1294	N = len(domain) # number_of_triangles
1295	self.explicit_update = zeros(N, Float )
1296	self.semi_implicit_update = zeros(N, Float )
1297	self.centroid_backup_values = zeros(N, Float)
1298
1299
1300
1301
1302	def update(self, timestep):
1303	#Call correct module function
1304	#(either from this module or C-extension)
1305	return update(self, timestep)
1306
1307
1308	def compute_gradients(self):
1309	#Call correct module function
1310	#(either from this module or C-extension)
1311	return compute_gradients(self)
1312
1313	def limit(self):
1314	#Call correct module function
1315	#(either from this module or C-extension)
1316	limit(self)
1317
1318	def extrapolate_second_order(self):
1319	#Call correct module function
1320	#(either from this module or C-extension)
1321	extrapolate_second_order(self)
1322
1323	def backup_centroid_values(self):
1324	#Call correct module function
1325	#(either from this module or C-extension)
1326	backup_centroid_values(self)
1327
1328	def saxpy_centroid_values(self,a,b):
1329	#Call correct module function
1330	#(either from this module or C-extension)
1331	saxpy_centroid_values(self,a,b)
1332
1333
1334	def update(quantity, timestep):
1335	"""Update centroid values based on values stored in
1336	explicit_update and semi_implicit_update as well as given timestep
1337
1338	Function implementing forcing terms must take on argument
1339	which is the domain and they must update either explicit
1340	or implicit updates, e,g,:
1341
1342	def gravity(domain):
1343	....
1344	domain.quantities['xmomentum'].explicit_update = ...
1345	domain.quantities['ymomentum'].explicit_update = ...
1346
1347
1348
1349	Explicit terms must have the form
1350
1351	G(q, t)
1352
1353	and explicit scheme is
1354
1355	q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
1356
1357
1358	Semi implicit forcing terms are assumed to have the form
1359
1360	G(q, t) = H(q, t) q
1361
1362	and the semi implicit scheme will then be
1363
1364	q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
1365
1366
1367	"""
1368
1369	from Numeric import sum, equal, ones, exp, Float
1370
1371	N = quantity.centroid_values.shape[0]
1372
1373
1374	# Divide H by conserved quantity to obtain G (see docstring above)
1375
1376
1377	for k in range(N):
1378	x = quantity.centroid_values[k]
1379	if x == 0.0:
1380	#FIXME: Is this right
1381	quantity.semi_implicit_update[k] = 0.0
1382	else:
1383	quantity.semi_implicit_update[k] /= x
1384
1385
1386	# Semi implicit updates
1387	denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1388
1389	if sum(less(denominator, 1.0)) > 0.0:
1390	msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
1391	raise msg
1392
1393	if sum(equal(denominator, 0.0)) > 0.0:
1394	msg = 'Zero division in semi implicit update. Call Stephen :-)'
1395	raise msg
1396	else:
1397	#Update conserved_quantities from semi implicit updates
1398	quantity.centroid_values /= denominator
1399
1400	# quantity.centroid_values = exp(timestepquantity.semi_implicit_update)quantity.centroid_values
1401
1402	#Explicit updates
1403	quantity.centroid_values += timestep*quantity.explicit_update
1404
1405	def interpolate_from_vertices_to_edges(quantity):
1406	"""Compute edge values from vertex values using linear interpolation
1407	"""
1408
1409	for k in range(quantity.vertex_values.shape[0]):
1410	q0 = quantity.vertex_values[k, 0]
1411	q1 = quantity.vertex_values[k, 1]
1412	q2 = quantity.vertex_values[k, 2]
1413
1414	quantity.edge_values[k, 0] = 0.5*(q1+q2)
1415	quantity.edge_values[k, 1] = 0.5*(q0+q2)
1416	quantity.edge_values[k, 2] = 0.5*(q0+q1)
1417
1418
1419
1420	def backup_centroid_values(quantity):
1421	"""Copy centroid values to backup array"""
1422
1423	qc = quantity.centroid_values
1424	qb = quantity.centroid_backup_values
1425
1426	#Check each triangle
1427	for k in range(len(quantity.domain)):
1428	qb[k] = qc[k]
1429
1430
1431	def saxpy_centroid_values(quantity,a,b):
1432	"""saxpy operation between centroid value and backup"""
1433
1434	qc = quantity.centroid_values
1435	qb = quantity.centroid_backup_values
1436
1437	#Check each triangle
1438	for k in range(len(quantity.domain)):
1439	qc[k] = aqc[k]+bqb[k]
1440
1441
1442	def compute_gradients(quantity):
1443	"""Compute gradients of triangle surfaces defined by centroids of
1444	neighbouring volumes.
1445	If one edge is on the boundary, use own centroid as neighbour centroid.
1446	If two or more are on the boundary, fall back to first order scheme.
1447	"""
1448
1449	from Numeric import zeros, Float
1450	from utilitites.numerical_tools import gradient
1451
1452	centroid_coordinates = quantity.domain.centroid_coordinates
1453	surrogate_neighbours = quantity.domain.surrogate_neighbours
1454	centroid_values = quantity.centroid_values
1455	number_of_boundaries = quantity.domain.number_of_boundaries
1456
1457	N = centroid_values.shape[0]
1458
1459	a = zeros(N, Float)
1460	b = zeros(N, Float)
1461
1462	for k in range(N):
1463	if number_of_boundaries[k] < 2:
1464	#Two or three true neighbours
1465
1466	#Get indices of neighbours (or self when used as surrogate)
1467	k0, k1, k2 = surrogate_neighbours[k,:]
1468
1469	#Get data
1470	q0 = centroid_values[k0]
1471	q1 = centroid_values[k1]
1472	q2 = centroid_values[k2]
1473
1474	x0, y0 = centroid_coordinates[k0] #V0 centroid
1475	x1, y1 = centroid_coordinates[k1] #V1 centroid
1476	x2, y2 = centroid_coordinates[k2] #V2 centroid
1477
1478	#Gradient
1479	a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1480
1481	elif number_of_boundaries[k] == 2:
1482	#One true neighbour
1483
1484	#Get index of the one neighbour
1485	for k0 in surrogate_neighbours[k,:]:
1486	if k0 != k: break
1487	assert k0 != k
1488
1489	k1 = k #self
1490
1491	#Get data
1492	q0 = centroid_values[k0]
1493	q1 = centroid_values[k1]
1494
1495	x0, y0 = centroid_coordinates[k0] #V0 centroid
1496	x1, y1 = centroid_coordinates[k1] #V1 centroid
1497
1498	#Gradient
1499	a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1500	else:
1501	#No true neighbours -
1502	#Fall back to first order scheme
1503	pass
1504
1505
1506	return a, b
1507
1508
1509
1510	def limit(quantity):
1511	"""Limit slopes for each volume to eliminate artificial variance
1512	introduced by e.g. second order extrapolator
1513
1514	This is an unsophisticated limiter as it does not take into
1515	account dependencies among quantities.
1516
1517	precondition:
1518	vertex values are estimated from gradient
1519	postcondition:
1520	vertex values are updated
1521	"""
1522
1523	from Numeric import zeros, Float
1524
1525	N = quantity.domain.number_of_nodes
1526
1527	beta_w = quantity.domain.beta_w
1528
1529	qc = quantity.centroid_values
1530	qv = quantity.vertex_values
1531
1532	#Find min and max of this and neighbour's centroid values
1533	qmax = zeros(qc.shape, Float)
1534	qmin = zeros(qc.shape, Float)
1535
1536	for k in range(N):
1537	qmax[k] = qc[k]
1538	qmin[k] = qc[k]
1539	for i in range(3):
1540	n = quantity.domain.neighbours[k,i]
1541	if n >= 0:
1542	qn = qc[n] #Neighbour's centroid value
1543
1544	qmin[k] = min(qmin[k], qn)
1545	qmax[k] = max(qmax[k], qn)
1546	qmax[k] = min(qmax[k], 2.0*qc[k])
1547	qmin[k] = max(qmin[k], 0.5*qc[k])
1548
1549
1550	#Diffences between centroids and maxima/minima
1551	dqmax = qmax - qc
1552	dqmin = qmin - qc
1553
1554	#Deltas between vertex and centroid values
1555	dq = zeros(qv.shape, Float)
1556	for i in range(3):
1557	dq[:,i] = qv[:,i] - qc
1558
1559	#Phi limiter
1560	for k in range(N):
1561
1562	#Find the gradient limiter (phi) across vertices
1563	phi = 1.0
1564	for i in range(3):
1565	r = 1.0
1566	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1567	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1568
1569	phi = min( min(r*beta_w, 1), phi )
1570
1571	#Then update using phi limiter
1572	for i in range(3):
1573	qv[k,i] = qc[k] + phi*dq[k,i]
1574
1575
1576
1577	from anuga.utilities import compile
1578	if compile.can_use_C_extension('quantity_ext.c'):
1579	#Replace python version with c implementations
1580
1581	from quantity_ext import average_vertex_values, backup_centroid_values, \
1582	saxpy_centroid_values
1583
1584	from quantity_ext import compute_gradients, limit,\
1585	extrapolate_second_order, interpolate_from_vertices_to_edges, update

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

Download in other formats: