Context Navigation

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

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

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

Download in other formats: