Context Navigation

source: inundation-numpy-branch/pyvolution/quantity.py @ 2608

Last change on this file since 2608 was 2608, checked in by ole, 18 years ago
Played with custom exceptions for ANUGA
File size: 43.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 numpy import concatenate, zeros, Float, Int, array, reshape, ArrayType, allclose, take, sum, equal
18
19
20	class Quantity:
21
22	def __init__(self, domain, vertex_values=None):
23
24	from mesh import Mesh
25	#from Numeric import array, zeros, Float
26
27	msg = 'First argument in Quantity.__init__ '
28	msg += 'must be of class Mesh (or a subclass thereof)'
29	assert isinstance(domain, Mesh), msg
30
31	if vertex_values is None:
32	N = domain.number_of_elements
33	self.vertex_values = zeros((N, 3), Float)
34	else:
35	self.vertex_values = array(vertex_values).astype(Float)
36
37	N, V = self.vertex_values.shape
38	assert V == 3,\
39	'Three vertex values per element must be specified'
40
41
42	msg = 'Number of vertex values (%d) must be consistent with'\
43	%N
44	msg += 'number of elements in specified domain (%d).'\
45	%domain.number_of_elements
46
47	assert N == domain.number_of_elements, msg
48
49	self.domain = domain
50
51	#Allocate space for other quantities
52	self.centroid_values = zeros(N, Float)
53	self.edge_values = zeros((N, 3), Float)
54
55	#Intialise centroid and edge_values
56	self.interpolate()
57
58
59
60	#Methods for operator overloading
61	def __len__(self):
62	return self.centroid_values.shape[0]
63
64
65	def __neg__(self):
66	"""Negate all values in this quantity giving meaning to the
67	expression -Q where Q is an instance of class Quantity
68	"""
69
70	Q = Quantity(self.domain)
71	Q.set_values(-self.vertex_values)
72	return Q
73
74
75	def __add__(self, other):
76	"""Add to self anything that could populate a quantity
77
78	E.g other can be a constant, an array, a function, another quantity
79	(except for a filename or points, attributes (for now))
80	- see set_values for details
81	"""
82
83	Q = Quantity(self.domain)
84	Q.set_values(other)
85
86	result = Quantity(self.domain)
87	result.set_values(self.vertex_values + Q.vertex_values)
88	return result
89
90	def __radd__(self, other):
91	"""Handle cases like 7+Q, where Q is an instance of class Quantity
92	"""
93	return self + other
94
95
96	def __sub__(self, other):
97	return self + -other #Invoke __neg__
98
99	def __mul__(self, other):
100	"""Multiply self with anything that could populate a quantity
101
102	E.g other can be a constant, an array, a function, another quantity
103	(except for a filename or points, attributes (for now))
104	- see set_values for details
105
106	Note that if two quantitites q1 and q2 are multiplied,
107	vertex values are multiplied entry by entry
108	while centroid and edge values are re-interpolated.
109	Hence they won't be the product of centroid or edge values
110	from q1 and q2.
111	"""
112
113	Q = Quantity(self.domain)
114	Q.set_values(other)
115
116	result = Quantity(self.domain)
117	result.set_values(self.vertex_values * Q.vertex_values)
118	return result
119
120	def __rmul__(self, other):
121	"""Handle cases like 3*Q, where Q is an instance of class Quantity
122	"""
123	return self * other
124
125	def __pow__(self, other):
126	"""Raise quantity to (numerical) power
127
128	As with __mul__ vertex values are processed entry by entry
129	while centroid and edge values are re-interpolated.
130
131	Example using __pow__:
132	Q = (Q12 + Q22)**0.5
133
134	"""
135
136	result = Quantity(self.domain)
137	result.set_values(self.vertex_values**other)
138	return result
139
140
141
142	def interpolate(self):
143	"""Compute interpolated values at edges and centroid
144	Pre-condition: vertex_values have been set
145	"""
146
147	N = self.vertex_values.shape[0]
148	for i in range(N):
149	v0 = self.vertex_values[i, 0]
150	v1 = self.vertex_values[i, 1]
151	v2 = self.vertex_values[i, 2]
152
153	self.centroid_values[i] = (v0 + v1 + v2)/3
154
155	self.interpolate_from_vertices_to_edges()
156
157
158	def interpolate_from_vertices_to_edges(self):
159	#Call correct module function
160	#(either from this module or C-extension)
161	interpolate_from_vertices_to_edges(self)
162
163
164
165
166	#New leaner interface to setting values
167	def set_values(self,
168	numeric = None, # List, numeric array or constant
169	quantity = None, # Another quantity
170	function = None, # Callable object: f(x,y)
171	geospatial_data = None, #Arbitrary dataset
172	points = None, values = None, data_georef = None, #Input for least squares (obsoleted by use of geo_spatial object)
173	filename = None, attribute_name = None, #Input from file
174	alpha = None,
175	location = 'vertices',
176	indices = None,
177	verbose = False,
178	use_cache = False):
179
180	"""Set values for quantity based on different sources.
181
182	numeric:
183	Compatible list, Numeric array (see below) or constant.
184	If callable it will treated as a function (see below)
185	If instance of another Quantity it will be treated as such.
186	If geo_spatial object it will be treated as such
187
188	quantity:
189	Another quantity (compatible quantity, e.g. obtained as a
190	linear combination of quantities)
191
192	function:
193	Any callable object that takes two 1d arrays x and y
194	each of length N and returns an array also of length N.
195	The function will be evaluated at points determined by
196	location and indices in the underlying mesh.
197
198	geospatial_data:
199	Arbitrary geo spatial dataset in the form of the class
200	Geospatial_data. Mesh points are populated using least squares
201	fitting
202
203	points:
204	Nx2 array of data points for use with least squares fit
205	If points are present, an N array of attribute
206	values corresponding to
207	each data point must be present.
208
209	values:
210	If points is specified, values is an array of length N containing
211	attribute values for each point.
212
213	data_georef:
214	If points is specified, geo_reference applies to each point.
215
216	filename:
217	Name of a .pts file containing data points and attributes for
218	use with least squares.
219
220	attribute_name:
221	If specified, any array matching that name
222	will be used. from file or geospatial_data.
223	Otherwise a default will be used.
224
225	alpha:
226	Smoothing parameter to be used with least squares fits.
227	See module least_squares for further details about alpha.
228	Alpha will only be used with points, values or filename.
229	Otherwise it will be ignored.
230
231
232	location: Where values are to be stored.
233	Permissible options are: vertices, edges, centroids
234	Default is 'vertices'
235
236	In case of location == 'centroids' the dimension values must
237	be a list of a Numerical array of length N,
238	N being the number of elements.
239	Otherwise it must be of dimension Nx3
240
241
242	The values will be stored in elements following their
243	internal ordering.
244
245	If location is not 'unique vertices' Indices is the
246	set of element ids that the operation applies to.
247	If location is 'unique vertices' Indices is the set
248	of vertex ids that the operation applies to.
249
250	If selected location is vertices, values for
251	centroid and edges will be assigned interpolated
252	values. In any other case, only values for the
253	specified locations will be assigned and the others
254	will be left undefined.
255
256	verbose: True means that output to stdout is generated
257
258	use_cache: True means that caching of intermediate results is
259	attempted for least squares fit.
260
261
262
263
264	Exactly one of the arguments
265	numeric, quantity, function, points, filename
266	must be present.
267	"""
268
269	from geospatial_data.geospatial_data import Geospatial_data
270	from types import FloatType, IntType, LongType, ListType, NoneType
271
272	#General input checks
273	L = [numeric, quantity, function, geospatial_data, points, filename]
274	msg = 'Exactly one of the arguments '+\
275	'numeric, quantity, function, geospatial_data, points, or filename '+\
276	'must be present.'
277	assert L.count(None) == len(L)-1, msg
278
279
280	if location not in ['vertices', 'centroids', 'edges',
281	'unique vertices']:
282	msg = 'Invalid location: %s' %location
283	raise msg
284
285
286	msg = 'Indices must be a list or None'
287	assert type(indices) in [ListType, NoneType, ArrayType], msg
288
289
290	if not(points is None and values is None and data_georef is None):
291	from warnings import warn
292
293	msg = 'Using points, values or data_georef with set_quantity '
294	msg += 'is obsolete. Please use a Geospatial_data object instead.'
295	warn(msg, DeprecationWarning)
296
297
298
299	#Determine which 'set_values_from_...' to use
300
301	if numeric is not None:
302	if type(numeric) in [FloatType, IntType, LongType]:
303	self.set_values_from_constant(numeric,
304	location, indices, verbose)
305	elif type(numeric) in [ArrayType, ListType]:
306	self.set_values_from_array(numeric,
307	location, indices, verbose)
308	elif callable(numeric):
309	self.set_values_from_function(numeric,
310	location, indices, verbose)
311	elif isinstance(numeric, Quantity):
312	self.set_values_from_quantity(numeric,
313	location, indices, verbose)
314	elif isinstance(numeric, Geospatial_data):
315	self.set_values_from_geospatial_data(numeric,
316	alpha,
317	location, indices,
318	verbose = verbose,
319	use_cache = use_cache)
320	else:
321	msg = 'Illegal type for argument numeric: %s' %str(numeric)
322	raise msg
323
324	elif quantity is not None:
325	self.set_values_from_quantity(quantity,
326	location, indices, verbose)
327	elif function is not None:
328	msg = 'Argument function must be callable'
329	assert callable(function), msg
330	self.set_values_from_function(function,
331	location, indices, verbose)
332	elif geospatial_data is not None:
333	self.set_values_from_geospatial_data(geospatial_data,
334	alpha,
335	location, indices,
336	verbose = verbose,
337	use_cache = use_cache)
338	elif points is not None:
339	print 'The usage of points in set_values will be deprecated.' +\
340	'Please use the geospatial_data object.'
341
342	msg = 'When points are specified, associated values must also be.'
343	assert values is not None, msg
344	self.set_values_from_points(points, values, alpha,
345	location, indices,
346	data_georef = data_georef,
347	verbose = verbose,
348	use_cache = use_cache)
349	elif filename is not None:
350	self.set_values_from_file(filename, attribute_name, alpha,
351	location, indices,
352	verbose = verbose,
353	use_cache = use_cache)
354	else:
355	raise 'This can\'t happen :-)'
356
357
358	#Update all locations in triangles
359	if location == 'vertices' or location == 'unique vertices':
360	#Intialise centroid and edge_values
361	self.interpolate()
362
363	if location == 'centroids':
364	#Extrapolate 1st order - to capture notion of area being specified
365	self.extrapolate_first_order()
366
367
368
369	#Specific functions for setting values
370	def set_values_from_constant(self, X,
371	location, indices, verbose):
372	"""Set quantity values from specified constant X
373	"""
374
375
376	if location == 'centroids':
377	if (indices == None):
378	self.centroid_values[:] = X
379	else:
380	#Brute force
381	for i in indices:
382	self.centroid_values[i,:] = X
383
384	elif location == 'edges':
385	if (indices == None):
386	self.edge_values[:] = X
387	else:
388	#Brute force
389	for i in indices:
390	self.edge_values[i,:] = X
391
392	elif location == 'unique vertices':
393	if (indices == None):
394	self.edge_values[:] = X
395	else:
396
397	#Go through list of unique vertices
398	for unique_vert_id in indices:
399	triangles = self.domain.vertexlist[unique_vert_id]
400
401	#In case there are unused points
402	if triangles is None: continue
403
404	#Go through all triangle, vertex pairs
405	#and set corresponding vertex value
406	for triangle_id, vertex_id in triangles:
407	self.vertex_values[triangle_id, vertex_id] = X
408
409	#Intialise centroid and edge_values
410	self.interpolate()
411	else:
412	if (indices == None):
413	self.vertex_values[:] = X
414	else:
415	#Brute force
416	for i_vertex in indices:
417	self.vertex_values[i_vertex,:] = X
418
419
420
421
422
423
424	def set_values_from_array(self, values,
425	location, indices, verbose):
426	"""Set values for quantity
427
428	values: Numeric array
429	location: Where values are to be stored.
430	Permissible options are: vertices, edges, centroid, unique vertices
431	Default is 'vertices'
432
433	indices - if this action is carried out on a subset of
434	elements or unique vertices
435	The element/unique vertex indices are specified here.
436
437	In case of location == 'centroid' the dimension values must
438	be a list of a Numerical array of length N, N being the number
439	of elements.
440
441	Otherwise it must be of dimension Nx3
442
443	The values will be stored in elements following their
444	internal ordering.
445
446	If selected location is vertices, values for centroid and edges
447	will be assigned interpolated values.
448	In any other case, only values for the specified locations
449	will be assigned and the others will be left undefined.
450	"""
451
452	values = array(values).astype(Float)
453
454	if indices is not None:
455	indices = array(indices).astype(Int)
456	msg = 'Number of values must match number of indices'
457	assert values.shape[0] == indices.shape[0], msg
458
459	N = self.centroid_values.shape[0]
460
461	if location == 'centroids':
462	assert len(values.shape) == 1, 'Values array must be 1d'
463
464	if indices is None:
465	msg = 'Number of values must match number of elements'
466	assert values.shape[0] == N, msg
467
468	self.centroid_values = values
469	else:
470	msg = 'Number of values must match number of indices'
471	assert values.shape[0] == indices.shape[0], msg
472
473	#Brute force
474	for i in range(len(indices)):
475	self.centroid_values[indices[i]] = values[i]
476
477	elif location == 'edges':
478	assert len(values.shape) == 2, 'Values array must be 2d'
479
480	msg = 'Number of values must match number of elements'
481	assert values.shape[0] == N, msg
482
483	msg = 'Array must be N x 3'
484	assert values.shape[1] == 3, msg
485
486	self.edge_values = values
487
488	elif location == 'unique vertices':
489	assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
490	'Values array must be 1d'
491
492	self.set_vertex_values(values.flat, indices = indices)
493	else:
494	if len(values.shape) == 1:
495	self.set_vertex_values(values, indices = indices)
496	#if indices == None:
497	#Values are being specified once for each unique vertex
498	# msg = 'Number of values must match number of vertices'
499	# assert values.shape[0] == self.domain.coordinates.shape[0], msg
500	# self.set_vertex_values(values)
501	#else:
502	# for element_index, value in map(None, indices, values):
503	# self.vertex_values[element_index, :] = value
504
505	elif len(values.shape) == 2:
506	#Vertex values are given as a triplet for each triangle
507
508	msg = 'Array must be N x 3'
509	assert values.shape[1] == 3, msg
510
511	if indices == None:
512	self.vertex_values = values
513	else:
514	for element_index, value in map(None, indices, values):
515	self.vertex_values[element_index] = value
516	else:
517	msg = 'Values array must be 1d or 2d'
518	raise msg
519
520	def set_values_from_quantity(self, q,
521	location, indices, verbose):
522	"""Set quantity values from specified quantity instance q
523
524	Location is ignored
525	"""
526
527
528	A = q.vertex_values
529
530	msg = 'Quantities are defined on different meshes. '+\
531	'This might be a case for implementing interpolation '+\
532	'between different meshes.'
533	assert allclose(A.shape, self.vertex_values.shape), msg
534
535	self.set_values(A, location='vertices',
536	indices=indices,
537	verbose=verbose)
538
539
540	def set_values_from_function(self, f,
541	location, indices, verbose):
542	"""Set values for quantity using specified function
543
544	f: x, y -> z Function where x, y and z are arrays
545	location: Where values are to be stored.
546	Permissible options are: vertices, centroid, edges,
547	unique vertices
548	Default is "vertices"
549	"""
550
551	#FIXME: Should check that function returns something sensible and
552	#raise a meaningfull exception if it returns None for example
553
554	#FIXME: Should supply absolute coordinates
555
556	if (indices is None):
557	indices = range(len(self))
558	is_subset = False
559	else:
560	is_subset = True
561
562	if location == 'centroids':
563	P = take(self.domain.centroid_coordinates, indices)
564	if is_subset:
565	self.set_values(f(P[:,0], P[:,1]),
566	location = location,
567	indices = indices)
568	else:
569	self.set_values(f(P[:,0], P[:,1]), location = location)
570	elif location == 'vertices':
571	P = self.domain.vertex_coordinates
572	if is_subset:
573	#Brute force
574	for e in indices:
575	for i in range(3):
576	self.vertex_values[e,i] = f(P[e,2i], P[e,2i+1])
577	else:
578	for i in range(3):
579	self.vertex_values[:,i] = f(P[:,2i], P[:,2i+1])
580	else:
581	raise 'Not implemented: %s' %location
582
583
584
585	def set_values_from_geospatial_data(self, geospatial_data, alpha,
586	location, indices,
587	verbose = False,
588	use_cache = False):
589
590	#FIXME: Use this function for the time being. Later move code in here
591
592	points = geospatial_data.get_data_points(absolute = False)
593	values = geospatial_data.get_attributes()
594	data_georef = geospatial_data.get_geo_reference()
595
596
597
598	self.set_values_from_points(points, values, alpha,
599	location, indices,
600	data_georef = data_georef,
601	verbose = verbose,
602	use_cache = use_cache)
603
604
605
606	def set_values_from_points(self, points, values, alpha,
607	location, indices,
608	data_georef = None,
609	verbose = False,
610	use_cache = False):
611	"""Set quantity values from arbitray data points using least squares
612	"""
613
614
615	from utilities.numerical_tools import ensure_numeric
616	from least_squares import fit_to_mesh
617	from coordinate_transforms.geo_reference import Geo_reference
618
619
620	points = ensure_numeric(points, Float)
621	values = ensure_numeric(values, Float)
622
623	if location != 'vertices':
624	msg = 'set_values_from_points is only defined for '+\
625	'location=\'vertices\''
626	raise msg
627
628	coordinates = self.domain.coordinates
629	triangles = self.domain.triangles
630
631
632	#Take care of georeferencing
633	if data_georef is None:
634	data_georef = Geo_reference()
635
636
637	mesh_georef = self.domain.geo_reference
638
639	#print mesh_georef
640	#print data_georef
641	#print points
642
643
644	#Call least squares method
645	args = (coordinates, triangles, points, values)
646	kwargs = {'data_origin': data_georef.get_origin(),
647	'mesh_origin': mesh_georef.get_origin(),
648	'alpha': alpha,
649	'precrop': True,
650	'verbose': verbose}
651
652	#print kwargs
653
654	if use_cache is True:
655	try:
656	from caching import cache
657	except:
658	msg = 'Caching was requested, but caching module'+\
659	'could not be imported'
660	raise msg
661
662	vertex_attributes = cache(fit_to_mesh,
663	args, kwargs,
664	verbose = verbose)
665	else:
666
667	vertex_attributes = apply(fit_to_mesh,
668	args, kwargs)
669
670	#Call underlying method using array values
671	self.set_values_from_array(vertex_attributes,
672	location, indices, verbose)
673
674
675
676	def set_values_from_file(self, filename, attribute_name, alpha,
677	location, indices,
678	verbose = False,
679	use_cache = False):
680	"""Set quantity based on arbitrary points in .pts file
681	using least_squares attribute_name selects name of attribute
682	present in file.
683	If not specified try to use whatever is available in file.
684	"""
685
686	from load_mesh.loadASCII import import_points_file
687	from geospatial_data.geospatial_data import points_dictionary2geospatial_data
688
689	from types import StringType
690	msg = 'Filename must be a text string'
691	assert type(filename) == StringType, msg
692
693
694	#Read from (NetCDF) file
695	#FIXME (Ole): This function should really return a Geospatial_data object.
696	points_dict = import_points_file(filename)
697	points = points_dict['pointlist']
698	attributes = points_dict['attributelist']
699
700	if attribute_name is None:
701	names = attributes.keys()
702	attribute_name = names[0]
703
704	msg = 'Attribute_name must be a text string'
705	assert type(attribute_name) == StringType, msg
706
707
708	if verbose:
709	print 'Using attribute %s from file %s' %(attribute_name, filename)
710	print 'Available attributes: %s' %(names)
711
712	#try:
713	# z = attributes[attribute_name]
714	#except:
715	# msg = 'Could not extract attribute %s from file %s'\
716	# %(attribute_name, filename)
717	# raise msg
718
719
720	#Take care of georeferencing
721	if points_dict.has_key('geo_reference') and \
722	points_dict['geo_reference'] is not None:
723	data_georef = points_dict['geo_reference']
724	else:
725	data_georef = None
726
727
728
729	#Call underlying method for geospatial data
730	geospatial_data = points_dictionary2geospatial_data(points_dict)
731	geospatial_data.set_default_attribute_name(attribute_name)
732
733	self.set_values_from_geospatial_data(geospatial_data,
734	alpha,
735	location, indices,
736	verbose = verbose,
737	use_cache = use_cache)
738
739	#Call underlying method for points
740	#self.set_values_from_points(points, z, alpha,
741	# location, indices,
742	# data_georef = data_georef,
743	# verbose = verbose,
744	# use_cache = use_cache)
745
746
747
748	def get_values(self, location='vertices', indices = None):
749	"""get values for quantity
750
751	return X, Compatible list, Numeric array (see below)
752	location: Where values are to be stored.
753	Permissible options are: vertices, edges, centroid
754	and unique vertices. Default is 'vertices'
755
756	In case of location == 'centroids' the dimension values must
757	be a list of a Numerical array of length N, N being the number
758	of elements. Otherwise it must be of dimension Nx3
759
760	The returned values with be a list the length of indices
761	(N if indices = None). Each value will be a list of the three
762	vertex values for this quantity.
763
764	Indices is the set of element ids that the operation applies to.
765
766	The values will be stored in elements following their
767	internal ordering.
768
769	"""
770
771	if location not in ['vertices', 'centroids', 'edges', 'unique vertices']:
772	msg = 'Invalid location: %s' %location
773	raise msg
774
775	import types
776	assert type(indices) in [types.ListType, types.NoneType,
777	ArrayType],\
778	'Indices must be a list or None'
779
780	if location == 'centroids':
781	if (indices == None):
782	indices = range(len(self))
783	return take(self.centroid_values,indices)
784	elif location == 'edges':
785	if (indices == None):
786	indices = range(len(self))
787	return take(self.edge_values,indices)
788	elif location == 'unique vertices':
789	if (indices == None):
790	indices=range(self.domain.coordinates.shape[0])
791	vert_values = []
792	#Go through list of unique vertices
793	for unique_vert_id in indices:
794	triangles = self.domain.vertexlist[unique_vert_id]
795
796	#In case there are unused points
797	if triangles is None:
798	msg = 'Unique vertex not associated with triangles'
799	raise msg
800
801	# Go through all triangle, vertex pairs
802	# Average the values
803	sum = 0
804	for triangle_id, vertex_id in triangles:
805	sum += self.vertex_values[triangle_id, vertex_id]
806	vert_values.append(sum/len(triangles))
807	return array(vert_values)
808	else:
809	if (indices == None):
810	indices = range(len(self))
811	return take(self.vertex_values,indices)
812
813
814
815	def set_vertex_values(self, A, indices = None):
816	"""Set vertex values for all unique vertices based on input array A
817	which has one entry per unique vertex, i.e.
818	one value for each row in array self.domain.coordinates or
819	one value for each row in vertexlist.
820
821	indices is the list of vertex_id's that will be set.
822
823	This function is used by set_values_from_array
824	"""
825
826	#Assert that A can be converted to a Numeric array of appropriate dim
827	A = array(A, Float)
828
829	#print 'SHAPE A', A.shape
830	assert len(A.shape) == 1
831
832	if indices == None:
833	assert A.shape[0] == self.domain.coordinates.shape[0]
834	vertex_list = range(A.shape[0])
835	else:
836	assert A.shape[0] == len(indices)
837	vertex_list = indices
838
839	#Go through list of unique vertices
840	for i_index, unique_vert_id in enumerate(vertex_list):
841	triangles = self.domain.vertexlist[unique_vert_id]
842
843	if triangles is None: continue #In case there are unused points
844
845	#Go through all triangle, vertex pairs
846	#touching vertex unique_vert_id and set corresponding vertex value
847	for triangle_id, vertex_id in triangles:
848	self.vertex_values[triangle_id, vertex_id] = A[i_index]
849
850	#Intialise centroid and edge_values
851	self.interpolate()
852
853
854	def smooth_vertex_values(self, value_array='field_values',
855	precision = None):
856	""" Smooths field_values or conserved_quantities data.
857	TODO: be able to smooth individual fields
858	NOTE: This function does not have a test.
859	FIXME: NOT DONE - do we need it?
860	FIXME: this function isn't called by anything.
861	Maybe it should be removed..-DSG
862	"""
863
864
865	A,V = self.get_vertex_values(xy=False,
866	value_array=value_array,
867	smooth = True,
868	precision = precision)
869
870	#Set some field values
871	for volume in self:
872	for i,v in enumerate(volume.vertices):
873	if value_array == 'field_values':
874	volume.set_field_values('vertex', i, A[v,:])
875	elif value_array == 'conserved_quantities':
876	volume.set_conserved_quantities('vertex', i, A[v,:])
877
878	if value_array == 'field_values':
879	self.precompute()
880	elif value_array == 'conserved_quantities':
881	Volume.interpolate_conserved_quantities()
882
883
884	#Method for outputting model results
885	#FIXME: Split up into geometric and numeric stuff.
886	#FIXME: Geometric (X,Y,V) should live in mesh.py
887	#FIXME: STill remember to move XY to mesh
888	def get_vertex_values(self,
889	xy=True,
890	smooth = None,
891	precision = None,
892	reduction = None):
893	"""Return vertex values like an OBJ format
894
895	The vertex values are returned as one sequence in the 1D float array A.
896	If requested the coordinates will be returned in 1D arrays X and Y.
897
898	The connectivity is represented as an integer array, V, of dimension
899	M x 3, where M is the number of volumes. Each row has three indices
900	into the X, Y, A arrays defining the triangle.
901
902	if smooth is True, vertex values corresponding to one common
903	coordinate set will be smoothed according to the given
904	reduction operator. In this case vertex coordinates will be
905	de-duplicated.
906
907	If no smoothings is required, vertex coordinates and values will
908	be aggregated as a concatenation of values at
909	vertices 0, vertices 1 and vertices 2
910
911
912	Calling convention
913	if xy is True:
914	X,Y,A,V = get_vertex_values
915	else:
916	A,V = get_vertex_values
917
918	"""
919
920	if smooth is None:
921	smooth = self.domain.smooth
922
923	if precision is None:
924	precision = Float
925
926	if reduction is None:
927	reduction = self.domain.reduction
928
929	#Create connectivity
930
931	if smooth == True:
932
933	V = self.domain.get_vertices()
934	N = len(self.domain.vertexlist)
935	A = zeros(N, precision)
936
937	#Smoothing loop
938	for k in range(N):
939	L = self.domain.vertexlist[k]
940
941	#Go through all triangle, vertex pairs
942	#contributing to vertex k and register vertex value
943
944	if L is None: continue #In case there are unused points
945
946	contributions = []
947	for volume_id, vertex_id in L:
948	v = self.vertex_values[volume_id, vertex_id]
949	contributions.append(v)
950
951	A[k] = reduction(contributions)
952
953
954	if xy is True:
955	X = self.domain.coordinates[:,0].astype(precision)
956	Y = self.domain.coordinates[:,1].astype(precision)
957
958	return X, Y, A, V
959	else:
960	return A, V
961	else:
962	#Don't smooth
963	#obj machinery moved to general_mesh
964
965	# Create a V like [[0 1 2], [3 4 5]....[3m-2 3m-1 3*m]]
966	# These vert_id's will relate to the verts created below
967	#m = len(self.domain) #Number of volumes
968	#M = 3*m #Total number of unique vertices
969	#V = reshape(array(range(M)).astype(Int), (m,3))
970
971	V = self.domain.get_triangles(obj=True)
972	#FIXME use get_vertices, when ready
973
974	A = self.vertex_values.flat
975
976	#Do vertex coordinates
977	if xy is True:
978	C = self.domain.get_vertex_coordinates()
979
980	X = C[:,0:6:2].copy()
981	Y = C[:,1:6:2].copy()
982
983	return X.flat, Y.flat, A, V
984	else:
985	return A, V
986
987
988	def extrapolate_first_order(self):
989	"""Extrapolate conserved quantities from centroid to
990	vertices for each volume using
991	first order scheme.
992	"""
993
994	qc = self.centroid_values
995	qv = self.vertex_values
996
997	for i in range(3):
998	qv[:,i] = qc
999
1000
1001	def get_integral(self):
1002	"""Compute the integral of quantity across entire domain
1003	"""
1004	integral = 0
1005	for k in range(self.domain.number_of_elements):
1006	area = self.domain.areas[k]
1007	qc = self.centroid_values[k]
1008	integral += qc*area
1009
1010	return integral
1011
1012
1013
1014
1015	class Conserved_quantity(Quantity):
1016	"""Class conserved quantity adds to Quantity:
1017
1018	boundary values, storage and method for updating, and
1019	methods for (second order) extrapolation from centroid to vertices inluding
1020	gradients and limiters
1021	"""
1022
1023	def __init__(self, domain, vertex_values=None):
1024	Quantity.__init__(self, domain, vertex_values)
1025
1026
1027	#Allocate space for boundary values
1028	L = len(domain.boundary)
1029	self.boundary_values = zeros(L, Float)
1030
1031	#Allocate space for updates of conserved quantities by
1032	#flux calculations and forcing functions
1033
1034	N = domain.number_of_elements
1035	self.explicit_update = zeros(N, Float )
1036	self.semi_implicit_update = zeros(N, Float )
1037
1038
1039	def update(self, timestep):
1040	#Call correct module function
1041	#(either from this module or C-extension)
1042	return update(self, timestep)
1043
1044
1045	def compute_gradients(self):
1046	#Call correct module function
1047	#(either from this module or C-extension)
1048	return compute_gradients(self)
1049
1050
1051	def limit(self):
1052	#Call correct module function
1053	#(either from this module or C-extension)
1054	limit(self)
1055
1056
1057	def extrapolate_second_order(self):
1058	#Call correct module function
1059	#(either from this module or C-extension)
1060	extrapolate_second_order(self)
1061
1062
1063	def update(quantity, timestep):
1064	"""Update centroid values based on values stored in
1065	explicit_update and semi_implicit_update as well as given timestep
1066
1067	Function implementing forcing terms must take on argument
1068	which is the domain and they must update either explicit
1069	or implicit updates, e,g,:
1070
1071	def gravity(domain):
1072	....
1073	domain.quantities['xmomentum'].explicit_update = ...
1074	domain.quantities['ymomentum'].explicit_update = ...
1075
1076
1077
1078	Explicit terms must have the form
1079
1080	G(q, t)
1081
1082	and explicit scheme is
1083
1084	q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
1085
1086
1087	Semi implicit forcing terms are assumed to have the form
1088
1089	G(q, t) = H(q, t) q
1090
1091	and the semi implicit scheme will then be
1092
1093	q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
1094
1095
1096	"""
1097
1098	#from Numeric import sum, equal, ones, Float
1099
1100	N = quantity.centroid_values.shape[0]
1101
1102
1103	#Divide H by conserved quantity to obtain G (see docstring above)
1104
1105
1106	for k in range(N):
1107	x = quantity.centroid_values[k]
1108	if x == 0.0:
1109	#FIXME: Is this right
1110	quantity.semi_implicit_update[k] = 0.0
1111	else:
1112	quantity.semi_implicit_update[k] /= x
1113
1114	#Explicit updates
1115	quantity.centroid_values += timestep*quantity.explicit_update
1116
1117	#Semi implicit updates
1118	denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
1119
1120	if sum(equal(denominator, 0.0)) > 0.0:
1121	msg = 'Zero division in semi implicit update. Call Stephen :-)'
1122	raise msg
1123	else:
1124	#Update conserved_quantities from semi implicit updates
1125	quantity.centroid_values /= denominator
1126
1127
1128	def interpolate_from_vertices_to_edges(quantity):
1129	"""Compute edge values from vertex values using linear interpolation
1130	"""
1131
1132	for k in range(quantity.vertex_values.shape[0]):
1133	q0 = quantity.vertex_values[k, 0]
1134	q1 = quantity.vertex_values[k, 1]
1135	q2 = quantity.vertex_values[k, 2]
1136
1137	quantity.edge_values[k, 0] = 0.5*(q1+q2)
1138	quantity.edge_values[k, 1] = 0.5*(q0+q2)
1139	quantity.edge_values[k, 2] = 0.5*(q0+q1)
1140
1141
1142
1143	def extrapolate_second_order(quantity):
1144	"""Extrapolate conserved quantities from centroid to
1145	vertices for each volume using
1146	second order scheme.
1147	"""
1148
1149	a, b = quantity.compute_gradients()
1150
1151	X = quantity.domain.get_vertex_coordinates()
1152	qc = quantity.centroid_values
1153	qv = quantity.vertex_values
1154
1155	#Check each triangle
1156	for k in range(quantity.domain.number_of_elements):
1157	#Centroid coordinates
1158	x, y = quantity.domain.centroid_coordinates[k]
1159
1160	#vertex coordinates
1161	x0, y0, x1, y1, x2, y2 = X[k,:]
1162
1163	#Extrapolate
1164	qv[k,0] = qc[k] + a[k](x0-x) + b[k](y0-y)
1165	qv[k,1] = qc[k] + a[k](x1-x) + b[k](y1-y)
1166	qv[k,2] = qc[k] + a[k](x2-x) + b[k](y2-y)
1167
1168
1169	def compute_gradients(quantity):
1170	"""Compute gradients of triangle surfaces defined by centroids of
1171	neighbouring volumes.
1172	If one edge is on the boundary, use own centroid as neighbour centroid.
1173	If two or more are on the boundary, fall back to first order scheme.
1174	"""
1175
1176	from utilitites.numerical_tools import gradient
1177
1178	centroid_coordinates = quantity.domain.centroid_coordinates
1179	surrogate_neighbours = quantity.domain.surrogate_neighbours
1180	centroid_values = quantity.centroid_values
1181	number_of_boundaries = quantity.domain.number_of_boundaries
1182
1183	N = centroid_values.shape[0]
1184
1185	a = zeros(N, Float)
1186	b = zeros(N, Float)
1187
1188	for k in range(N):
1189	if number_of_boundaries[k] < 2:
1190	#Two or three true neighbours
1191
1192	#Get indices of neighbours (or self when used as surrogate)
1193	k0, k1, k2 = surrogate_neighbours[k,:]
1194
1195	#Get data
1196	q0 = centroid_values[k0]
1197	q1 = centroid_values[k1]
1198	q2 = centroid_values[k2]
1199
1200	x0, y0 = centroid_coordinates[k0] #V0 centroid
1201	x1, y1 = centroid_coordinates[k1] #V1 centroid
1202	x2, y2 = centroid_coordinates[k2] #V2 centroid
1203
1204	#Gradient
1205	a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
1206
1207	elif number_of_boundaries[k] == 2:
1208	#One true neighbour
1209
1210	#Get index of the one neighbour
1211	for k0 in surrogate_neighbours[k,:]:
1212	if k0 != k: break
1213	assert k0 != k
1214
1215	k1 = k #self
1216
1217	#Get data
1218	q0 = centroid_values[k0]
1219	q1 = centroid_values[k1]
1220
1221	x0, y0 = centroid_coordinates[k0] #V0 centroid
1222	x1, y1 = centroid_coordinates[k1] #V1 centroid
1223
1224	#Gradient
1225	a[k], b[k] = gradient2(x0, y0, x1, y1, q0, q1)
1226	else:
1227	#No true neighbours -
1228	#Fall back to first order scheme
1229	pass
1230
1231
1232	return a, b
1233
1234
1235
1236	def limit(quantity):
1237	"""Limit slopes for each volume to eliminate artificial variance
1238	introduced by e.g. second order extrapolator
1239
1240	This is an unsophisticated limiter as it does not take into
1241	account dependencies among quantities.
1242
1243	precondition:
1244	vertex values are estimated from gradient
1245	postcondition:
1246	vertex values are updated
1247	"""
1248
1249	N = quantity.domain.number_of_elements
1250
1251	beta_w = quantity.domain.beta_w
1252
1253	qc = quantity.centroid_values
1254	qv = quantity.vertex_values
1255
1256	#Find min and max of this and neighbour's centroid values
1257	qmax = zeros(qc.shape, Float)
1258	qmin = zeros(qc.shape, Float)
1259
1260	for k in range(N):
1261	qmax[k] = qmin[k] = qc[k]
1262	for i in range(3):
1263	n = quantity.domain.neighbours[k,i]
1264	if n >= 0:
1265	qn = qc[n] #Neighbour's centroid value
1266
1267	qmin[k] = min(qmin[k], qn)
1268	qmax[k] = max(qmax[k], qn)
1269
1270
1271	#Diffences between centroids and maxima/minima
1272	dqmax = qmax - qc
1273	dqmin = qmin - qc
1274
1275	#Deltas between vertex and centroid values
1276	dq = zeros(qv.shape, Float)
1277	for i in range(3):
1278	dq[:,i] = qv[:,i] - qc
1279
1280	#Phi limiter
1281	for k in range(N):
1282
1283	#Find the gradient limiter (phi) across vertices
1284	phi = 1.0
1285	for i in range(3):
1286	r = 1.0
1287	if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
1288	if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
1289
1290	phi = min( min(r*beta_w, 1), phi )
1291
1292	#Then update using phi limiter
1293	for i in range(3):
1294	qv[k,i] = qc[k] + phi*dq[k,i]
1295
1296
1297
1298	from utilities import compile
1299	if compile.can_use_C_extension('quantity_ext.c'):
1300	#Replace python version with c implementations
1301
1302	from quantity_ext import limit, compute_gradients,\
1303	extrapolate_second_order, interpolate_from_vertices_to_edges, update

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

Download in other formats: