Context Navigation

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

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

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

Download in other formats: