Context Navigation

source: branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py @ 6689

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

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

Download in other formats: