Context Navigation

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

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

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

Download in other formats: