Context Navigation

source: anuga_core/source_numpy_conversion/anuga/abstract_2d_finite_volumes/quantity.py @ 5899

Last change on this file since 5899 was 5899, checked in by rwilson, 15 years ago
Initial NumPy? changes (again!).
File size: 54.0 KB

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

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

Download in other formats: