Context Navigation

source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py @ 7776

Last change on this file since 7776 was 7737, checked in by hudson, 14 years ago
Various refactorings, all unit tests pass. Domain renamed to generic domain.
File size: 59.1 KB

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

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

Download in other formats: