Context Navigation

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

Last change on this file since 5162 was 5162, checked in by steve, 16 years ago
Updated some methods for quantity. Looks like we can use old limiting system with larger values of beta.
File size: 50.4 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	# New leaner interface to setting values
262	def set_values(self,
263	numeric = None, # List, numeric array or constant
264	quantity = None, # Another quantity
265	function = None, # Callable object: f(x,y)
266	geospatial_data = None, # Arbitrary dataset
267	points = None, values = None, data_georef = None, # Obsoleted by use of geo_spatial object
268	filename = None, attribute_name = None, #Input from file
269	alpha = None,
270	location = 'vertices',
271	polygon = None,
272	indices = None,
273	smooth = False,
274	verbose = False,
275	use_cache = False):
276
277	"""Set values for quantity based on different sources.
278
279	numeric:
280	Compatible list, Numeric array (see below) or constant.
281	If callable it will treated as a function (see below)
282	If instance of another Quantity it will be treated as such.
283	If geo_spatial object it will be treated as such
284
285	quantity:
286	Another quantity (compatible quantity, e.g. obtained as a
287	linear combination of quantities)
288
289	function:
290	Any callable object that takes two 1d arrays x and y
291	each of length N and returns an array also of length N.
292	The function will be evaluated at points determined by
293	location and indices in the underlying mesh.
294
295	geospatial_data:
296	Arbitrary geo spatial dataset in the form of the class
297	Geospatial_data. Mesh points are populated using
298	fit_interpolate.fit fitting
299
300	points:
301	Nx2 array of data points for use with fit_interpolate.fit
302	If points are present, an N array of attribute
303	values corresponding to
304	each data point must be present.
305	(Obsoleted by geospatial_data)
306
307	values:
308	If points is specified, values is an array of length N containing
309	attribute values for each point.
310	(Obsoleted by geospatial_data)
311
312	data_georef:
313	If points is specified, geo_reference applies to each point.
314	(Obsoleted by geospatial_data)
315
316	filename:
317	Name of a points file containing data points and attributes for
318	use with fit_interpolate.fit.
319
320	attribute_name:
321	If specified, any array matching that name
322	will be used. from file or geospatial_data.
323	Otherwise a default will be used.
324
325	alpha:
326	Smoothing parameter to be used with fit_interpolate.fit.
327	See module fit_interpolate.fit for further details about alpha.
328	Alpha will only be used with points, values or filename.
329	Otherwise it will be ignored.
330
331
332	location: Where values are to be stored.
333	Permissible options are: vertices, edges, centroids
334	Default is 'vertices'
335
336	In case of location == 'centroids' the dimension values must
337	be a list of a Numerical array of length N,
338	N being the number of elements.
339	Otherwise it must be of dimension Nx3
340
341
342	The values will be stored in elements following their
343	internal ordering.
344
345	If location is not 'unique vertices' Indices is the
346	set of element ids that the operation applies to.
347	If location is 'unique vertices' Indices is the set
348	of vertex ids that the operation applies to.
349
350	If selected location is vertices, values for
351	centroid and edges will be assigned interpolated
352	values. In any other case, only values for the
353	specified locations will be assigned and the others
354	will be left undefined.
355
356
357	polygon: Restrict update of quantity to locations that fall
358	inside polygon. Polygon works by selecting indices
359	and calling set_values recursively.
360	Polygon mode has only been implemented for
361	constant values so far.
362
363	indices: Restrict update of quantity to locations that are
364	identified by indices (e.g. node ids if location
365	is 'vertices')
366
367	verbose: True means that output to stdout is generated
368
369	use_cache: True means that caching of intermediate results is
370	attempted for fit_interpolate.fit.
371
372
373
374
375	Exactly one of the arguments
376	numeric, quantity, function, points, filename
377	must be present.
378	"""
379
380	from anuga.geospatial_data.geospatial_data import Geospatial_data
381	from types import FloatType, IntType, LongType, ListType, NoneType
382	from Numeric import ArrayType
383
384
385	# Treat special case: Polygon situation
386	# Location will be ignored and set to 'centroids'
387	# FIXME (Ole): This needs to be generalised and
388	# perhaps the notion of location and indices simplified
389
390	if polygon is not None:
391	if indices is not None:
392	msg = 'Only one of polygon and indices can be specified'
393	raise Exception, msg
394
395	msg = 'With polygon selected, set_quantity must provide '
396	msg += 'the keyword numeric and it must (currently) be '
397	msg += 'a constant.'
398	if numeric is None:
399	raise Exception, msg
400	else:
401	# Check that numeric is as constant
402	assert type(numeric) in [FloatType, IntType, LongType], msg
403
404
405	location = 'centroids'
406
407
408	points = self.domain.get_centroid_coordinates(absolute=True)
409	indices = inside_polygon(points, polygon)
410
411	self.set_values_from_constant(numeric,
412	location, indices, verbose)
413
414
415	self.extrapolate_first_order()
416
417	if smooth:
418	self.smooth_vertex_values()
419
420
421	return
422
423
424
425
426
427
428	# General input checks
429	L = [numeric, quantity, function, geospatial_data, points, filename]
430	msg = 'Exactly one of the arguments '+\
431	'numeric, quantity, function, geospatial_data, points, '+\
432	'or filename must be present.'
433	assert L.count(None) == len(L)-1, msg
434
435
436	if location not in ['vertices', 'centroids', 'edges',
437	'unique vertices']:
438	msg = 'Invalid location: %s' %location
439	raise Exception, msg
440
441
442	msg = 'Indices must be a list or None'
443	assert type(indices) in [ListType, NoneType, ArrayType], msg
444
445
446	if not(points is None and values is None and data_georef is None):
447	from warnings import warn
448
449	msg = 'Using points, values or data_georef with set_quantity '
450	msg += 'is obsolete. Please use a Geospatial_data object instead.'
451	#warn(msg, DeprecationWarning)
452	raise Exception, msg
453
454
455
456	# Determine which 'set_values_from_...' to use
457
458	if numeric is not None:
459	if type(numeric) in [FloatType, IntType, LongType]:
460	self.set_values_from_constant(numeric,
461	location, indices, verbose)
462	elif type(numeric) in [ArrayType, ListType]:
463	self.set_values_from_array(numeric,
464	location, indices, verbose)
465	elif callable(numeric):
466	self.set_values_from_function(numeric,
467	location, indices, verbose)
468	elif isinstance(numeric, Quantity):
469	self.set_values_from_quantity(numeric,
470	location, indices, verbose)
471	elif isinstance(numeric, Geospatial_data):
472	self.set_values_from_geospatial_data(numeric,
473	alpha,
474	location, indices,
475	verbose = verbose,
476	use_cache = use_cache)
477	else:
478	msg = 'Illegal type for argument numeric: %s' %str(numeric)
479	raise msg
480
481	elif quantity is not None:
482	self.set_values_from_quantity(quantity,
483	location, indices, verbose)
484	elif function is not None:
485	msg = 'Argument function must be callable'
486	assert callable(function), msg
487	self.set_values_from_function(function,
488	location, indices, verbose)
489	elif geospatial_data is not None:
490	self.set_values_from_geospatial_data(geospatial_data,
491	alpha,
492	location, indices,
493	verbose = verbose,
494	use_cache = use_cache)
495	elif points is not None:
496	msg = 'The usage of points in set_values has been deprecated.' +\
497	'Please use the geospatial_data object instead.'
498	raise Exception, msg
499
500
501
502	elif filename is not None:
503	if hasattr(self.domain, 'points_file_block_line_size'):
504	max_read_lines = self.domain.points_file_block_line_size
505	else:
506	max_read_lines = default_block_line_size
507	self.set_values_from_file(filename, attribute_name, alpha,
508	location, indices,
509	verbose = verbose,
510	max_read_lines=max_read_lines,
511	use_cache = use_cache)
512	else:
513	raise Exception, 'This can\'t happen :-)'
514
515
516
517	# Update all locations in triangles
518	if location == 'vertices' or location == 'unique vertices':
519	# Intialise centroid and edge_values
520	self.interpolate()
521
522	if location == 'centroids':
523	# Extrapolate 1st order - to capture notion of area being specified
524	self.extrapolate_first_order()
525
526
527
528	#Specific functions for setting values
529	def set_values_from_constant(self, X,
530	location, indices, verbose):
531	"""Set quantity values from specified constant X
532	"""
533
534	# FIXME (Ole): Somehow indices refer to centroids
535	# rather than vertices as default. See unit test
536	# test_set_vertex_values_using_general_interface_with_subset(self):
537
538
539	if location == 'centroids':
540	if indices is None:
541	self.centroid_values[:] = X
542	else:
543	#Brute force
544	for i in indices:
545	self.centroid_values[i] = X
546
547	elif location == 'edges':
548	if indices is None:
549	self.edge_values[:] = X
550	else:
551	#Brute force
552	for i in indices:
553	self.edge_values[i] = X
554
555	elif location == 'unique vertices':
556	if indices is None:
557	self.edge_values[:] = X #FIXME (Ole): Shouldn't this be vertex_values?
558	else:
559
560	#Go through list of unique vertices
561	for unique_vert_id in indices:
562
563	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
564
565	#In case there are unused points
566	if len(triangles) == 0:
567	continue
568
569	#Go through all triangle, vertex pairs
570	#and set corresponding vertex value
571	for triangle_id, vertex_id in triangles:
572	self.vertex_values[triangle_id, vertex_id] = X
573
574	#Intialise centroid and edge_values
575	self.interpolate()
576	else:
577	if indices is None:
578	self.vertex_values[:] = X
579	else:
580	#Brute force
581	for i_vertex in indices:
582	self.vertex_values[i_vertex] = X
583
584
585
586
587	def set_values_from_array(self, values,
588	location='vertices',
589	indices=None,
590	verbose=False):
591	"""Set values for quantity
592
593	values: Numeric array
594	location: Where values are to be stored.
595	Permissible options are: vertices, edges, centroid, unique vertices
596	Default is 'vertices'
597
598	indices - if this action is carried out on a subset of
599	elements or unique vertices
600	The element/unique vertex indices are specified here.
601
602	In case of location == 'centroid' the dimension values must
603	be a list of a Numerical array of length N, N being the number
604	of elements.
605
606	Otherwise it must be of dimension Nx3
607
608	The values will be stored in elements following their
609	internal ordering.
610
611	If selected location is vertices, values for centroid and edges
612	will be assigned interpolated values.
613	In any other case, only values for the specified locations
614	will be assigned and the others will be left undefined.
615	"""
616
617	from Numeric import array, Float, Int, allclose
618
619	values = array(values).astype(Float)
620
621	if indices is not None:
622	indices = array(indices).astype(Int)
623	msg = 'Number of values must match number of indices:'
624	msg += 'You specified %d values and %d indices'\
625	%(values.shape[0], indices.shape[0])
626	assert values.shape[0] == indices.shape[0], msg
627
628	N = self.centroid_values.shape[0]
629
630	if location == 'centroids':
631	assert len(values.shape) == 1, 'Values array must be 1d'
632
633	if indices is None:
634	msg = 'Number of values must match number of elements'
635	assert values.shape[0] == N, msg
636
637	self.centroid_values = values
638	else:
639	msg = 'Number of values must match number of indices'
640	assert values.shape[0] == indices.shape[0], msg
641
642	#Brute force
643	for i in range(len(indices)):
644	self.centroid_values[indices[i]] = values[i]
645
646	elif location == 'edges':
647	# FIXME (Ole): No mention of indices here. However, I don't
648	# think we ever need to set values at edges anyway
649	assert len(values.shape) == 2, 'Values array must be 2d'
650
651	msg = 'Number of values must match number of elements'
652	assert values.shape[0] == N, msg
653
654	msg = 'Array must be N x 3'
655	assert values.shape[1] == 3, msg
656
657	self.edge_values = values
658
659	elif location == 'unique vertices':
660	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
661	'Values array must be 1d'
662
663	self.set_vertex_values(values.flat, indices=indices)
664
665	else:
666	# Location vertices
667	if len(values.shape) == 1:
668	self.set_vertex_values(values, indices=indices)
669
670	elif len(values.shape) == 2:
671	#Vertex values are given as a triplet for each triangle
672
673	msg = 'Array must be N x 3'
674	assert values.shape[1] == 3, msg
675
676	if indices is None:
677	self.vertex_values = values
678	else:
679	for element_index, value in map(None, indices, values):
680	self.vertex_values[element_index] = value
681	else:
682	msg = 'Values array must be 1d or 2d'
683	raise msg
684
685
686	def set_values_from_quantity(self, q,
687	location, indices, verbose):
688	"""Set quantity values from specified quantity instance q
689
690	Location is ignored - vertices will always be used here.
691	"""
692
693
694	A = q.vertex_values
695
696	from Numeric import allclose
697	msg = 'Quantities are defined on different meshes. '+\
698	'This might be a case for implementing interpolation '+\
699	'between different meshes.'
700	assert allclose(A.shape, self.vertex_values.shape), msg
701
702	self.set_values(A, location='vertices',
703	indices=indices,
704	verbose=verbose)
705
706
707	def set_values_from_function(self, f,
708	location='vertices',
709	indices=None,
710	verbose=False):
711	"""Set values for quantity using specified function
712
713	Input
714
715	f: x, y -> z Function where x, y and z are arrays
716	location: Where values are to be stored.
717	Permissible options are: vertices, centroid, edges,
718	unique vertices
719	Default is "vertices"
720	indices:
721
722
723	"""
724
725	#FIXME: Should check that function returns something sensible and
726	#raise a meaningfull exception if it returns None for example
727
728	#FIXME: Should supply absolute coordinates
729
730
731	# Compute the function values and call set_values again
732	if location == 'centroids':
733	if indices is None:
734	indices = range(len(self))
735
736	V = take(self.domain.get_centroid_coordinates(), indices)
737	self.set_values(f(V[:,0], V[:,1]),
738	location=location,
739	indices=indices)
740
741	elif location == 'vertices':
742
743	M = self.domain.number_of_triangles
744	V = self.domain.get_vertex_coordinates()
745
746	x = V[:,0]; y = V[:,1];
747	values = f(x, y)
748
749
750	# FIXME (Ole): This code should replace all the
751	# rest of this function and it would work, except
752	# one unit test in test_region fails.
753	# If that could be resolved this one will be
754	# more robust and simple.
755
756	#values = reshape(values, (M,3))
757	#self.set_values(values,
758	# location='vertices',
759	# indices=indices)
760
761
762	# This should be removed
763	if is_scalar(values):
764	# Function returned a constant value
765	self.set_values_from_constant(values,
766	location, indices, verbose)
767	return
768
769	# This should be removed
770	if indices is None:
771	for j in range(3):
772	self.vertex_values[:,j] = values[j::3]
773	else:
774	#Brute force
775	for i in indices:
776	for j in range(3):
777	self.vertex_values[i,j] = values[3*i+j]
778
779
780	else:
781	raise 'Not implemented: %s' %location
782
783
784
785	def set_values_from_geospatial_data(self, geospatial_data, alpha,
786	location, indices,
787	verbose = False,
788	use_cache = False):
789	#FIXME: Use this function for the time being. Later move code in here
790
791	points = geospatial_data.get_data_points(absolute = False)
792	values = geospatial_data.get_attributes()
793	data_georef = geospatial_data.get_geo_reference()
794
795
796
797	from anuga.coordinate_transforms.geo_reference import Geo_reference
798
799
800	points = ensure_numeric(points, Float)
801	values = ensure_numeric(values, Float)
802
803	if location != 'vertices':
804	msg = 'set_values_from_points is only defined for '+\
805	'location=\'vertices\''
806	raise ms
807
808	coordinates = self.domain.get_nodes()
809	triangles = self.domain.triangles #FIXME
810
811
812	#Take care of georeferencing
813	if data_georef is None:
814	data_georef = Geo_reference()
815
816
817	mesh_georef = self.domain.geo_reference
818
819
820	# Call fit_interpolate.fit function
821	#args = (coordinates, triangles, points, values)
822	args = (points, )
823	kwargs = {'vertex_coordinates': coordinates,
824	'triangles': triangles,
825	'mesh': None,
826	'point_attributes': values,
827	'data_origin': data_georef.get_origin(),
828	'mesh_origin': mesh_georef.get_origin(),
829	'alpha': alpha,
830	'verbose': verbose}
831
832	vertex_attributes = apply(fit_to_mesh,
833	args, kwargs)
834
835	# Call underlying method using array values
836	self.set_values_from_array(vertex_attributes,
837	location, indices, verbose)
838
839
840
841	def set_values_from_points(self, points, values, alpha,
842	location, indices,
843	data_georef = None,
844	verbose = False,
845	use_cache = False):
846	"""
847	Set quantity values from arbitray data points using
848	fit_interpolate.fit
849	"""
850
851	raise Exception, 'set_values_from_points is obsolete, use geospatial data object instead'
852
853
854	def set_values_from_file(self, filename, attribute_name, alpha,
855	location, indices,
856	verbose = False,
857	use_cache = False,
858	max_read_lines=None):
859	"""Set quantity based on arbitrary points in a points file
860	using attribute_name selects name of attribute
861	present in file.
862	If attribute_name is not specified, use first available attribute
863	as defined in geospatial_data.
864	"""
865
866	from types import StringType
867	msg = 'Filename must be a text string'
868	assert type(filename) == StringType, msg
869
870
871	if location != 'vertices':
872	msg = 'set_values_from_file is only defined for '+\
873	'location=\'vertices\''
874	raise msg
875
876	if True:
877	vertex_attributes = fit_to_mesh(filename,
878	mesh=self.domain,
879	alpha=alpha,
880	attribute_name=attribute_name,
881	use_cache=use_cache,
882	verbose=verbose,
883	max_read_lines=max_read_lines)
884	else:
885
886	coordinates = self.domain.get_nodes(absolute=True)
887	triangles = self.domain.triangles #FIXME
888	vertex_attributes = fit_to_mesh(filename,
889	coordinates, triangles,
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	def get_extremum_index(self, mode=None, indices=None):
902	"""Return index for maximum or minimum value of quantity (on centroids)
903
904	Optional arguments:
905	mode is either 'max'(default) or 'min'.
906	indices is the set of element ids that the operation applies to.
907
908	Usage:
909	i = get_extreme_index()
910
911	Notes:
912	We do not seek the extremum at vertices as each vertex can
913	have multiple values - one for each triangle sharing it.
914
915	If there are multiple cells with same maximum value, the
916	first cell encountered in the triangle array is returned.
917	"""
918
919	V = self.get_values(location='centroids', indices=indices)
920
921	# Always return absolute indices
922	if mode is None or mode == 'max':
923	i = argmax(V)
924	elif mode == 'min':
925	i = argmin(V)
926
927
928	if indices is None:
929	return i
930	else:
931	return indices[i]
932
933
934	def get_maximum_index(self, indices=None):
935	"""See get extreme index for details
936	"""
937
938	return self.get_extremum_index(mode='max',
939	indices=indices)
940
941
942
943	def get_maximum_value(self, indices=None):
944	"""Return maximum value of quantity (on centroids)
945
946	Optional argument:
947	indices is the set of element ids that the operation applies to.
948
949	Usage:
950	v = get_maximum_value()
951
952	Note, we do not seek the maximum at vertices as each vertex can
953	have multiple values - one for each triangle sharing it
954	"""
955
956
957	i = self.get_maximum_index(indices)
958	V = self.get_values(location='centroids') #, indices=indices)
959
960	return V[i]
961
962
963	def get_maximum_location(self, indices=None):
964	"""Return location of maximum 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	x, y = get_maximum_location()
971
972
973	Notes:
974	We do not seek the maximum at vertices as each vertex can
975	have multiple values - one for each triangle sharing it.
976
977	If there are multiple cells with same maximum value, the
978	first cell encountered in the triangle array is returned.
979	"""
980
981	i = self.get_maximum_index(indices)
982	x, y = self.domain.get_centroid_coordinates()[i]
983
984	return x, y
985
986
987	def get_minimum_index(self, indices=None):
988	"""See get extreme index for details
989	"""
990
991	return self.get_extremum_index(mode='min',
992	indices=indices)
993
994
995	def get_minimum_value(self, indices=None):
996	"""Return minimum value of quantity (on centroids)
997
998	Optional argument:
999	indices is the set of element ids that the operation applies to.
1000
1001	Usage:
1002	v = get_minimum_value()
1003
1004	See get_maximum_value for more details.
1005	"""
1006
1007
1008	i = self.get_minimum_index(indices)
1009	V = self.get_values(location='centroids')
1010
1011	return V[i]
1012
1013
1014	def get_minimum_location(self, indices=None):
1015	"""Return location of minimum value of quantity (on centroids)
1016
1017	Optional argument:
1018	indices is the set of element ids that the operation applies to.
1019
1020	Usage:
1021	x, y = get_minimum_location()
1022
1023
1024	Notes:
1025	We do not seek the maximum at vertices as each vertex can
1026	have multiple values - one for each triangle sharing it.
1027
1028	If there are multiple cells with same maximum value, the
1029	first cell encountered in the triangle array is returned.
1030	"""
1031
1032	i = self.get_minimum_index(indices)
1033	x, y = self.domain.get_centroid_coordinates()[i]
1034
1035	return x, y
1036
1037
1038
1039
1040	def get_interpolated_values(self, interpolation_points):
1041
1042	# Interpolation object based on internal (discontinuous triangles)
1043	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
1044	smooth=False)
1045	# FIXME: 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, I._point_coordinates):
1055	can_reuse = True
1056
1057
1058	if can_reuse is True:
1059	# Use absence of points to indicate reuse in I.interpolate
1060	result = I.interpolate(vertex_values)
1061	else:
1062	from anuga.fit_interpolate.interpolate import Interpolate
1063
1064	# Create interpolation object with matrix
1065	I = Interpolate(vertex_coordinates, triangles)
1066	self.interpolation_object = I
1067
1068	# Call interpolate with points the first time
1069	interpolation_points = ensure_numeric(interpolation_points, Float)
1070	result = I.interpolate(vertex_values, interpolation_points)
1071
1072	return result
1073
1074
1075	def get_values(self, interpolation_points=None,
1076	location='vertices',
1077	indices = None):
1078	"""get values for quantity
1079
1080	return X, Compatible list, Numeric array (see below)
1081	interpolation_points: List of x, y coordinates where value is
1082	sought (using interpolation). If points are given, values of
1083	location and indices are ignored
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	from Numeric import take
1110
1111	# FIXME (Ole): I reckon we should have the option of passing a
1112	# polygon into get_values. The question becomes how
1113	# resulting values should be ordered.
1114
1115	if interpolation_points is not None:
1116	return self.get_interpolated_values(interpolation_points)
1117
1118
1119
1120	if location not in ['vertices', 'centroids', 'edges',
1121	'unique vertices']:
1122	msg = 'Invalid location: %s' %location
1123	raise msg
1124
1125	import types, Numeric
1126	assert type(indices) in [types.ListType, types.NoneType,
1127	Numeric.ArrayType],\
1128	'Indices must be a list or None'
1129
1130	if location == 'centroids':
1131	if (indices == None):
1132	indices = range(len(self))
1133	return take(self.centroid_values,indices)
1134	elif location == 'edges':
1135	if (indices == None):
1136	indices = range(len(self))
1137	return take(self.edge_values,indices)
1138	elif location == 'unique vertices':
1139	if (indices == None):
1140	indices=range(self.domain.number_of_nodes)
1141	vert_values = []
1142	#Go through list of unique vertices
1143	for unique_vert_id in indices:
1144	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1145
1146	#In case there are unused points
1147	if len(triangles) == 0:
1148	msg = 'Unique vertex not associated with triangles'
1149	raise msg
1150
1151	# Go through all triangle, vertex pairs
1152	# Average the values
1153
1154	# FIXME (Ole): Should we merge this with get_vertex_values
1155	sum = 0
1156	for triangle_id, vertex_id in triangles:
1157	sum += self.vertex_values[triangle_id, vertex_id]
1158	vert_values.append(sum/len(triangles))
1159	return Numeric.array(vert_values)
1160	else:
1161	if (indices is None):
1162	indices = range(len(self))
1163	return take(self.vertex_values, indices)
1164
1165
1166
1167	def set_vertex_values(self, A, indices = None):
1168	"""Set vertex values for all unique vertices based on input array A
1169	which has one entry per unique vertex, i.e.
1170	one value for each row in array self.domain.nodes.
1171
1172	indices is the list of vertex_id's that will be set.
1173
1174	This function is used by set_values_from_array
1175	"""
1176
1177	from Numeric import array, Float
1178
1179	#Assert that A can be converted to a Numeric array of appropriate dim
1180	A = array(A, Float)
1181
1182	#print 'SHAPE A', A.shape
1183	assert len(A.shape) == 1
1184
1185	if indices is None:
1186	assert A.shape[0] == self.domain.get_nodes().shape[0]
1187	vertex_list = range(A.shape[0])
1188	else:
1189	assert A.shape[0] == len(indices)
1190	vertex_list = indices
1191
1192	#Go through list of unique vertices
1193
1194	for i_index, unique_vert_id in enumerate(vertex_list):
1195
1196
1197	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1198
1199	#In case there are unused points
1200	if len(triangles) == 0: continue
1201
1202	#Go through all triangle, vertex pairs
1203	#touching vertex unique_vert_id and set corresponding vertex value
1204	for triangle_id, vertex_id in triangles:
1205	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1206
1207	#Intialise centroid and edge_values
1208	self.interpolate()
1209
1210
1211	def smooth_vertex_values(self):
1212	""" Smooths vertex values.
1213	"""
1214
1215	A,V = self.get_vertex_values(xy=False, smooth = True)
1216	self.set_vertex_values(A)
1217
1218
1219	# Methods for outputting model results
1220	def get_vertex_values(self,
1221	xy=True,
1222	smooth=None,
1223	precision=None):
1224	"""Return vertex values like an OBJ format i.e. one value per node.
1225
1226	The vertex values are returned as one sequence in the 1D float array A.
1227	If requested the coordinates will be returned in 1D arrays X and Y.
1228
1229	The connectivity is represented as an integer array, V, of dimension
1230	Mx3, where M is the number of triangles. Each row has three indices
1231	defining the triangle and they correspond to elements in the arrays
1232	X, Y and A.
1233
1234	if smooth is True, vertex values corresponding to one common
1235	coordinate set will be smoothed by taking the average of vertex values for each node.
1236	In this case vertex coordinates will be
1237	de-duplicated corresponding to the original nodes as obtained from
1238	the method general_mesh.get_nodes()
1239
1240	If no smoothings is required, vertex coordinates and values will
1241	be aggregated as a concatenation of values at
1242	vertices 0, vertices 1 and vertices 2. This corresponds to
1243	the node coordinates obtained from the method
1244	general_mesh.get_vertex_coordinates()
1245
1246
1247	Calling convention
1248	if xy is True:
1249	X,Y,A,V = get_vertex_values
1250	else:
1251	A,V = get_vertex_values
1252
1253	"""
1254
1255	from Numeric import concatenate, zeros, Float, Int, array, reshape
1256
1257
1258	if smooth is None:
1259	# Take default from domain
1260	try:
1261	smooth = self.domain.smooth
1262	except:
1263	smooth = False
1264
1265	if precision is None:
1266	precision = Float
1267
1268
1269	if smooth is True:
1270	# Ensure continuous vertex values by averaging
1271	# values at each node
1272
1273	V = self.domain.get_triangles()
1274	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1275	A = zeros(N, Float)
1276	points = self.domain.get_nodes()
1277
1278	if 1:
1279	# Fast C version
1280	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1281	ensure_numeric(self.domain.number_of_triangles_per_node),
1282	ensure_numeric(self.vertex_values),
1283	A)
1284	A = A.astype(precision)
1285	else:
1286
1287	# Slow Python version
1288
1289	current_node = 0
1290	k = 0 # Track triangles touching on node
1291	total = 0.0
1292	for index in self.domain.vertex_value_indices:
1293	if current_node == N:
1294	msg = 'Current node exceeding number of nodes (%d) ' %(N)
1295	raise msg
1296
1297
1298
1299	k += 1
1300
1301	volume_id = index / 3
1302	vertex_id = index % 3
1303
1304	#assert V[volume_id, vertex_id] == current_node
1305
1306	v = self.vertex_values[volume_id, vertex_id]
1307	total += v
1308
1309	#print 'current_node=%d, index=%d, k=%d, total=%f' %(current_node, index, k, total)
1310	if self.domain.number_of_triangles_per_node[current_node] == k:
1311	A[current_node] = total/k
1312
1313
1314	# Move on to next node
1315	total = 0.0
1316	k = 0
1317	current_node += 1
1318
1319
1320
1321	else:
1322	# Allow discontinuous vertex values
1323	V = self.domain.get_disconnected_triangles()
1324	points = self.domain.get_vertex_coordinates()
1325	A = self.vertex_values.flat.astype(precision)
1326
1327
1328	# Return
1329	if xy is True:
1330	X = points[:,0].astype(precision)
1331	Y = points[:,1].astype(precision)
1332
1333	return X, Y, A, V
1334	else:
1335	return A, V
1336
1337
1338
1339	def extrapolate_first_order(self):
1340	"""Extrapolate conserved quantities from centroid to
1341	vertices and edges for each volume using
1342	first order scheme.
1343	"""
1344
1345	qc = self.centroid_values
1346	qv = self.vertex_values
1347	qe = self.edge_values
1348
1349	for i in range(3):
1350	qv[:,i] = qc
1351	qe[:,i] = qc
1352
1353	self.x_gradient *= 0.0
1354	self.y_gradient *= 0.0
1355
1356
1357	def get_integral(self):
1358	"""Compute the integral of quantity across entire domain
1359	"""
1360	integral = 0
1361	for k in range(len(self.domain)):
1362	area = self.domain.areas[k]
1363	qc = self.centroid_values[k]
1364	integral += qc*area
1365
1366	return integral
1367
1368	def get_gradients(self):
1369	"""Provide gradients. Use compute_gradients first
1370	"""
1371
1372	return self.x_gradient, self.y_gradient
1373
1374
1375	def update(self, timestep):
1376	#Call correct module function
1377	#(either from this module or C-extension)
1378	return update(self, timestep)
1379
1380	def compute_gradients(self):
1381	#Call correct module function
1382	#(either from this module or C-extension)
1383	return compute_gradients(self)
1384
1385	def limit(self):
1386	#Call correct module depending on whether
1387	#basing limit calculations on edges or vertices
1388	limit_old(self)
1389
1390	def limit_vertices_by_all_neighbours(self):
1391	#Call correct module function
1392	#(either from this module or C-extension)
1393	limit_vertices_by_all_neighbours(self)
1394
1395	def limit_edges_by_all_neighbours(self):
1396	#Call correct module function
1397	#(either from this module or C-extension)
1398	limit_edges_by_all_neighbours(self)
1399
1400	def limit_edges_by_neighbour(self):
1401	#Call correct module function
1402	#(either from this module or C-extension)
1403	limit_edges_by_neighbour(self)
1404
1405	def extrapolate_second_order(self):
1406	#Call correct module function
1407	#(either from this module or C-extension)
1408	compute_gradients(self)
1409	extrapolate_from_gradient(self)
1410
1411	def extrapolate_second_order_and_limit(self):
1412	#Call correct module function
1413	#(either from this module or C-extension)
1414	extrapolate_second_order_and_limit(self)
1415
1416	def bound_vertices_below_by_constant(self, bound):
1417	#Call correct module function
1418	#(either from this module or C-extension)
1419	bound_vertices_below_by_constant(self, bound)
1420
1421	def bound_vertices_below_by_quantity(self, quantity):
1422	#Call correct module function
1423	#(either from this module or C-extension)
1424
1425	#check consistency
1426	assert self.domain == quantity.domain
1427	bound_vertices_below_by_quantity(self, quantity)
1428
1429	def backup_centroid_values(self):
1430	#Call correct module function
1431	#(either from this module or C-extension)
1432	backup_centroid_values(self)
1433
1434	def saxpy_centroid_values(self,a,b):
1435	#Call correct module function
1436	#(either from this module or C-extension)
1437	saxpy_centroid_values(self,a,b)
1438
1439	#Conserved_quantity = Quantity
1440
1441	class Conserved_quantity(Quantity):
1442	"""Class conserved quantity being removed, use Quantity
1443
1444	"""
1445
1446	def __init__(self, domain, vertex_values=None):
1447	#Quantity.__init__(self, domain, vertex_values)
1448
1449	msg = 'ERROR: Use Quantity instead of Conserved_quantity'
1450
1451	raise Exception, msg
1452
1453
1454
1455	from anuga.utilities import compile
1456	if compile.can_use_C_extension('quantity_ext.c'):
1457	# Underlying C implementations can be accessed
1458
1459	from quantity_ext import \
1460	average_vertex_values,\
1461	backup_centroid_values,\
1462	saxpy_centroid_values,\
1463	compute_gradients,\
1464	limit_old,\
1465	limit_vertices_by_all_neighbours,\
1466	limit_edges_by_all_neighbours,\
1467	limit_edges_by_neighbour,\
1468	limit_gradient_by_neighbour,\
1469	extrapolate_from_gradient,\
1470	extrapolate_second_order_and_limit,\
1471	bound_vertices_below_by_constant,\
1472	bound_vertices_below_by_quantity,\
1473	interpolate_from_vertices_to_edges,\
1474	interpolate_from_edges_to_vertices,\
1475	update
1476	else:
1477	msg = 'C implementations could not be accessed by %s.\n ' %__file__
1478	msg += 'Make sure compile_all.py has been run as described in '
1479	msg += 'the ANUGA installation guide.'
1480	raise Exception, msg
1481
1482
1483

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

Download in other formats: