Context Navigation

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

Last change on this file since 4806 was 4806, checked in by duncan, 16 years ago
Removing fit using the domain mesh. There's a bug somewhere.
File size: 46.5 KB

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

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

Download in other formats: