Context Navigation

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

Last change on this file since 3945 was 3945, checked in by ole, 17 years ago
One large step towards major cleanup. This has mainly to do with the way vertex coordinates are handled internally.
File size: 49.1 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, is_scalar
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
561	# FIXME (Ole): Now we can compute the arrays once and for all
562	# for both centroids and vertices and the use set_values_from_array
563
564	if location == 'centroids':
565	V = take(self.domain.centroid_coordinates, indices)
566	if is_subset:
567	self.set_values(f(V[:,0], V[:,1]),
568	location = location,
569	indices = indices)
570	else:
571	self.set_values(f(V[:,0], V[:,1]), location = location)
572	elif location == 'vertices':
573	V = self.domain.get_vertex_coordinates()
574
575	x = V[:,0]; y = V[:,1];
576	values = f(x, y)
577
578	if is_scalar(values):
579	# Function returned a constant value
580	self.set_values_from_constant(values,
581	location, indices, verbose)
582	return
583
584
585	if is_subset:
586	#Brute force
587	for i in indices:
588	for j in range(3):
589	self.vertex_values[i,j] = values[3*i+j]
590	else:
591	for j in range(3):
592	self.vertex_values[:,j] = values[j::3]
593	else:
594	raise 'Not implemented: %s' %location
595
596
597
598	def set_values_from_geospatial_data(self, geospatial_data, alpha,
599	location, indices,
600	verbose = False,
601	use_cache = False):
602
603	#FIXME: Use this function for the time being. Later move code in here
604
605	points = geospatial_data.get_data_points(absolute = False)
606	values = geospatial_data.get_attributes()
607	data_georef = geospatial_data.get_geo_reference()
608
609
610
611	self.set_values_from_points(points, values, alpha,
612	location, indices,
613	data_georef = data_georef,
614	verbose = verbose,
615	use_cache = use_cache)
616
617
618
619	def set_values_from_points(self, points, values, alpha,
620	location, indices,
621	data_georef = None,
622	verbose = False,
623	use_cache = False):
624	"""
625	Set quantity values from arbitray data points using
626	fit_interpolate.fit
627	"""
628
629
630	from anuga.fit_interpolate.fit import fit_to_mesh
631	from anuga.coordinate_transforms.geo_reference import Geo_reference
632
633
634	points = ensure_numeric(points, Float)
635	values = ensure_numeric(values, Float)
636
637	if location != 'vertices':
638	msg = 'set_values_from_points is only defined for '+\
639	'location=\'vertices\''
640	raise ms
641
642	coordinates = self.domain.get_nodes()
643	triangles = self.domain.triangles #FIXME
644
645
646	#Take care of georeferencing
647	if data_georef is None:
648	data_georef = Geo_reference()
649
650
651	mesh_georef = self.domain.geo_reference
652
653	#print mesh_georef
654	#print data_georef
655	#print points
656
657
658	#Call fit_interpolate.fit function
659	args = (coordinates, triangles, points, values)
660	kwargs = {'data_origin': data_georef.get_origin(),
661	'mesh_origin': mesh_georef.get_origin(),
662	'alpha': alpha,
663	'verbose': verbose}
664
665	#print kwargs
666
667	if use_cache is True:
668	try:
669	from caching import cache
670	except:
671	msg = 'Caching was requested, but caching module'+\
672	'could not be imported'
673	raise msg
674
675	vertex_attributes = cache(fit_to_mesh,
676	args, kwargs,
677	verbose=verbose,
678	compression=False)
679	else:
680
681	vertex_attributes = apply(fit_to_mesh,
682	args, kwargs)
683
684	#Call underlying method using array values
685	self.set_values_from_array(vertex_attributes,
686	location, indices, verbose)
687
688
689
690	def set_values_from_file(self, filename, attribute_name, alpha,
691	location, indices,
692	verbose = False,
693	use_cache = False):
694	"""Set quantity based on arbitrary points in .pts file
695	using attribute_name selects name of attribute
696	present in file.
697	If not specified try to use whatever is available in file.
698	"""
699
700	from load_mesh.loadASCII import import_points_file
701	from anuga.geospatial_data.geospatial_data import\
702	points_dictionary2geospatial_data
703
704	from types import StringType
705	msg = 'Filename must be a text string'
706	assert type(filename) == StringType, msg
707
708
709	# Read from (NetCDF) file
710	# FIXME (Ole): This function should really return a
711	# Geospatial_data object.
712	points_dict = import_points_file(filename)
713	points = points_dict['pointlist']
714	attributes = points_dict['attributelist']
715
716	if attribute_name is None:
717	names = attributes.keys()
718	attribute_name = names[0]
719
720	msg = 'Attribute_name must be a text string'
721	assert type(attribute_name) == StringType, msg
722
723
724	if verbose:
725	print 'Using attribute %s from file %s' %(attribute_name, filename)
726	print 'Available attributes: %s' %(names)
727
728	#try:
729	# z = attributes[attribute_name]
730	#except:
731	# msg = 'Could not extract attribute %s from file %s'\
732	# %(attribute_name, filename)
733	# raise msg
734
735
736	#Take care of georeferencing
737	if points_dict.has_key('geo_reference') and \
738	points_dict['geo_reference'] is not None:
739	data_georef = points_dict['geo_reference']
740	else:
741	data_georef = None
742
743
744
745	#Call underlying method for geospatial data
746	geospatial_data = points_dictionary2geospatial_data(points_dict)
747	geospatial_data.set_default_attribute_name(attribute_name)
748
749	self.set_values_from_geospatial_data(geospatial_data,
750	alpha,
751	location, indices,
752	verbose = verbose,
753	use_cache = use_cache)
754
755
756	def get_maximum_index(self, indices=None):
757	"""Return index for maximum value of quantity (on centroids)
758
759	Optional argument:
760	indices is the set of element ids that the operation applies to.
761
762	Usage:
763	i = get_maximum_index()
764
765	Notes:
766	We do not seek the maximum at vertices as each vertex can
767	have multiple values - one for each triangle sharing it.
768
769	If there are multiple cells with same maximum value, the
770	first cell encountered in the triangle array is returned.
771	"""
772
773	V = self.get_values(location='centroids', indices=indices)
774
775	# Always return absolute indices
776	i = argmax(V)
777
778	if indices is None:
779	return i
780	else:
781	return indices[i]
782
783
784	def get_maximum_value(self, indices=None):
785	"""Return maximum value of quantity (on centroids)
786
787	Optional argument:
788	indices is the set of element ids that the operation applies to.
789
790	Usage:
791	v = get_maximum_value()
792
793	Note, we do not seek the maximum at vertices as each vertex can
794	have multiple values - one for each triangle sharing it
795	"""
796
797
798	i = self.get_maximum_index(indices)
799	V = self.get_values(location='centroids') #, indices=indices)
800
801	return V[i]
802
803
804	def get_maximum_location(self, indices=None):
805	"""Return location of maximum value of quantity (on centroids)
806
807	Optional argument:
808	indices is the set of element ids that the operation applies to.
809
810	Usage:
811	x, y = get_maximum_location()
812
813
814	Notes:
815	We do not seek the maximum at vertices as each vertex can
816	have multiple values - one for each triangle sharing it.
817
818	If there are multiple cells with same maximum value, the
819	first cell encountered in the triangle array is returned.
820	"""
821
822	i = self.get_maximum_index(indices)
823	x, y = self.domain.get_centroid_coordinates()[i]
824
825	return x, y
826
827
828
829
830	def get_interpolated_values(self, interpolation_points):
831
832	# Interpolation object based on internal (discontinuous triangles)
833	x, y, vertex_values, triangles = self.get_vertex_values(xy=True,
834	smooth=False)
835	# FIXME: This concat should roll into get_vertex_values
836	vertex_coordinates = concatenate((x[:, NewAxis], y[:, NewAxis]),
837	axis=1)
838
839	can_reuse = False
840	if hasattr(self, 'interpolation_object'):
841	# Reuse to save time
842	I = self.interpolation_object
843
844	if allclose(interpolation_points, I._point_coordinates):
845	can_reuse = True
846
847
848	if can_reuse is True:
849	# Use absence of points to indicate reuse in I.interpolate
850	result = I.interpolate(vertex_values)
851	else:
852	from anuga.fit_interpolate.interpolate import Interpolate
853
854	# Create interpolation object with matrix
855	I = Interpolate(vertex_coordinates, triangles)
856	self.interpolation_object = I
857
858	# Call interpolate with points the first time
859	interpolation_points = ensure_numeric(interpolation_points, Float)
860	result = I.interpolate(vertex_values, interpolation_points)
861
862	return result
863
864
865	def get_values(self, interpolation_points=None,
866	location='vertices',
867	indices = None):
868	"""get values for quantity
869
870	return X, Compatible list, Numeric array (see below)
871	interpolation_points: List of x, y coordinates where value is
872	sought (using interpolation). If points are given, values of
873	location and indices are ignored
874
875	location: Where values are to be stored.
876	Permissible options are: vertices, edges, centroid
877	and unique vertices. Default is 'vertices'
878
879
880	The returned values with be a list the length of indices
881	(N if indices = None).
882
883	In case of location == 'centroids' the dimension of returned
884	values will be a list or a Numerical array of length N, N being
885	the number of elements.
886
887	In case of location == 'vertices' or 'edges' the dimension of
888	returned values will be of dimension Nx3
889
890	In case of location == 'unique vertices' the average value at
891	each vertex will be returned and the dimension of returned values
892	will be a 1d array of length "number of vertices"
893
894	Indices is the set of element ids that the operation applies to.
895
896	The values will be stored in elements following their
897	internal ordering.
898
899	"""
900	from Numeric import take
901
902	if interpolation_points is not None:
903	return self.get_interpolated_values(interpolation_points)
904
905
906
907	if location not in ['vertices', 'centroids', 'edges',
908	'unique vertices']:
909	msg = 'Invalid location: %s' %location
910	raise msg
911
912	import types, Numeric
913	assert type(indices) in [types.ListType, types.NoneType,
914	Numeric.ArrayType],\
915	'Indices must be a list or None'
916
917	if location == 'centroids':
918	if (indices == None):
919	indices = range(len(self))
920	return take(self.centroid_values,indices)
921	elif location == 'edges':
922	if (indices == None):
923	indices = range(len(self))
924	return take(self.edge_values,indices)
925	elif location == 'unique vertices':
926	if (indices == None):
927	indices=range(self.domain.number_of_nodes)
928	vert_values = []
929	#Go through list of unique vertices
930	for unique_vert_id in indices:
931	triangles = self.domain.vertexlist[unique_vert_id]
932
933	#In case there are unused points
934	if triangles is None:
935	msg = 'Unique vertex not associated with triangles'
936	raise msg
937
938	# Go through all triangle, vertex pairs
939	# Average the values
940
941	# FIXME (Ole): Should we merge this with get_vertex_values
942	# and use the concept of a reduction operator?
943	sum = 0
944	for triangle_id, vertex_id in triangles:
945	sum += self.vertex_values[triangle_id, vertex_id]
946	vert_values.append(sum/len(triangles))
947	return Numeric.array(vert_values)
948	else:
949	if (indices == None):
950	indices = range(len(self))
951	return take(self.vertex_values,indices)
952
953
954
955	def set_vertex_values(self, A, indices = None):
956	"""Set vertex values for all unique vertices based on input array A
957	which has one entry per unique vertex, i.e.
958	one value for each row in array self.domain.coordinates or
959	one value for each row in vertexlist.
960
961	indices is the list of vertex_id's that will be set.
962
963	This function is used by set_values_from_array
964	"""
965
966	from Numeric import array, Float
967
968	#Assert that A can be converted to a Numeric array of appropriate dim
969	A = array(A, Float)
970
971	#print 'SHAPE A', A.shape
972	assert len(A.shape) == 1
973
974	if indices == None:
975	assert A.shape[0] == self.domain.get_nodes().shape[0]
976	vertex_list = range(A.shape[0])
977	else:
978	assert A.shape[0] == len(indices)
979	vertex_list = indices
980
981	#Go through list of unique vertices
982	for i_index, unique_vert_id in enumerate(vertex_list):
983	triangles = self.domain.vertexlist[unique_vert_id]
984
985	if triangles is None: continue #In case there are unused points
986
987	#Go through all triangle, vertex pairs
988	#touching vertex unique_vert_id and set corresponding vertex value
989	for triangle_id, vertex_id in triangles:
990	self.vertex_values[triangle_id, vertex_id] = A[i_index]
991
992	#Intialise centroid and edge_values
993	self.interpolate()
994
995
996	def smooth_vertex_values(self, value_array='field_values',
997	precision = None):
998	""" Smooths field_values or conserved_quantities data.
999	TODO: be able to smooth individual fields
1000	NOTE: This function does not have a test.
1001	FIXME: NOT DONE - do we need it?
1002	FIXME: this function isn't called by anything.
1003	Maybe it should be removed..-DSG
1004	"""
1005
1006	from Numeric import concatenate, zeros, Float, Int, array, reshape
1007
1008
1009	A,V = self.get_vertex_values(xy=False,
1010	value_array=value_array,
1011	smooth = True,
1012	precision = precision)
1013
1014	#Set some field values
1015	for volume in self:
1016	for i,v in enumerate(volume.vertices):
1017	if value_array == 'field_values':
1018	volume.set_field_values('vertex', i, A[v,:])
1019	elif value_array == 'conserved_quantities':
1020	volume.set_conserved_quantities('vertex', i, A[v,:])
1021
1022	if value_array == 'field_values':
1023	self.precompute()
1024	elif value_array == 'conserved_quantities':
1025	Volume.interpolate_conserved_quantities()
1026
1027
1028	#Method for outputting model results
1029	#FIXME: Split up into geometric and numeric stuff.
1030	#FIXME: Geometric (X,Y,V) should live in mesh.py
1031	#FIXME: STill remember to move XY to mesh
1032	def get_vertex_values(self,
1033	xy=True,
1034	smooth = None,
1035	precision = None,
1036	reduction = None):
1037	"""Return vertex values like an OBJ format
1038
1039	The vertex values are returned as one sequence in the 1D float array A.
1040	If requested the coordinates will be returned in 1D arrays X and Y.
1041
1042	The connectivity is represented as an integer array, V, of dimension
1043	M x 3, where M is the number of volumes. Each row has three indices
1044	into the X, Y, A arrays defining the triangle.
1045
1046	if smooth is True, vertex values corresponding to one common
1047	coordinate set will be smoothed according to the given
1048	reduction operator. In this case vertex coordinates will be
1049	de-duplicated.
1050
1051	If no smoothings is required, vertex coordinates and values will
1052	be aggregated as a concatenation of values at
1053	vertices 0, vertices 1 and vertices 2
1054
1055
1056	Calling convention
1057	if xy is True:
1058	X,Y,A,V = get_vertex_values
1059	else:
1060	A,V = get_vertex_values
1061
1062	"""
1063
1064	from Numeric import concatenate, zeros, Float, Int, array, reshape
1065
1066
1067	if smooth is None:
1068	smooth = self.domain.smooth
1069
1070	if precision is None:
1071	precision = Float
1072
1073	#Create connectivity
1074
1075	if smooth == True:
1076
1077	if reduction is None:
1078	reduction = self.domain.reduction
1079
1080	V = self.domain.get_vertices()
1081	N = len(self.domain.vertexlist)
1082	A = zeros(N, precision)
1083
1084	#Smoothing loop
1085	for k in range(N):
1086	L = self.domain.vertexlist[k]
1087
1088	#Go through all triangle, vertex pairs
1089	#contributing to vertex k and register vertex value
1090
1091	if L is None: continue #In case there are unused points
1092
1093	contributions = []
1094	for volume_id, vertex_id in L:
1095	v = self.vertex_values[volume_id, vertex_id]
1096	contributions.append(v)
1097
1098	A[k] = reduction(contributions)
1099
1100
1101	if xy is True:
1102	X = self.domain.get_nodes()[:,0].astype(precision)
1103	Y = self.domain.get_nodes()[:,1].astype(precision)
1104
1105	return X, Y, A, V
1106	else:
1107	return A, V
1108	else:
1109	#Don't smooth
1110	#obj machinery moved to general_mesh
1111
1112	# Create a V like [[0 1 2], [3 4 5]....[3m-2 3m-1 3*m]]
1113	# These vert_id's will relate to the verts created below
1114	#m = len(self.domain) #Number of volumes
1115	#M = 3*m #Total number of unique vertices
1116	#V = reshape(array(range(M)).astype(Int), (m,3))
1117
1118	V = self.domain.get_triangles(obj=True)
1119	#FIXME use get_vertices, when ready
1120
1121	A = self.vertex_values.flat
1122
1123	#Do vertex coordinates
1124	if xy is True:
1125	C = self.domain.get_vertex_coordinates()
1126
1127	X = C[:,0:6:2].copy()
1128	Y = C[:,1:6:2].copy()
1129
1130	return X.flat, Y.flat, A, V
1131	else:
1132	return A, V
1133
1134
1135	def extrapolate_first_order(self):
1136	"""Extrapolate conserved quantities from centroid to
1137	vertices for each volume using
1138	first order scheme.
1139	"""
1140
1141	qc = self.centroid_values
1142	qv = self.vertex_values
1143
1144	for i in range(3):
1145	qv[:,i] = qc
1146
1147
1148	def get_integral(self):
1149	"""Compute the integral of quantity across entire domain
1150	"""
1151	integral = 0
1152	for k in range(len(self.domain)):
1153	area = self.domain.areas[k]
1154	qc = self.centroid_values[k]
1155	integral += qc*area
1156
1157	return integral
1158
1159
1160
1161
1162	class Conserved_quantity(Quantity):
1163	"""Class conserved quantity adds to Quantity:
1164
1165	boundary values, storage and method for updating, and
1166	methods for (second order) extrapolation from centroid to vertices inluding
1167	gradients and limiters
1168	"""
1169
1170	def __init__(self, domain, vertex_values=None):
1171	Quantity.__init__(self, domain, vertex_values)
1172
1173	from Numeric import zeros, Float
1174
1175	#Allocate space for boundary values
1176	L = len(domain.boundary)
1177	self.boundary_values = zeros(L, Float)
1178
1179	#Allocate space for updates of conserved quantities by
1180	#flux calculations and forcing functions
1181
1182	N = len(domain) # number_of_triangles
1183	self.explicit_update = zeros(N, Float )
1184	self.semi_implicit_update = zeros(N, Float )
1185
1186
1187	def update(self, timestep):
1188	#Call correct module function
1189	#(either from this module or C-extension)
1190	return update(self, timestep)
1191
1192
1193	def compute_gradients(self):
1194	#Call correct module function
1195	#(either from this module or C-extension)
1196	return compute_gradients(self)
1197
1198
1199	def limit(self):
1200	#Call correct module function
1201	#(either from this module or C-extension)
1202	limit(self)
1203
1204
1205	def extrapolate_second_order(self):
1206	#Call correct module function
1207	#(either from this module or C-extension)
1208	extrapolate_second_order(self)
1209
1210
1211	def update(quantity, timestep):
1212	"""Update centroid values based on values stored in
1213	explicit_update and semi_implicit_update as well as given timestep
1214
1215	Function implementing forcing terms must take on argument
1216	which is the domain and they must update either explicit
1217	or implicit updates, e,g,:
1218
1219	def gravity(domain):
1220	....
1221	domain.quantities['xmomentum'].explicit_update = ...
1222	domain.quantities['ymomentum'].explicit_update = ...
1223
1224
1225
1226	Explicit terms must have the form
1227
1228	G(q, t)
1229
1230	and explicit scheme is
1231
1232	q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
1233
1234
1235	Semi implicit forcing terms are assumed to have the form
1236
1237	G(q, t) = H(q, t) q
1238
1239	and the semi implicit scheme will then be
1240
1241	q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
1242
1243
1244	"""
1245
1246	from Numeric import sum, equal, ones, exp, Float
1247
1248	N = quantity.centroid_values.shape[0]
1249
1250
1251	#Divide H by conserved quantity to obtain G (see docstring above)
1252
1253
1254	for k in range(N):
1255	x = quantity.centroid_values[k]
1256	if x == 0.0:
1257	#FIXME: Is this right
1258	quantity.semi_implicit_update[k] = 0.0
1259	else:
1260	quantity.semi_implicit_update[k] /= x
1261
1262
1263	#Semi implicit updates
1264	denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1265
1266	if sum(less(denominator, 1.0)) > 0.0:
1267	msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
1268	raise msg
1269
1270	if sum(equal(denominator, 0.0)) > 0.0:
1271	msg = 'Zero division in semi implicit update. Call Stephen :-)'
1272	raise msg
1273	else:
1274	#Update conserved_quantities from semi implicit updates
1275	quantity.centroid_values /= denominator
1276
1277	# quantity.centroid_values = exp(timestepquantity.semi_implicit_update)quantity.centroid_values
1278
1279	#Explicit updates
1280	quantity.centroid_values += timestep*quantity.explicit_update
1281
1282	def interpolate_from_vertices_to_edges(quantity):
1283	"""Compute edge values from vertex values using linear interpolation
1284	"""
1285
1286	for k in range(quantity.vertex_values.shape[0]):
1287	q0 = quantity.vertex_values[k, 0]
1288	q1 = quantity.vertex_values[k, 1]
1289	q2 = quantity.vertex_values[k, 2]
1290
1291	quantity.edge_values[k, 0] = 0.5*(q1+q2)
1292	quantity.edge_values[k, 1] = 0.5*(q0+q2)
1293	quantity.edge_values[k, 2] = 0.5*(q0+q1)
1294
1295
1296
1297	def extrapolate_second_order(quantity):
1298	"""Extrapolate conserved quantities from centroid to
1299	vertices for each volume using
1300	second order scheme.
1301	"""
1302
1303	a, b = quantity.compute_gradients()
1304
1305	X = quantity.domain.get_vertex_coordinates()
1306	qc = quantity.centroid_values
1307	qv = quantity.vertex_values
1308
1309	#Check each triangle
1310	for k in range(len(quantity.domain)):
1311	#Centroid coordinates
1312	x, y = quantity.domain.centroid_coordinates[k]
1313
1314	#vertex coordinates
1315	x0, y0, x1, y1, x2, y2 = X[k,:]
1316
1317	#Extrapolate
1318	qv[k,0] = qc[k] + a[k](x0-x) + b[k](y0-y)
1319	qv[k,1] = qc[k] + a[k](x1-x) + b[k](y1-y)
1320	qv[k,2] = qc[k] + a[k](x2-x) + b[k](y2-y)
1321
1322
1323	def compute_gradients(quantity):
1324	"""Compute gradients of triangle surfaces defined by centroids of
1325	neighbouring volumes.
1326	If one edge is on the boundary, use own centroid as neighbour centroid.
1327	If two or more are on the boundary, fall back to first order scheme.
1328	"""
1329
1330	from Numeric import zeros, Float
1331	from utilitites.numerical_tools import gradient
1332
1333	centroid_coordinates = quantity.domain.centroid_coordinates
1334	surrogate_neighbours = quantity.domain.surrogate_neighbours
1335	centroid_values = quantity.centroid_values
1336	number_of_boundaries = quantity.domain.number_of_boundaries
1337
1338	N = centroid_values.shape[0]
1339
1340	a = zeros(N, Float)
1341	b = zeros(N, Float)
1342
1343	for k in range(N):
1344	if number_of_boundaries[k] < 2:
1345	#Two or three true neighbours
1346
1347	#Get indices of neighbours (or self when used as surrogate)
1348	k0, k1, k2 = surrogate_neighbours[k,:]
1349
1350	#Get data
1351	q0 = centroid_values[k0]
1352	q1 = centroid_values[k1]
1353	q2 = centroid_values[k2]
1354
1355	x0, y0 = centroid_coordinates[k0] #V0 centroid
1356	x1, y1 = centroid_coordinates[k1] #V1 centroid
1357	x2, y2 = centroid_coordinates[k2] #V2 centroid
1358
1359	#Gradient
1360	a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1361
1362	elif number_of_boundaries[k] == 2:
1363	#One true neighbour
1364
1365	#Get index of the one neighbour
1366	for k0 in surrogate_neighbours[k,:]:
1367	if k0 != k: break
1368	assert k0 != k
1369
1370	k1 = k #self
1371
1372	#Get data
1373	q0 = centroid_values[k0]
1374	q1 = centroid_values[k1]
1375
1376	x0, y0 = centroid_coordinates[k0] #V0 centroid
1377	x1, y1 = centroid_coordinates[k1] #V1 centroid
1378
1379	#Gradient
1380	a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1381	else:
1382	#No true neighbours -
1383	#Fall back to first order scheme
1384	pass
1385
1386
1387	return a, b
1388
1389
1390
1391	def limit(quantity):
1392	"""Limit slopes for each volume to eliminate artificial variance
1393	introduced by e.g. second order extrapolator
1394
1395	This is an unsophisticated limiter as it does not take into
1396	account dependencies among quantities.
1397
1398	precondition:
1399	vertex values are estimated from gradient
1400	postcondition:
1401	vertex values are updated
1402	"""
1403
1404	from Numeric import zeros, Float
1405
1406	N = quantity.domain.number_of_nodes
1407
1408	beta_w = quantity.domain.beta_w
1409
1410	qc = quantity.centroid_values
1411	qv = quantity.vertex_values
1412
1413	#Find min and max of this and neighbour's centroid values
1414	qmax = zeros(qc.shape, Float)
1415	qmin = zeros(qc.shape, Float)
1416
1417	for k in range(N):
1418	qmax[k] = qc[k]
1419	qmin[k] = qc[k]
1420	for i in range(3):
1421	n = quantity.domain.neighbours[k,i]
1422	if n >= 0:
1423	qn = qc[n] #Neighbour's centroid value
1424
1425	qmin[k] = min(qmin[k], qn)
1426	qmax[k] = max(qmax[k], qn)
1427	qmax[k] = min(qmax[k], 2.0*qc[k])
1428	qmin[k] = max(qmin[k], 0.5*qc[k])
1429
1430
1431	#Diffences between centroids and maxima/minima
1432	dqmax = qmax - qc
1433	dqmin = qmin - qc
1434
1435	#Deltas between vertex and centroid values
1436	dq = zeros(qv.shape, Float)
1437	for i in range(3):
1438	dq[:,i] = qv[:,i] - qc
1439
1440	#Phi limiter
1441	for k in range(N):
1442
1443	#Find the gradient limiter (phi) across vertices
1444	phi = 1.0
1445	for i in range(3):
1446	r = 1.0
1447	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1448	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1449
1450	phi = min( min(r*beta_w, 1), phi )
1451
1452	#Then update using phi limiter
1453	for i in range(3):
1454	qv[k,i] = qc[k] + phi*dq[k,i]
1455
1456
1457
1458	from anuga.utilities import compile
1459	if compile.can_use_C_extension('quantity_ext.c'):
1460	#Replace python version with c implementations
1461
1462	from quantity_ext import compute_gradients, limit,\
1463	extrapolate_second_order, interpolate_from_vertices_to_edges, update

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

Download in other formats: