Context Navigation

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

Last change on this file since 5736 was 5736, checked in by ole, 15 years ago
Implemented run-time version of get_energy_through_cross_section and tested it.
File size: 51.6 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
1015	list of absolute UTM coordinates or a geospatial data object.
1016	"""
1017
1018	# FIXME (Ole): Could do with an input check (should be generalised
1019	# and used widely)
1020	# That will check that interpolation points is either a list of
1021	# points, Nx2 array, or geospatial
1022
1023	# Ensure points are converted to coordinates relative to mesh origin
1024	# FIXME (Ole): This could all be refactored using the
1025	# 'change_points_geo_ref' method of Class geo_reference.
1026	# The purpose is to make interpolation points relative
1027	# to the mesh origin.
1028	#
1029	# Speed is also a consideration here.
1030
1031	if isinstance(interpolation_points, Geospatial_data):
1032	# Ensure interpolation points are in absolute UTM coordinates
1033	interpolation_points = interpolation_points.get_data_points(absolute=True)
1034
1035	# Reconcile interpolation points with georeference of domain
1036	interpolation_points = self.domain.geo_reference.get_relative(interpolation_points)
1037	interpolation_points = ensure_numeric(interpolation_points)
1038
1039
1040
1041	# Get internal (discontinuous) triangles for use with the
1042	# interpolation object.
1043	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
1044	smooth=False)
1045	# FIXME (Ole): This concat should roll into get_vertex_values
1046	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
1047	axis=1)
1048
1049	can_reuse = False
1050	if hasattr(self, 'interpolation_object'):
1051	# Reuse to save time
1052	I = self.interpolation_object
1053
1054	if allclose(interpolation_points.shape,
1055	I._point_coordinates.shape):
1056	if allclose(interpolation_points, I._point_coordinates):
1057	can_reuse = True
1058
1059
1060	if can_reuse is True:
1061	# Use absence of points to indicate reuse in I.interpolate
1062	result = I.interpolate(vertex_values)
1063	else:
1064	from anuga.fit_interpolate.interpolate import Interpolate
1065
1066	# Create interpolation object with matrix
1067	I = Interpolate(vertex_coordinates, triangles)
1068	self.interpolation_object = I
1069
1070	# Call interpolate with points the first time
1071	interpolation_points = ensure_numeric(interpolation_points, Float)
1072	result = I.interpolate(vertex_values, interpolation_points)
1073
1074	return result
1075
1076
1077	def get_values(self, interpolation_points=None,
1078	location='vertices',
1079	indices = None):
1080	"""get values for quantity
1081
1082	return X, Compatible list, Numeric array (see below)
1083
1084	Inputs:
1085	interpolation_points: List of x, y coordinates where value is
1086	sought (using interpolation). If points
1087	are given, values of location and indices
1088	are ignored. Assume either absolute UTM
1089	coordinates or geospatial data object.
1090
1091	location: Where values are to be stored.
1092	Permissible options are: vertices, edges, centroids
1093	and unique vertices. Default is 'vertices'
1094
1095
1096	The returned values with be a list the length of indices
1097	(N if indices = None).
1098
1099	In case of location == 'centroids' the dimension of returned
1100	values will be a list or a Numerical array of length N, N being
1101	the number of elements.
1102
1103	In case of location == 'vertices' or 'edges' the dimension of
1104	returned values will be of dimension Nx3
1105
1106	In case of location == 'unique vertices' the average value at
1107	each vertex will be returned and the dimension of returned values
1108	will be a 1d array of length "number of vertices"
1109
1110	Indices is the set of element ids that the operation applies to.
1111
1112	The values will be stored in elements following their
1113	internal ordering.
1114	"""
1115
1116	from Numeric import take
1117
1118	# FIXME (Ole): I reckon we should have the option of passing a
1119	# polygon into get_values. The question becomes how
1120	# resulting values should be ordered.
1121
1122	if interpolation_points is not None:
1123	return self.get_interpolated_values(interpolation_points)
1124
1125
1126	# FIXME (Ole): Consider deprecating 'edges' - but not if it is used
1127	# elsewhere in ANUGA.
1128	# Edges have already been deprecated in set_values, see changeset:5521,
1129	# but might be useful in get_values. Any thoughts anyone?
1130
1131	if location not in ['vertices', 'centroids', 'edges',
1132	'unique vertices']:
1133	msg = 'Invalid location: %s' %location
1134	raise msg
1135
1136	import types, Numeric
1137	assert type(indices) in [types.ListType, types.NoneType,
1138	Numeric.ArrayType],\
1139	'Indices must be a list or None'
1140
1141	if location == 'centroids':
1142	if (indices == None):
1143	indices = range(len(self))
1144	return take(self.centroid_values,indices)
1145	elif location == 'edges':
1146	if (indices == None):
1147	indices = range(len(self))
1148	return take(self.edge_values,indices)
1149	elif location == 'unique vertices':
1150	if (indices == None):
1151	indices=range(self.domain.number_of_nodes)
1152	vert_values = []
1153
1154	# Go through list of unique vertices
1155	for unique_vert_id in indices:
1156	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1157
1158	# In case there are unused points
1159	if len(triangles) == 0:
1160	msg = 'Unique vertex not associated with triangles'
1161	raise msg
1162
1163	# Go through all triangle, vertex pairs
1164	# Average the values
1165
1166	# FIXME (Ole): Should we merge this with get_vertex_values
1167	sum = 0
1168	for triangle_id, vertex_id in triangles:
1169	sum += self.vertex_values[triangle_id, vertex_id]
1170	vert_values.append(sum/len(triangles))
1171	return Numeric.array(vert_values)
1172	else:
1173	if (indices is None):
1174	indices = range(len(self))
1175	return take(self.vertex_values, indices)
1176
1177
1178
1179	def set_vertex_values(self, A, indices = None):
1180	"""Set vertex values for all unique vertices based on input array A
1181	which has one entry per unique vertex, i.e.
1182	one value for each row in array self.domain.nodes.
1183
1184	indices is the list of vertex_id's that will be set.
1185
1186	This function is used by set_values_from_array
1187	"""
1188
1189	from Numeric import array, Float
1190
1191	# Assert that A can be converted to a Numeric array of appropriate dim
1192	A = ensure_numeric(A, Float)
1193
1194	# print 'SHAPE A', A.shape
1195	assert len(A.shape) == 1
1196
1197	if indices is None:
1198	assert A.shape[0] == self.domain.get_nodes().shape[0]
1199	vertex_list = range(A.shape[0])
1200	else:
1201	assert A.shape[0] == len(indices)
1202	vertex_list = indices
1203
1204	# Go through list of unique vertices
1205	for i_index, unique_vert_id in enumerate(vertex_list):
1206
1207
1208	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1209
1210	# In case there are unused points
1211	if len(triangles) == 0: continue
1212
1213	# Go through all triangle, vertex pairs
1214	# touching vertex unique_vert_id and set corresponding vertex value
1215	for triangle_id, vertex_id in triangles:
1216	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1217
1218	# Intialise centroid and edge_values
1219	self.interpolate()
1220
1221
1222	def smooth_vertex_values(self):
1223	""" Smooths vertex values.
1224	"""
1225
1226	A,V = self.get_vertex_values(xy=False, smooth=True)
1227	self.set_vertex_values(A)
1228
1229
1230	# Methods for outputting model results
1231	def get_vertex_values(self,
1232	xy=True,
1233	smooth=None,
1234	precision=None):
1235	"""Return vertex values like an OBJ format i.e. one value per node.
1236
1237	The vertex values are returned as one sequence in the 1D float array A.
1238	If requested the coordinates will be returned in 1D arrays X and Y.
1239
1240	The connectivity is represented as an integer array, V, of dimension
1241	Mx3, where M is the number of triangles. Each row has three indices
1242	defining the triangle and they correspond to elements in the arrays
1243	X, Y and A.
1244
1245	if smooth is True, vertex values corresponding to one common
1246	coordinate set will be smoothed by taking the average of vertex values for each node.
1247	In this case vertex coordinates will be
1248	de-duplicated corresponding to the original nodes as obtained from
1249	the method general_mesh.get_nodes()
1250
1251	If no smoothings is required, vertex coordinates and values will
1252	be aggregated as a concatenation of values at
1253	vertices 0, vertices 1 and vertices 2. This corresponds to
1254	the node coordinates obtained from the method
1255	general_mesh.get_vertex_coordinates()
1256
1257
1258	Calling convention
1259	if xy is True:
1260	X,Y,A,V = get_vertex_values
1261	else:
1262	A,V = get_vertex_values
1263
1264	"""
1265
1266	from Numeric import concatenate, zeros, Float, Int, array, reshape
1267
1268
1269	if smooth is None:
1270	# Take default from domain
1271	try:
1272	smooth = self.domain.smooth
1273	except:
1274	smooth = False
1275
1276	if precision is None:
1277	precision = Float
1278
1279
1280	if smooth is True:
1281	# Ensure continuous vertex values by averaging
1282	# values at each node
1283
1284	V = self.domain.get_triangles()
1285	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1286	A = zeros(N, Float)
1287	points = self.domain.get_nodes()
1288
1289	if 1:
1290	# Fast C version
1291	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1292	ensure_numeric(self.domain.number_of_triangles_per_node),
1293	ensure_numeric(self.vertex_values),
1294	A)
1295	A = A.astype(precision)
1296	else:
1297
1298	# Slow Python version
1299
1300	current_node = 0
1301	k = 0 # Track triangles touching on node
1302	total = 0.0
1303	for index in self.domain.vertex_value_indices:
1304	if current_node == N:
1305	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1306	raise msg
1307
1308
1309
1310	k += 1
1311
1312	volume_id = index / 3
1313	vertex_id = index % 3
1314
1315	#assert V[volume_id, vertex_id] == current_node
1316
1317	v = self.vertex_values[volume_id, vertex_id]
1318	total += v
1319
1320	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1321	if self.domain.number_of_triangles_per_node[current_node] == k:
1322	A[current_node] = total/k
1323
1324
1325	# Move on to next node
1326	total = 0.0
1327	k = 0
1328	current_node += 1
1329
1330
1331
1332	else:
1333	# Allow discontinuous vertex values
1334	V = self.domain.get_disconnected_triangles()
1335	points = self.domain.get_vertex_coordinates()
1336	A = self.vertex_values.flat.astype(precision)
1337
1338
1339	# Return
1340	if xy is True:
1341	X = points[:,0].astype(precision)
1342	Y = points[:,1].astype(precision)
1343
1344	return X, Y, A, V
1345	else:
1346	return A, V
1347
1348
1349
1350	def extrapolate_first_order(self):
1351	"""Extrapolate conserved quantities from centroid to
1352	vertices and edges for each volume using
1353	first order scheme.
1354	"""
1355
1356	qc = self.centroid_values
1357	qv = self.vertex_values
1358	qe = self.edge_values
1359
1360	for i in range(3):
1361	qv[:,i] = qc
1362	qe[:,i] = qc
1363
1364	self.x_gradient *= 0.0
1365	self.y_gradient *= 0.0
1366
1367
1368	def get_integral(self):
1369	"""Compute the integral of quantity across entire domain
1370	"""
1371	integral = 0
1372	for k in range(len(self.domain)):
1373	area = self.domain.areas[k]
1374	qc = self.centroid_values[k]
1375	integral += qc*area
1376
1377	return integral
1378
1379	def get_gradients(self):
1380	"""Provide gradients. Use compute_gradients first
1381	"""
1382
1383	return self.x_gradient, self.y_gradient
1384
1385
1386	def update(self, timestep):
1387	# Call correct module function
1388	# (either from this module or C-extension)
1389	return update(self, timestep)
1390
1391	def compute_gradients(self):
1392	# Call correct module function
1393	# (either from this module or C-extension)
1394	return compute_gradients(self)
1395
1396	def limit(self):
1397	# Call correct module depending on whether
1398	# basing limit calculations on edges or vertices
1399	limit_old(self)
1400
1401	def limit_vertices_by_all_neighbours(self):
1402	# Call correct module function
1403	# (either from this module or C-extension)
1404	limit_vertices_by_all_neighbours(self)
1405
1406	def limit_edges_by_all_neighbours(self):
1407	# Call correct module function
1408	# (either from this module or C-extension)
1409	limit_edges_by_all_neighbours(self)
1410
1411	def limit_edges_by_neighbour(self):
1412	# Call correct module function
1413	# (either from this module or C-extension)
1414	limit_edges_by_neighbour(self)
1415
1416	def extrapolate_second_order(self):
1417	# Call correct module function
1418	# (either from this module or C-extension)
1419	compute_gradients(self)
1420	extrapolate_from_gradient(self)
1421
1422	def extrapolate_second_order_and_limit_by_edge(self):
1423	# Call correct module function
1424	# (either from this module or C-extension)
1425	extrapolate_second_order_and_limit_by_edge(self)
1426
1427	def extrapolate_second_order_and_limit_by_vertex(self):
1428	# Call correct module function
1429	# (either from this module or C-extension)
1430	extrapolate_second_order_and_limit_by_vertex(self)
1431
1432	def bound_vertices_below_by_constant(self, bound):
1433	# Call correct module function
1434	# (either from this module or C-extension)
1435	bound_vertices_below_by_constant(self, bound)
1436
1437	def bound_vertices_below_by_quantity(self, quantity):
1438	# Call correct module function
1439	# (either from this module or C-extension)
1440
1441	# check consistency
1442	assert self.domain == quantity.domain
1443	bound_vertices_below_by_quantity(self, quantity)
1444
1445	def backup_centroid_values(self):
1446	# Call correct module function
1447	# (either from this module or C-extension)
1448	backup_centroid_values(self)
1449
1450	def saxpy_centroid_values(self,a,b):
1451	# Call correct module function
1452	# (either from this module or C-extension)
1453	saxpy_centroid_values(self,a,b)
1454
1455	#Conserved_quantity = Quantity
1456
1457	class Conserved_quantity(Quantity):
1458	"""Class conserved quantity being removed, use Quantity
1459
1460	"""
1461
1462	def __init__(self, domain, vertex_values=None):
1463	#Quantity.__init__(self, domain, vertex_values)
1464
1465	msg = 'ERROR: Use Quantity instead of Conserved_quantity'
1466
1467	raise Exception, msg
1468
1469
1470
1471	from anuga.utilities import compile
1472	if compile.can_use_C_extension('quantity_ext.c'):
1473	# Underlying C implementations can be accessed
1474
1475	from quantity_ext import \
1476	average_vertex_values,\
1477	backup_centroid_values,\
1478	saxpy_centroid_values,\
1479	compute_gradients,\
1480	limit_old,\
1481	limit_vertices_by_all_neighbours,\
1482	limit_edges_by_all_neighbours,\
1483	limit_edges_by_neighbour,\
1484	limit_gradient_by_neighbour,\
1485	extrapolate_from_gradient,\
1486	extrapolate_second_order_and_limit_by_edge,\
1487	extrapolate_second_order_and_limit_by_vertex,\
1488	bound_vertices_below_by_constant,\
1489	bound_vertices_below_by_quantity,\
1490	interpolate_from_vertices_to_edges,\
1491	interpolate_from_edges_to_vertices,\
1492	update
1493	else:
1494	msg = 'C implementations could not be accessed by %s.\n ' %__file__
1495	msg += 'Make sure compile_all.py has been run as described in '
1496	msg += 'the ANUGA installation guide.'
1497	raise Exception, msg
1498
1499
1500

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

Download in other formats: