Context Navigation

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

Last change on this file since 7711 was 7711, checked in by hudson, 14 years ago
Refactored geometry classes to live in their own folder.
File size: 59.0 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 types import FloatType, IntType, LongType, NoneType
18
19	from anuga.utilities.numerical_tools import ensure_numeric, is_scalar
20	from anuga.geometry.polygon import inside_polygon
21	from anuga.geospatial_data.geospatial_data import Geospatial_data
22	from anuga.fit_interpolate.fit import fit_to_mesh
23	from anuga.config import points_file_block_line_size as default_block_line_size
24	from anuga.config import epsilon
25	from anuga.caching import cache
26	import anuga.utilities.log as log
27
28	import anuga.utilities.numerical_tools as aunt
29
30	import numpy as num
31
32
33	##
34	# @brief Implement values at each triangular element.
35	class Quantity:
36
37	##
38	# @brief Construct values art each triangular element.
39	# @param domain ??
40	# @param vertex_values ??
41	def __init__(self, domain, vertex_values=None):
42	from anuga.abstract_2d_finite_volumes.domain import Domain
43
44	msg = ('First argument in Quantity.__init__() must be of class Domain '
45	'(or a subclass thereof). I got %s.' % str(domain.__class__))
46	assert isinstance(domain, Domain), msg
47
48	if vertex_values is None:
49	N = len(domain) # number_of_elements
50	self.vertex_values = num.zeros((N, 3), num.float)
51	else:
52	self.vertex_values = num.array(vertex_values, num.float)
53
54	N, V = self.vertex_values.shape
55	assert V == 3, 'Three vertex values per element must be specified'
56
57	msg = 'Number of vertex values (%d) must be consistent with' % N
58	msg += 'number of elements in specified domain (%d).' % len(domain)
59	assert N == len(domain), msg
60
61	self.domain = domain
62
63	# Allocate space for other quantities
64	self.centroid_values = num.zeros(N, num.float)
65	self.edge_values = num.zeros((N, 3), num.float)
66
67	# Allocate space for Gradient
68	self.x_gradient = num.zeros(N, num.float)
69	self.y_gradient = num.zeros(N, num.float)
70
71	# Allocate space for Limiter Phi
72	self.phi = num.zeros(N, num.float)
73
74	# Intialise centroid and edge_values
75	self.interpolate()
76
77	# Allocate space for boundary values
78	L = len(domain.boundary)
79	self.boundary_values = num.zeros(L, num.float)
80
81	# Allocate space for updates of conserved quantities by
82	# flux calculations and forcing functions
83
84	# Allocate space for update fields
85	self.explicit_update = num.zeros(N, num.float )
86	self.semi_implicit_update = num.zeros(N, num.float )
87	self.centroid_backup_values = num.zeros(N, num.float)
88
89	self.set_beta(1.0)
90
91	############################################################################
92	# Methods for operator overloading
93	############################################################################
94
95	def __len__(self):
96	return self.centroid_values.shape[0]
97
98	def __neg__(self):
99	"""Negate all values in this quantity giving meaning to the
100	expression -Q where Q is an instance of class Quantity
101	"""
102
103	Q = Quantity(self.domain)
104	Q.set_values(-self.vertex_values)
105	return Q
106
107	def __add__(self, other):
108	"""Add to self anything that could populate a quantity
109
110	E.g other can be a constant, an array, a function, another quantity
111	(except for a filename or points, attributes (for now))
112	- see set_values for details
113	"""
114
115	Q = Quantity(self.domain)
116	Q.set_values(other)
117
118	result = Quantity(self.domain)
119	result.set_values(self.vertex_values + Q.vertex_values)
120	return result
121
122	def __radd__(self, other):
123	"""Handle cases like 7+Q, where Q is an instance of class Quantity
124	"""
125
126	return self + other
127
128	def __sub__(self, other):
129	return self + -other # Invoke self.__neg__()
130
131	def __mul__(self, other):
132	"""Multiply self with anything that could populate a quantity
133
134	E.g other can be a constant, an array, a function, another quantity
135	(except for a filename or points, attributes (for now))
136	- see set_values for details
137	"""
138
139	if isinstance(other, Quantity):
140	Q = other
141	else:
142	Q = Quantity(self.domain)
143	Q.set_values(other)
144
145	result = Quantity(self.domain)
146
147	# The product of vertex_values, edge_values and centroid_values
148	# are calculated and assigned directly without using
149	# set_values (which calls interpolate). Otherwise
150	# edge and centroid values wouldn't be products from q1 and q2
151	result.vertex_values = self.vertex_values * Q.vertex_values
152	result.edge_values = self.edge_values * Q.edge_values
153	result.centroid_values = self.centroid_values * Q.centroid_values
154
155	return result
156
157	def __rmul__(self, other):
158	"""Handle cases like 3*Q, where Q is an instance of class Quantity
159	"""
160
161	return self * other
162
163	def __div__(self, other):
164	"""Divide self with anything that could populate a quantity
165
166	E.g other can be a constant, an array, a function, another quantity
167	(except for a filename or points, attributes (for now))
168	- see set_values for details
169
170	Zero division is dealt with by adding an epsilon to the divisore
171	FIXME (Ole): Replace this with native INF once we migrate to NumPy
172	"""
173
174	if isinstance(other, Quantity):
175	Q = other
176	else:
177	Q = Quantity(self.domain)
178	Q.set_values(other)
179
180	result = Quantity(self.domain)
181
182	# The quotient of vertex_values, edge_values and centroid_values
183	# are calculated and assigned directly without using
184	# set_values (which calls interpolate). Otherwise
185	# edge and centroid values wouldn't be quotient of q1 and q2
186	result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
187	result.edge_values = self.edge_values/(Q.edge_values + epsilon)
188	result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
189
190	return result
191
192	def __rdiv__(self, other):
193	"""Handle cases like 3/Q, where Q is an instance of class Quantity
194	"""
195
196	return self / other
197
198	def __pow__(self, other):
199	"""Raise quantity to (numerical) power
200
201	As with __mul__ vertex values are processed entry by entry
202	while centroid and edge values are re-interpolated.
203
204	Example using __pow__:
205	Q = (Q12 + Q22)**0.5
206	"""
207
208	if isinstance(other, Quantity):
209	Q = other
210	else:
211	Q = Quantity(self.domain)
212	Q.set_values(other)
213
214	result = Quantity(self.domain)
215
216	# The power of vertex_values, edge_values and centroid_values
217	# are calculated and assigned directly without using
218	# set_values (which calls interpolate). Otherwise
219	# edge and centroid values wouldn't be correct
220	result.vertex_values = self.vertex_values ** other
221	result.edge_values = self.edge_values ** other
222	result.centroid_values = self.centroid_values ** other
223
224	return result
225
226	############################################################################
227	# Setters/Getters
228	############################################################################
229
230	##
231	# @brief Set default beta value for limiting.
232	# @param beta ??
233	def set_beta(self, beta):
234	"""Set default beta value for limiting """
235
236	if beta < 0.0:
237	log.critical('WARNING: setting beta < 0.0')
238	if beta > 2.0:
239	log.critical('WARNING: setting beta > 2.0')
240
241	self.beta = beta
242
243	##
244	# @brief Get the current beta value.
245	# @return The current beta value.
246	def get_beta(self):
247	"""Get default beta value for limiting"""
248
249	return self.beta
250
251	##
252	# @brief Compute interpolated values at edges and centroid.
253	# @note vertex_values must be set before calling this.
254	def interpolate(self):
255	"""Compute interpolated values at edges and centroid
256	Pre-condition: vertex_values have been set
257	"""
258
259	# FIXME (Ole): Maybe this function
260	# should move to the C-interface?
261	# However, it isn't called by validate_all.py, so it
262	# may not be that important to optimise it?
263
264	N = self.vertex_values.shape[0]
265	for i in range(N):
266	v0 = self.vertex_values[i, 0]
267	v1 = self.vertex_values[i, 1]
268	v2 = self.vertex_values[i, 2]
269
270	self.centroid_values[i] = (v0 + v1 + v2)/3
271
272	self.interpolate_from_vertices_to_edges()
273
274	##
275	# @brief ??
276	def interpolate_from_vertices_to_edges(self):
277	# Call correct module function (either from this module or C-extension)
278	interpolate_from_vertices_to_edges(self)
279
280	##
281	# @brief ??
282	def interpolate_from_edges_to_vertices(self):
283	# Call correct module function (either from this module or C-extension)
284	interpolate_from_edges_to_vertices(self)
285
286	#---------------------------------------------
287	# Public interface for setting quantity values
288	#---------------------------------------------
289
290	##
291	# @brief Set values for quantity based on different sources.
292	# @param numeric A num array, list or constant value.
293	# @param quantity Another Quantity.
294	# @param function Any callable object that takes two 1d arrays.
295	# @param geospatial_data Arbitrary instance of class Geospatial_data
296	# @param filename Path to a points file.
297	# @param attribute_name If specified any array using that name will be used.
298	# @param alpha Smoothing parameter to be used with fit_interpolate.fit.
299	# @param location Where to store values (vertices, edges, centroids).
300	# @param polygon Restrict update to locations that fall inside polygon.
301	# @param indices Restrict update to locations specified by this.
302	# @param smooth If True, smooth vertex values.
303	# @param verbose True if this method is to be verbose.
304	# @param use_cache If True cache results for fit_interpolate.fit.
305	# @note Exactly one of 'numeric', 'quantity', 'function', 'filename'
306	# must be present.
307	def set_values(self, numeric=None, # List, numeric array or constant
308	quantity=None, # Another quantity
309	function=None, # Callable object: f(x,y)
310	geospatial_data=None, # Arbitrary dataset
311	filename=None,
312	attribute_name=None, # Input from file
313	alpha=None,
314	location='vertices',
315	polygon=None,
316	indices=None,
317	smooth=False,
318	verbose=False,
319	use_cache=False):
320	"""Set values for quantity based on different sources.
321
322	numeric:
323	Compatible list, numeric array (see below) or constant.
324	If callable it will treated as a function (see below)
325	If instance of another Quantity it will be treated as such.
326	If geo_spatial object it will be treated as such
327
328	quantity:
329	Another quantity (compatible quantity, e.g. obtained as a
330	linear combination of quantities)
331
332	function:
333	Any callable object that takes two 1d arrays x and y
334	each of length N and returns an array also of length N.
335	The function will be evaluated at points determined by
336	location and indices in the underlying mesh.
337
338	geospatial_data:
339	Arbitrary geo spatial dataset in the form of the class
340	Geospatial_data. Mesh points are populated using
341	fit_interpolate.fit fitting
342
343	filename:
344	Name of a points file containing data points and attributes for
345	use with fit_interpolate.fit.
346
347	attribute_name:
348	If specified, any array matching that name
349	will be used. from file or geospatial_data.
350	Otherwise a default will be used.
351
352	alpha:
353	Smoothing parameter to be used with fit_interpolate.fit.
354	See module fit_interpolate.fit for further details about alpha.
355	Alpha will only be used with points, values or filename.
356	Otherwise it will be ignored.
357
358
359	location: Where values are to be stored.
360	Permissible options are: vertices, edges, centroids
361	Default is 'vertices'
362
363	In case of location == 'centroids' the dimension values must
364	be a list of a numerical array of length N,
365	N being the number of elements.
366	Otherwise it must be of dimension Nx3
367
368
369	The values will be stored in elements following their
370	internal ordering.
371
372	If location is 'unique vertices' indices refers the set
373	of node ids that the operation applies to.
374	If location is not 'unique vertices' indices refers the
375	set of triangle ids that the operation applies to.
376
377
378	If selected location is vertices, values for
379	centroid and edges will be assigned interpolated
380	values. In any other case, only values for the
381	specified locations will be assigned and the others
382	will be left undefined.
383
384
385	polygon: Restrict update of quantity to locations that fall
386	inside polygon. Polygon works by selecting indices
387	and calling set_values recursively.
388	Polygon mode has only been implemented for
389	constant values so far.
390
391	indices: Restrict update of quantity to locations that are
392	identified by indices (e.g. node ids if location
393	is 'unique vertices' or triangle ids otherwise).
394
395	verbose: True means that output to stdout is generated
396
397	use_cache: True means that caching of intermediate results is
398	attempted for fit_interpolate.fit.
399
400
401
402
403	Exactly one of the arguments
404	numeric, quantity, function, filename
405	must be present.
406	"""
407
408	from anuga.geospatial_data.geospatial_data import Geospatial_data
409
410	# Treat special case: Polygon situation
411	# Location will be ignored and set to 'centroids'
412	# FIXME (Ole): This needs to be generalised and
413	# perhaps the notion of location and indices simplified
414
415	# FIXME (Ole): Need to compute indices based on polygon
416	# (and location) and use existing code after that.
417
418	# See ticket:275, ticket:250, ticeket:254 for refactoring plan
419
420	if polygon is not None:
421	if indices is not None:
422	msg = 'Only one of polygon and indices can be specified'
423	raise Exception, msg
424
425	msg = 'With polygon selected, set_quantity must provide '
426	msg += 'the keyword numeric and it must (currently) be '
427	msg += 'a constant.'
428	if numeric is None:
429	raise Exception, msg
430	else:
431	# Check that numeric is as constant
432	assert type(numeric) in [FloatType, IntType, LongType], msg
433
434	location = 'centroids'
435
436	points = self.domain.get_centroid_coordinates(absolute=True)
437	indices = inside_polygon(points, polygon)
438
439	self.set_values_from_constant(numeric, location, indices, verbose)
440
441	self.extrapolate_first_order()
442
443	if smooth:
444	self.smooth_vertex_values(use_cache=use_cache,
445	verbose=verbose)
446
447	return
448
449	# General input checks
450	L = [numeric, quantity, function, geospatial_data, filename]
451	msg = ('Exactly one of the arguments numeric, quantity, function, '
452	'geospatial_data, or filename must be present.')
453	assert L.count(None) == len(L)-1, msg
454
455	if location == 'edges':
456	msg = 'edges has been deprecated as valid location'
457	raise Exception, msg
458
459	if location not in ['vertices', 'centroids', 'unique vertices']:
460	msg = 'Invalid location: %s' % location
461	raise Exception, msg
462
463	msg = 'Indices must be a list, array or None'
464	assert isinstance(indices, (NoneType, list, num.ndarray)), msg
465
466	# Determine which 'set_values_from_...' to use
467	if numeric is not None:
468	if isinstance(numeric, (list, num.ndarray)):
469	self.set_values_from_array(numeric, location, indices,
470	use_cache=use_cache, verbose=verbose)
471	elif callable(numeric):
472	self.set_values_from_function(numeric, location, indices,
473	use_cache=use_cache,
474	verbose=verbose)
475	elif isinstance(numeric, Quantity):
476	self.set_values_from_quantity(numeric, location, indices,
477	verbose=verbose)
478	elif isinstance(numeric, Geospatial_data):
479	self.set_values_from_geospatial_data(numeric, alpha, location,
480	indices, verbose=verbose,
481	use_cache=use_cache)
482	else: # see if it's coercible to a float (float, int or long, etc)
483	try:
484	numeric = float(numeric)
485	except ValueError:
486	msg = ("Illegal type for variable 'numeric': %s"
487	% type(numeric))
488	raise Exception, msg
489	self.set_values_from_constant(numeric, location,
490	indices, verbose)
491	elif quantity is not None:
492	self.set_values_from_quantity(quantity, location, indices, verbose)
493	elif function is not None:
494	msg = 'Argument function must be callable'
495	assert callable(function), msg
496	self.set_values_from_function(function, location, indices,
497	use_cache=use_cache, verbose=verbose)
498	elif geospatial_data is not None:
499	self.set_values_from_geospatial_data(geospatial_data, alpha,
500	location, indices,
501	verbose=verbose,
502	use_cache=use_cache)
503	elif filename is not None:
504	if hasattr(self.domain, 'points_file_block_line_size'):
505	max_read_lines = self.domain.points_file_block_line_size
506	else:
507	max_read_lines = default_block_line_size
508	self.set_values_from_file(filename, attribute_name, alpha, location,
509	indices, verbose=verbose,
510	max_read_lines=max_read_lines,
511	use_cache=use_cache)
512	else:
513	raise Exception, "This can't happen :-)"
514
515	# Update all locations in triangles
516	if location == 'vertices' or location == 'unique vertices':
517	# Intialise centroid and edge_values
518	self.interpolate()
519
520	if location == 'centroids':
521	# Extrapolate 1st order - to capture notion of area being specified
522	self.extrapolate_first_order()
523
524	############################################################################
525	# Specific internal functions for setting values based on type
526	############################################################################
527
528	##
529	# @brief Set quantity values from specified constant.
530	# @param X The constant to set quantity values to.
531	# @param location
532	# @param indices
533	# @param verbose
534	def set_values_from_constant(self, X, location, indices, verbose):
535	"""Set quantity values from specified constant X"""
536
537	# FIXME (Ole): Somehow indices refer to centroids
538	# rather than vertices as default. See unit test
539	# test_set_vertex_values_using_general_interface_with_subset(self):
540
541	if location == 'centroids':
542	if indices is None:
543	self.centroid_values[:] = X
544	else:
545	# Brute force
546	for i in indices:
547	self.centroid_values[i] = X
548	elif location == 'unique vertices':
549	if indices is None:
550	self.edge_values[:] = X #FIXME (Ole): Shouldn't this be vertex_values?
551	else:
552	# Go through list of unique vertices
553	for unique_vert_id in indices:
554	triangles = \
555	self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
556
557	# In case there are unused points
558	if len(triangles) == 0:
559	continue
560
561	# Go through all triangle, vertex pairs
562	# and set corresponding vertex value
563	for triangle_id, vertex_id in triangles:
564	self.vertex_values[triangle_id, vertex_id] = X
565
566	# Intialise centroid and edge_values
567	self.interpolate()
568	else:
569	if indices is None:
570	self.vertex_values[:] = X
571	else:
572	# Brute force
573	for i_vertex in indices:
574	self.vertex_values[i_vertex] = X
575
576	##
577	# @brief Set values for a quantity.
578	# @param values Array of values.
579	# @param location Where values are to be stored.
580	# @param indices Limit update to these indices.
581	# @param use_cache ??
582	# @param verbose True if this method is to be verbose.
583	def set_values_from_array(self, values,
584	location='vertices',
585	indices=None,
586	use_cache=False,
587	verbose=False):
588	"""Set values for quantity
589
590	values: numeric array
591	location: Where values are to be stored.
592	Permissible options are: vertices, centroid, unique vertices
593	Default is 'vertices'
594
595	indices - if this action is carried out on a subset of
596	elements or unique vertices
597	The element/unique vertex indices are specified here.
598
599	In case of location == 'centroid' the dimension values must
600	be a list of a numerical array of length N, N being the number
601	of elements.
602
603	Otherwise it must be of dimension Nx3
604
605	The values will be stored in elements following their
606	internal ordering.
607
608	If selected location is vertices, values for centroid and edges
609	will be assigned interpolated values.
610	In any other case, only values for the specified locations
611	will be assigned and the others will be left undefined.
612	"""
613
614	values = num.array(values, num.float)
615
616	if indices is not None:
617	indices = num.array(indices, num.int)
618	msg = ('Number of values must match number of indices: You '
619	'specified %d values and %d indices'
620	% (values.shape[0], indices.shape[0]))
621	assert values.shape[0] == indices.shape[0], msg
622
623	N = self.centroid_values.shape[0]
624
625	if location == 'centroids':
626	assert len(values.shape) == 1, 'Values array must be 1d'
627
628	if indices is None:
629	msg = 'Number of values must match number of elements'
630	assert values.shape[0] == N, msg
631
632	self.centroid_values = values
633	else:
634	msg = 'Number of values must match number of indices'
635	assert values.shape[0] == indices.shape[0], msg
636
637	# Brute force
638	for i in range(len(indices)):
639	self.centroid_values[indices[i]] = values[i]
640	elif location == 'unique vertices':
641	msg = 'Values array must be 1d'
642	assert len(values.shape) == 1 or num.allclose(values.shape[1:], 1), msg
643
644	self.set_vertex_values(values.flatten(), indices=indices,
645	use_cache=use_cache, verbose=verbose)
646	else:
647	# Location vertices
648	if len(values.shape) == 1:
649	# This is the common case arising from fitted
650	# values (e.g. from pts file).
651	self.set_vertex_values(values, indices=indices,
652	use_cache=use_cache, verbose=verbose)
653	elif len(values.shape) == 2:
654	# Vertex values are given as a triplet for each triangle
655	msg = 'Array must be N x 3'
656	assert values.shape[1] == 3, msg
657
658	if indices is None:
659	self.vertex_values = values
660	else:
661	for element_index, value in map(None, indices, values):
662	self.vertex_values[element_index] = value
663	else:
664	msg = 'Values array must be 1d or 2d'
665	raise Exception, msg
666
667	##
668	# @brief Set quantity values from a specified quantity instance.
669	# @param q The quantity instance to take values from.
670	# @param location IGNORED, 'vertices' ALWAYS USED!
671	# @param indices ??
672	# @param verbose True if this method is to be verbose.
673	def set_values_from_quantity(self, q, location, indices, verbose):
674	"""Set quantity values from specified quantity instance q
675
676	Location is ignored - vertices will always be used here.
677	"""
678
679
680	A = q.vertex_values
681
682	msg = 'Quantities are defined on different meshes. '+\
683	'This might be a case for implementing interpolation '+\
684	'between different meshes.'
685	assert num.allclose(A.shape, self.vertex_values.shape), msg
686
687	self.set_values(A, location='vertices',
688	indices=indices, verbose=verbose)
689
690	##
691	# @brief Set quantity values from a specified quantity instance.
692	# @param f Callable that takes two 1d array -> 1d array.
693	# @param location Where values are to be stored.
694	# @param indices ??
695	# @param use_cache ??
696	# @param verbose True if this method is to be verbose.
697	def set_values_from_function(self,
698	f,
699	location='vertices',
700	indices=None,
701	use_cache=False,
702	verbose=False):
703	"""Set values for quantity using specified function
704
705	Input
706	f: x, y -> z Function where x, y and z are arrays
707	location: Where values are to be stored.
708	Permissible options are: vertices, centroid,
709	unique vertices
710	Default is "vertices"
711	indices:
712	"""
713
714	# FIXME: Should check that function returns something sensible and
715	# raise a meaningful exception if it returns None for example
716
717	# FIXME: Should supply absolute coordinates
718
719	# Compute the function values and call set_values again
720	if location == 'centroids':
721	if indices is None:
722	indices = range(len(self))
723
724	V = num.take(self.domain.get_centroid_coordinates(), indices, axis=0)
725	x = V[:,0]; y = V[:,1]
726	if use_cache is True:
727	res = cache(f, (x, y), verbose=verbose)
728	else:
729	res = f(x, y)
730
731	self.set_values(res, location=location, indices=indices)
732	elif location == 'vertices':
733	# This is the default branch taken by set_quantity
734	M = self.domain.number_of_triangles
735	V = self.domain.get_vertex_coordinates()
736
737	x = V[:,0];
738	y = V[:,1]
739	if use_cache is True:
740	values = cache(f, (x, y), verbose=verbose)
741	else:
742	if verbose is True:
743	log.critical('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	log.critical('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
1145	length of the indices list or 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	If indices is None (or omitted) all elements are returned as
1160	a copy of the relevant array. If performance is critical,
1161	use arrays domain.centroid_values, domain.vertex_values and
1162	domain.edge_values directly.
1163
1164	The values will be stored in elements following their
1165	internal ordering.
1166	"""
1167
1168	# FIXME (Ole): I reckon we should have the option of passing a
1169	# polygon into get_values. The question becomes how
1170	# resulting values should be ordered.
1171
1172	if verbose is True:
1173	log.critical('Getting values from %s' % location)
1174
1175	if interpolation_points is not None:
1176	return self.get_interpolated_values(interpolation_points,
1177	use_cache=use_cache,
1178	verbose=verbose)
1179
1180	# FIXME (Ole): Consider deprecating 'edges' - but not if it is used
1181	# elsewhere in ANUGA.
1182	# Edges have already been deprecated in set_values, see changeset:5521,
1183	# but might be useful in get_values. Any thoughts anyone?
1184	# YES (Ole): Edge values are necessary for volumetric balance
1185	# check and inflow boundary. Keep them!
1186
1187	if location not in ['vertices', 'centroids',
1188	'edges', 'unique vertices']:
1189	msg = 'Invalid location: %s' % location
1190	raise Exception, msg
1191
1192
1193	msg = '\'indices\' must be a list, array or None'
1194	assert isinstance(indices, (NoneType, list, num.ndarray)), msg
1195
1196	if location == 'centroids':
1197	if indices is None:
1198	return self.centroid_values.copy()
1199	else:
1200	return num.take(self.centroid_values, indices, axis=0)
1201	elif location == 'edges':
1202	if indices is None:
1203	return self.edge_values.copy()
1204	else:
1205	return num.take(self.edge_values, indices, axis=0)
1206	elif location == 'unique vertices':
1207	if indices is None:
1208	indices=range(self.domain.get_number_of_nodes())
1209	vert_values = []
1210
1211	# Go through list of unique vertices
1212	for unique_vert_id in indices:
1213	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1214
1215	# In case there are unused points
1216	if len(triangles) == 0:
1217	msg = 'Unique vertex not associated with triangles'
1218	raise Exception, msg
1219
1220	# Go through all triangle, vertex pairs
1221	# Average the values
1222	# FIXME (Ole): Should we merge this with get_vertex_values
1223	sum = 0
1224	for triangle_id, vertex_id in triangles:
1225	sum += self.vertex_values[triangle_id, vertex_id]
1226	vert_values.append(sum / len(triangles))
1227	return num.array(vert_values, num.float)
1228	else:
1229	if indices is None:
1230	return self.vertex_values.copy()
1231	else:
1232	return num.take(self.vertex_values, indices, axis=0)
1233
1234	##
1235	# @brief Set vertex values for all unique vertices based on array.
1236	# @param A Array to set values with.
1237	# @param indices Set of IDs of elements to work on.
1238	# @param use_cache ??
1239	# @param verbose??
1240	def set_vertex_values(self,
1241	A,
1242	indices=None,
1243	use_cache=False,
1244	verbose=False):
1245	"""Set vertex values for all unique vertices based on input array A
1246	which has one entry per unique vertex, i.e. one value for each row in
1247	array self.domain.nodes.
1248
1249	indices is the list of vertex_id's that will be set.
1250
1251	This function is used by set_values_from_array
1252	"""
1253
1254	# Check that A can be converted to array and is of appropriate dim
1255	A = ensure_numeric(A, num.float)
1256	assert len(A.shape) == 1
1257
1258	if indices is None:
1259	assert A.shape[0] == self.domain.get_nodes().shape[0]
1260	vertex_list = range(A.shape[0])
1261	else:
1262	assert A.shape[0] == len(indices)
1263	vertex_list = indices
1264
1265	#FIXME(Ole): This function ought to be faster.
1266	# We need to get the triangles_and_vertices list
1267	# from domain in one hit, then cache the computation of the
1268	# Nx3 array of vertex values that can then be assigned using
1269	# set_values_from_array.
1270	#
1271	# Alternatively, some C code would be handy
1272	#
1273	self._set_vertex_values(vertex_list, A)
1274
1275	##
1276	# @brief Go through list of unique vertices.
1277	# @param vertex_list ??
1278	# @param A ??
1279	def _set_vertex_values(self, vertex_list, A):
1280	"""Go through list of unique vertices
1281	This is the common case e.g. when values
1282	are obtained from a pts file through fitting
1283	"""
1284
1285	# Go through list of unique vertices
1286	for i_index, unique_vert_id in enumerate(vertex_list):
1287	triangles = self.domain.get_triangles_and_vertices_per_node(node=unique_vert_id)
1288
1289	# In case there are unused points
1290	if len(triangles) == 0:
1291	continue
1292
1293	# Go through all triangle, vertex pairs
1294	# touching vertex unique_vert_id and set corresponding vertex value
1295	for triangle_id, vertex_id in triangles:
1296	self.vertex_values[triangle_id, vertex_id] = A[i_index]
1297
1298	# Intialise centroid and edge_values
1299	self.interpolate()
1300
1301	##
1302	# @brief Smooth vertex values.
1303	def smooth_vertex_values(self, use_cache=False, verbose=False):
1304	"""Smooths vertex values."""
1305
1306	A, V = self.get_vertex_values(xy=False, smooth=True)
1307	self.set_vertex_values(A, use_cache=use_cache, verbose=verbose)
1308
1309	############################################################################
1310	# Methods for outputting model results
1311	############################################################################
1312
1313	##
1314	# @brief Get vertex values like an OBJ format i.e. one value per node.
1315	# @param xy True if we return X and Y as well as A and V.
1316	# @param smooth True if vertex values are to be smoothed.
1317	# @param precision The type of the result values (default float).
1318	def get_vertex_values(self, xy=True, smooth=None, precision=None):
1319	"""Return vertex values like an OBJ format i.e. one value per node.
1320
1321	The vertex values are returned as one sequence in the 1D float array A.
1322	If requested the coordinates will be returned in 1D arrays X and Y.
1323
1324	The connectivity is represented as an integer array, V, of dimension
1325	Mx3, where M is the number of triangles. Each row has three indices
1326	defining the triangle and they correspond to elements in the arrays
1327	X, Y and A.
1328
1329	If smooth is True, vertex values corresponding to one common coordinate
1330	set will be smoothed by taking the average of vertex values for each
1331	node. In this case vertex coordinates will be de-duplicated
1332	corresponding to the original nodes as obtained from the method
1333	general_mesh.get_nodes()
1334
1335	If no smoothings is required, vertex coordinates and values will be
1336	aggregated as a concatenation of values at vertices 0, vertices 1 and
1337	vertices 2. This corresponds to the node coordinates obtained from the
1338	method general_mesh.get_vertex_coordinates()
1339
1340	Calling convention
1341	if xy is True:
1342	X, Y, A, V = get_vertex_values
1343	else:
1344	A, V = get_vertex_values
1345	"""
1346
1347	if smooth is None:
1348	# Take default from domain
1349	try:
1350	smooth = self.domain.smooth
1351	except:
1352	smooth = False
1353
1354	if precision is None:
1355	precision = num.float
1356
1357	if smooth is True:
1358	# Ensure continuous vertex values by averaging values at each node
1359	V = self.domain.get_triangles()
1360	N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
1361	A = num.zeros(N, num.float)
1362	points = self.domain.get_nodes()
1363
1364	if 1:
1365	# Fast C version
1366	average_vertex_values(ensure_numeric(self.domain.vertex_value_indices),
1367	ensure_numeric(self.domain.number_of_triangles_per_node),
1368	ensure_numeric(self.vertex_values),
1369	A)
1370	A = A.astype(precision)
1371	else:
1372	# Slow Python version
1373	current_node = 0
1374	k = 0 # Track triangles touching on node
1375	total = 0.0
1376	for index in self.domain.vertex_value_indices:
1377	if current_node == N:
1378	msg = 'Current node exceeding number of nodes (%d) ' % N
1379	raise Exception, msg
1380
1381	k += 1
1382
1383	volume_id = index / 3
1384	vertex_id = index % 3
1385
1386	v = self.vertex_values[volume_id, vertex_id]
1387	total += v
1388
1389	if self.domain.number_of_triangles_per_node[current_node] == k:
1390	A[current_node] = total/k
1391
1392	# Move on to next node
1393	total = 0.0
1394	k = 0
1395	current_node += 1
1396	else:
1397	# Return disconnected internal vertex values
1398	V = self.domain.get_disconnected_triangles()
1399	points = self.domain.get_vertex_coordinates()
1400	A = self.vertex_values.flatten().astype(precision)
1401
1402	# Return
1403	if xy is True:
1404	X = points[:,0].astype(precision)
1405	Y = points[:,1].astype(precision)
1406
1407	return X, Y, A, V
1408	else:
1409	return A, V
1410
1411	##
1412	# @brief Extrapolate conserved quantities from centroid.
1413	def extrapolate_first_order(self):
1414	"""Extrapolate conserved quantities from centroid to vertices and edges
1415	for each volume using first order scheme.
1416	"""
1417
1418	qc = self.centroid_values
1419	qv = self.vertex_values
1420	qe = self.edge_values
1421
1422	for i in range(3):
1423	qv[:,i] = qc
1424	qe[:,i] = qc
1425
1426	self.x_gradient *= 0.0
1427	self.y_gradient *= 0.0
1428
1429	##
1430	# @brief Compute the integral of quantity across entire domain.
1431	# @return The integral.
1432	def get_integral(self):
1433	"""Compute the integral of quantity across entire domain."""
1434
1435
1436	areas = self.domain.get_areas()
1437	"""
1438	integral = 0
1439	for k in range(len(self.domain)):
1440	area = areas[k]
1441	qc = self.centroid_values[k]
1442	integral += qc*area
1443	"""
1444	return num.sum(areas*self.centroid_values)
1445
1446	#return integral
1447
1448	##
1449	# @brief get the gradients.
1450	def get_gradients(self):
1451	"""Provide gradients. Use compute_gradients first."""
1452
1453	return self.x_gradient, self.y_gradient
1454
1455	##
1456	# @brief ??
1457	# @param timestep ??
1458	def update(self, timestep):
1459	# Call correct module function
1460	# (either from this module or C-extension)
1461	return update(self, timestep)
1462
1463	##
1464	# @brief ??
1465	def compute_gradients(self):
1466	# Call correct module function
1467	# (either from this module or C-extension)
1468	return compute_gradients(self)
1469
1470	##
1471	# @brief ??
1472	def limit(self):
1473	# Call correct module depending on whether
1474	# basing limit calculations on edges or vertices
1475	limit_old(self)
1476
1477	##
1478	# @brief ??
1479	def limit_vertices_by_all_neighbours(self):
1480	# Call correct module function
1481	# (either from this module or C-extension)
1482	limit_vertices_by_all_neighbours(self)
1483
1484	##
1485	# @brief ??
1486	def limit_edges_by_all_neighbours(self):
1487	# Call correct module function
1488	# (either from this module or C-extension)
1489	limit_edges_by_all_neighbours(self)
1490
1491	##
1492	# @brief ??
1493	def limit_edges_by_neighbour(self):
1494	# Call correct module function
1495	# (either from this module or C-extension)
1496	limit_edges_by_neighbour(self)
1497
1498	##
1499	# @brief ??
1500	def extrapolate_second_order(self):
1501	# Call correct module function
1502	# (either from this module or C-extension)
1503	compute_gradients(self)
1504	extrapolate_from_gradient(self)
1505
1506	##
1507	# @brief ??
1508	def extrapolate_second_order_and_limit_by_edge(self):
1509	# Call correct module function
1510	# (either from this module or C-extension)
1511	extrapolate_second_order_and_limit_by_edge(self)
1512
1513	##
1514	# @brief ??
1515	def extrapolate_second_order_and_limit_by_vertex(self):
1516	# Call correct module function
1517	# (either from this module or C-extension)
1518	extrapolate_second_order_and_limit_by_vertex(self)
1519
1520	##
1521	# @brief ??
1522	# @param bound ??
1523	def bound_vertices_below_by_constant(self, bound):
1524	# Call correct module function
1525	# (either from this module or C-extension)
1526	bound_vertices_below_by_constant(self, bound)
1527
1528	##
1529	# @brief ??
1530	# @param quantity ??
1531	def bound_vertices_below_by_quantity(self, quantity):
1532	# Call correct module function
1533	# (either from this module or C-extension)
1534
1535	# check consistency
1536	assert self.domain == quantity.domain
1537	bound_vertices_below_by_quantity(self, quantity)
1538
1539	##
1540	# @brief ??
1541	def backup_centroid_values(self):
1542	# Call correct module function
1543	# (either from this module or C-extension)
1544	backup_centroid_values(self)
1545
1546	##
1547	# @brief ??
1548	# @param a ??
1549	# @param b ??
1550	def saxpy_centroid_values(self, a, b):
1551	# Call correct module function
1552	# (either from this module or C-extension)
1553	saxpy_centroid_values(self, a, b)
1554
1555
1556	##
1557	# @brief OBSOLETE!
1558	class Conserved_quantity(Quantity):
1559	"""Class conserved quantity being removed, use Quantity."""
1560
1561	def __init__(self, domain, vertex_values=None):
1562	msg = 'ERROR: Use Quantity instead of Conserved_quantity'
1563	raise Exception, msg
1564
1565
1566	######
1567	# Prepare the C extensions.
1568	######
1569
1570	from anuga.utilities import compile
1571
1572	if compile.can_use_C_extension('quantity_ext.c'):
1573	# Underlying C implementations can be accessed
1574
1575	from quantity_ext import \
1576	average_vertex_values,\
1577	backup_centroid_values,\
1578	saxpy_centroid_values,\
1579	compute_gradients,\
1580	limit_old,\
1581	limit_vertices_by_all_neighbours,\
1582	limit_edges_by_all_neighbours,\
1583	limit_edges_by_neighbour,\
1584	limit_gradient_by_neighbour,\
1585	extrapolate_from_gradient,\
1586	extrapolate_second_order_and_limit_by_edge,\
1587	extrapolate_second_order_and_limit_by_vertex,\
1588	bound_vertices_below_by_constant,\
1589	bound_vertices_below_by_quantity,\
1590	interpolate_from_vertices_to_edges,\
1591	interpolate_from_edges_to_vertices,\
1592	update
1593	else:
1594	msg = 'C implementations could not be accessed by %s.\n ' % __file__
1595	msg += 'Make sure compile_all.py has been run as described in '
1596	msg += 'the ANUGA installation guide.'
1597	raise Exception, msg

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

Download in other formats: