Context Navigation

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

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

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

Download in other formats: