Context Navigation

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

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

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

Download in other formats: