Context Navigation

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

Last change on this file since 3928 was 3928, checked in by ole, 17 years ago
Parallel domains now store only full triangles in sww files. Still need to remove ghost nodes.
File size: 48.6 KB

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

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

Download in other formats: