Context Navigation

source: anuga_core/source/anuga/fit_interpolate/fit.py @ 6537

Last change on this file since 6537 was 6537, checked in by ole, 15 years ago
Formatting
File size: 26.6 KB

Line
1	"""Least squares fitting.
2
3	Implements a penalised least-squares fit.
4	putting point data onto the mesh.
5
6	The penalty term (or smoothing term) is controlled by the smoothing
7	parameter alpha.
8	With a value of alpha=0, the fit function will attempt
9	to interpolate as closely as possible in the least-squares sense.
10	With values alpha > 0, a certain amount of smoothing will be applied.
11	A positive alpha is essential in cases where there are too few
12	data points.
13	A negative alpha is not allowed.
14	A typical value of alpha is 1.0e-6
15
16
17	Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
18	Geoscience Australia, 2004.
19
20	TO DO
21	* test geo_ref, geo_spatial
22
23	IDEAS
24	* (DSG-) Change the interface of fit, so a domain object can
25	be passed in. (I don't know if this is feasible). If could
26	save time/memory.
27	"""
28	import types
29
30	from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
31	from anuga.caching import cache
32	from anuga.geospatial_data.geospatial_data import Geospatial_data, \
33	ensure_absolute
34	from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
35	from anuga.utilities.sparse import Sparse, Sparse_CSR
36	from anuga.utilities.polygon import inside_polygon, is_inside_polygon
37	from anuga.fit_interpolate.search_functions import search_tree_of_vertices
38
39	from anuga.utilities.cg_solve import conjugate_gradient
40	from anuga.utilities.numerical_tools import ensure_numeric, gradient
41	from anuga.config import default_smoothing_parameter as DEFAULT_ALPHA
42
43	import exceptions
44	class TooFewPointsError(exceptions.Exception): pass
45	class VertsWithNoTrianglesError(exceptions.Exception): pass
46
47	import Numeric as num
48
49
50	class Fit(FitInterpolate):
51
52	def __init__(self,
53	vertex_coordinates=None,
54	triangles=None,
55	mesh=None,
56	mesh_origin=None,
57	alpha = None,
58	verbose=False,
59	max_vertices_per_cell=None):
60
61
62	"""
63	Fit data at points to the vertices of a mesh.
64
65	Inputs:
66
67	vertex_coordinates: List of coordinate pairs [xi, eta] of
68	points constituting a mesh (or an m x 2 Numeric array or
69	a geospatial object)
70	Points may appear multiple times
71	(e.g. if vertices have discontinuities)
72
73	triangles: List of 3-tuples (or a Numeric array) of
74	integers representing indices of all vertices in the mesh.
75
76	mesh_origin: A geo_reference object or 3-tuples consisting of
77	UTM zone, easting and northing.
78	If specified vertex coordinates are assumed to be
79	relative to their respective origins.
80
81	max_vertices_per_cell: Number of vertices in a quad tree cell
82	at which the cell is split into 4.
83
84	Note: Don't supply a vertex coords as a geospatial object and
85	a mesh origin, since geospatial has its own mesh origin.
86
87
88	Usage,
89	To use this in a blocking way, call build_fit_subset, with z info,
90	and then fit, with no point coord, z info.
91
92	"""
93	# Initialise variabels
94	if alpha is None:
95	self.alpha = DEFAULT_ALPHA
96	else:
97	self.alpha = alpha
98
99	FitInterpolate.__init__(self,
100	vertex_coordinates,
101	triangles,
102	mesh,
103	mesh_origin,
104	verbose,
105	max_vertices_per_cell)
106
107	m = self.mesh.number_of_nodes # Nbr of basis functions (vertices)
108
109	self.AtA = None
110	self.Atz = None
111
112	self.point_count = 0
113	if self.alpha <> 0:
114	if verbose: print 'Building smoothing matrix'
115	self._build_smoothing_matrix_D()
116
117	bd_poly = self.mesh.get_boundary_polygon()
118	self.mesh_boundary_polygon = ensure_numeric(bd_poly)
119
120	def _build_coefficient_matrix_B(self,
121	verbose = False):
122	"""
123	Build final coefficient matrix
124
125	Precon
126	If alpha is not zero, matrix D has been built
127	Matrix Ata has been built
128	"""
129
130	if self.alpha <> 0:
131	#if verbose: print 'Building smoothing matrix'
132	#self._build_smoothing_matrix_D()
133	self.B = self.AtA + self.alpha*self.D
134	else:
135	self.B = self.AtA
136
137	# Convert self.B matrix to CSR format for faster matrix vector
138	self.B = Sparse_CSR(self.B)
139
140	def _build_smoothing_matrix_D(self):
141	"""Build m x m smoothing matrix, where
142	m is the number of basis functions phi_k (one per vertex)
143
144	The smoothing matrix is defined as
145
146	D = D1 + D2
147
148	where
149
150	[D1]_{k,l} = \int_\Omega
151	\frac{\partial \phi_k}{\partial x}
152	\frac{\partial \phi_l}{\partial x}\,
153	dx dy
154
155	[D2]_{k,l} = \int_\Omega
156	\frac{\partial \phi_k}{\partial y}
157	\frac{\partial \phi_l}{\partial y}\,
158	dx dy
159
160
161	The derivatives \frac{\partial \phi_k}{\partial x},
162	\frac{\partial \phi_k}{\partial x} for a particular triangle
163	are obtained by computing the gradient a_k, b_k for basis function k
164	"""
165
166	# FIXME: algorithm might be optimised by computing local 9x9
167	# "element stiffness matrices:
168
169	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
170
171	self.D = Sparse(m,m)
172
173	# For each triangle compute contributions to D = D1+D2
174	for i in range(len(self.mesh)):
175
176	# Get area
177	area = self.mesh.areas[i]
178
179	# Get global vertex indices
180	v0 = self.mesh.triangles[i,0]
181	v1 = self.mesh.triangles[i,1]
182	v2 = self.mesh.triangles[i,2]
183
184	# Get the three vertex_points
185	xi0 = self.mesh.get_vertex_coordinate(i, 0)
186	xi1 = self.mesh.get_vertex_coordinate(i, 1)
187	xi2 = self.mesh.get_vertex_coordinate(i, 2)
188
189	# Compute gradients for each vertex
190	a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
191	1, 0, 0)
192
193	a1, b1 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
194	0, 1, 0)
195
196	a2, b2 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
197	0, 0, 1)
198
199	# Compute diagonal contributions
200	self.D[v0,v0] += (a0a0 + b0b0)*area
201	self.D[v1,v1] += (a1a1 + b1b1)*area
202	self.D[v2,v2] += (a2a2 + b2b2)*area
203
204	# Compute contributions for basis functions sharing edges
205	e01 = (a0a1 + b0b1)*area
206	self.D[v0,v1] += e01
207	self.D[v1,v0] += e01
208
209	e12 = (a1a2 + b1b2)*area
210	self.D[v1,v2] += e12
211	self.D[v2,v1] += e12
212
213	e20 = (a2a0 + b2b0)*area
214	self.D[v2,v0] += e20
215	self.D[v0,v2] += e20
216
217	def get_D(self):
218	return self.D.todense()
219
220
221
222	def _build_matrix_AtA_Atz(self,
223	point_coordinates,
224	z,
225	verbose = False):
226	"""Build:
227	AtA m x m interpolation matrix, and,
228	Atz m x a interpolation matrix where,
229	m is the number of basis functions phi_k (one per vertex)
230	a is the number of data attributes
231
232	This algorithm uses a quad tree data structure for fast binning of
233	data points.
234
235	If Ata is None, the matrices AtA and Atz are created.
236
237	This function can be called again and again, with sub-sets of
238	the point coordinates. Call fit to get the results.
239
240	Preconditions
241	z and points are numeric
242	Point_coordindates and mesh vertices have the same origin.
243
244	The number of attributes of the data points does not change
245	"""
246
247	# Build n x m interpolation matrix
248	if self.AtA == None:
249	# AtA and Atz need to be initialised.
250	m = self.mesh.number_of_nodes
251	if len(z.shape) > 1:
252	att_num = z.shape[1]
253	self.Atz = num.zeros((m,att_num), num.Float)
254	else:
255	att_num = 1
256	self.Atz = num.zeros((m,), num.Float)
257	assert z.shape[0] == point_coordinates.shape[0]
258
259	AtA = Sparse(m,m)
260	# The memory damage has been done by now.
261	else:
262	AtA = self.AtA # Did this for speed, did ~nothing
263	self.point_count += point_coordinates.shape[0]
264
265
266	inside_indices = inside_polygon(point_coordinates,
267	self.mesh_boundary_polygon,
268	closed=True,
269	verbose=False) # Suppress output
270
271	n = len(inside_indices)
272
273	# Compute matrix elements for points inside the mesh
274	triangles = self.mesh.triangles # Shorthand
275	for d, i in enumerate(inside_indices):
276	# For each data_coordinate point
277	# if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
278	x = point_coordinates[i]
279
280	element_found, sigma0, sigma1, sigma2, k = \
281	search_tree_of_vertices(self.root, self.mesh, x)
282
283	if element_found is True:
284	j0 = triangles[k,0] # Global vertex id for sigma0
285	j1 = triangles[k,1] # Global vertex id for sigma1
286	j2 = triangles[k,2] # Global vertex id for sigma2
287
288	sigmas = {j0:sigma0, j1:sigma1, j2:sigma2}
289	js = [j0,j1,j2]
290
291	for j in js:
292	self.Atz[j] += sigmas[j]*z[i]
293	#print "self.Atz building", self.Atz
294	#print "self.Atz[j]", self.Atz[j]
295	#print " sigmas[j]", sigmas[j]
296	#print "z[i]",z[i]
297	#print "result", sigmas[j]*z[i]
298
299	for k in js:
300	AtA[j,k] += sigmas[j]*sigmas[k]
301	else:
302	flag = is_inside_polygon(x,
303	self.mesh_boundary_polygon,
304	closed=True,
305	verbose=False) # Suppress output
306	msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
307	assert flag is True, msg
308
309	# FIXME(Ole): This is the message referred to in ticket:314
310	minx = min(self.mesh_boundary_polygon[:,0])
311	maxx = max(self.mesh_boundary_polygon[:,0])
312	miny = min(self.mesh_boundary_polygon[:,1])
313	maxy = max(self.mesh_boundary_polygon[:,1])
314	msg = 'Could not find triangle for point %s. ' % str(x)
315	msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
316	% (minx, maxx, miny, maxy)
317	#msg += '\nBoundary polygon = %s' %str(self.mesh_boundary_polygon)
318	raise RuntimeError, msg
319
320	self.AtA = AtA
321
322
323	def fit(self, point_coordinates_or_filename=None, z=None,
324	verbose=False,
325	point_origin=None,
326	attribute_name=None,
327	max_read_lines=500):
328	"""Fit a smooth surface to given 1d array of data points z.
329
330	The smooth surface is computed at each vertex in the underlying
331	mesh using the formula given in the module doc string.
332
333	Inputs:
334	point_coordinates: The co-ordinates of the data points.
335	List of coordinate pairs [x, y] of
336	data points or an nx2 Numeric array or a Geospatial_data object
337	or points file filename
338	z: Single 1d vector or array of data at the point_coordinates.
339
340	"""
341	# Use blocking to load in the point info
342	if type(point_coordinates_or_filename) == types.StringType:
343	msg = "Don't set a point origin when reading from a file"
344	assert point_origin is None, msg
345	filename = point_coordinates_or_filename
346	G_data = Geospatial_data(filename,
347	max_read_lines=max_read_lines,
348	load_file_now=False,
349	verbose=verbose)
350
351	for i, geo_block in enumerate(G_data):
352	if verbose is True and 0 == i%200:
353	# The time this will take
354	# is dependant on the # of Triangles
355
356	print 'Processing Block %d' %i
357	# FIXME (Ole): It would be good to say how many blocks
358	# there are here. But this is no longer necessary
359	# for pts files as they are reported in geospatial_data
360	# I suggest deleting this verbose output and make
361	# Geospatial_data more informative for txt files.
362	#
363	# I still think so (12/12/7, Ole).
364
365
366
367	# Build the array
368
369	points = geo_block.get_data_points(absolute=True)
370	z = geo_block.get_attributes(attribute_name=attribute_name)
371	self.build_fit_subset(points, z, verbose=verbose)
372
373	# FIXME(Ole): I thought this test would make sense here
374	# See test_fitting_example_that_crashed_2 in test_shallow_water_domain.py
375	# Committed 11 March 2009
376	msg = 'Matrix AtA was not built'
377	assert self.AtA is not None, msg
378
379	#print 'Matrix was built OK'
380
381
382	point_coordinates = None
383	else:
384	point_coordinates = point_coordinates_or_filename
385
386	if point_coordinates is None:
387	if verbose: print 'Warning: no data points in fit'
388	msg = 'No interpolation matrix.'
389	assert self.AtA is not None, msg
390	assert self.Atz is not None
391
392	# FIXME (DSG) - do a message
393	else:
394	point_coordinates = ensure_absolute(point_coordinates,
395	geo_reference=point_origin)
396	# if isinstance(point_coordinates,Geospatial_data) and z is None:
397	# z will come from the geo-ref
398	self.build_fit_subset(point_coordinates, z, verbose)
399
400	# Check sanity
401	m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
402	n = self.point_count
403	if n<m and self.alpha == 0.0:
404	msg = 'ERROR (least_squares): Too few data points\n'
405	msg += 'There are only %d data points and alpha == 0. ' %n
406	msg += 'Need at least %d\n' %m
407	msg += 'Alternatively, set smoothing parameter alpha to a small '
408	msg += 'positive value,\ne.g. 1.0e-3.'
409	raise TooFewPointsError(msg)
410
411	self._build_coefficient_matrix_B(verbose)
412	loners = self.mesh.get_lone_vertices()
413	# FIXME - make this as error message.
414	# test with
415	# Not_yet_test_smooth_att_to_mesh_with_excess_verts.
416	if len(loners)>0:
417	msg = 'WARNING: (least_squares): \nVertices with no triangles\n'
418	msg += 'All vertices should be part of a triangle.\n'
419	msg += 'In the future this will be inforced.\n'
420	msg += 'The following vertices are not part of a triangle;\n'
421	msg += str(loners)
422	print msg
423	#raise VertsWithNoTrianglesError(msg)
424
425
426	return conjugate_gradient(self.B, self.Atz, self.Atz,
427	imax=2*len(self.Atz) )
428
429
430	def build_fit_subset(self, point_coordinates, z=None, attribute_name=None,
431	verbose=False):
432	"""Fit a smooth surface to given 1d array of data points z.
433
434	The smooth surface is computed at each vertex in the underlying
435	mesh using the formula given in the module doc string.
436
437	Inputs:
438	point_coordinates: The co-ordinates of the data points.
439	List of coordinate pairs [x, y] of
440	data points or an nx2 Numeric array or a Geospatial_data object
441	z: Single 1d vector or array of data at the point_coordinates.
442	attribute_name: Used to get the z values from the
443	geospatial object if no attribute_name is specified,
444	it's a bit of a lucky dip as to what attributes you get.
445	If there is only one attribute it will be that one.
446
447	"""
448
449	# FIXME(DSG-DSG): Check that the vert and point coords
450	# have the same zone.
451	if isinstance(point_coordinates,Geospatial_data):
452	point_coordinates = point_coordinates.get_data_points( \
453	absolute = True)
454
455	# Convert input to Numeric arrays
456	if z is not None:
457	z = ensure_numeric(z, num.Float)
458	else:
459	msg = 'z not specified'
460	assert isinstance(point_coordinates,Geospatial_data), msg
461	z = point_coordinates.get_attributes(attribute_name)
462
463	point_coordinates = ensure_numeric(point_coordinates, num.Float)
464	self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
465
466
467	############################################################################
468
469	def fit_to_mesh(point_coordinates, # this can also be a points file name
470	vertex_coordinates=None,
471	triangles=None,
472	mesh=None,
473	point_attributes=None,
474	alpha=DEFAULT_ALPHA,
475	verbose=False,
476	acceptable_overshoot=1.01,
477	# FIXME: Move to config - this value is assumed in caching test
478	# FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.
479	mesh_origin=None,
480	data_origin=None,
481	max_read_lines=None,
482	attribute_name=None,
483	use_cache=False):
484	"""Wrapper around internal function _fit_to_mesh for use with caching.
485
486	"""
487
488	args = (point_coordinates, )
489	kwargs = {'vertex_coordinates': vertex_coordinates,
490	'triangles': triangles,
491	'mesh': mesh,
492	'point_attributes': point_attributes,
493	'alpha': alpha,
494	'verbose': verbose,
495	'acceptable_overshoot': acceptable_overshoot,
496	'mesh_origin': mesh_origin,
497	'data_origin': data_origin,
498	'max_read_lines': max_read_lines,
499	'attribute_name': attribute_name
500	}
501
502	if use_cache is True:
503	if isinstance(point_coordinates, basestring):
504	# We assume that point_coordinates is the name of a .csv/.txt
505	# file which must be passed onto caching as a dependency
506	# (in case it has changed on disk)
507	dep = [point_coordinates]
508	else:
509	dep = None
510
511
512	#from caching import myhash
513	#import copy
514	#print args
515	#print kwargs
516	#print 'hashing:'
517	#print 'args', myhash( (args, kwargs) )
518	#print 'again', myhash( copy.deepcopy( (args, kwargs)) )
519
520	#print 'mesh hash', myhash( kwargs['mesh'] )
521
522	#print '-------------------------'
523	#print 'vertices hash', myhash( kwargs['mesh'].nodes )
524	#print 'triangles hash', myhash( kwargs['mesh'].triangles )
525	#print '-------------------------'
526
527	#for key in mesh.__dict__:
528	# print key, myhash(mesh.__dict__[key])
529
530	#for key in mesh.quantities.keys():
531	# print key, myhash(mesh.quantities[key])
532
533	#import sys; sys.exit()
534
535	return cache(_fit_to_mesh,
536	args, kwargs,
537	verbose=verbose,
538	compression=False,
539	dependencies=dep)
540	else:
541	return apply(_fit_to_mesh,
542	args, kwargs)
543
544	def _fit_to_mesh(point_coordinates, # this can also be a points file name
545	vertex_coordinates=None,
546	triangles=None,
547	mesh=None,
548	point_attributes=None,
549	alpha=DEFAULT_ALPHA,
550	verbose=False,
551	acceptable_overshoot=1.01,
552	mesh_origin=None,
553	data_origin=None,
554	max_read_lines=None,
555	attribute_name=None):
556	"""
557	Fit a smooth surface to a triangulation,
558	given data points with attributes.
559
560
561	Inputs:
562	vertex_coordinates: List of coordinate pairs [xi, eta] of
563	points constituting a mesh (or an m x 2 Numeric array or
564	a geospatial object)
565	Points may appear multiple times
566	(e.g. if vertices have discontinuities)
567
568	triangles: List of 3-tuples (or a Numeric array) of
569	integers representing indices of all vertices in the mesh.
570
571	point_coordinates: List of coordinate pairs [x, y] of data points
572	(or an nx2 Numeric array). This can also be a .csv/.txt/.pts
573	file name.
574
575	alpha: Smoothing parameter.
576
577	acceptable overshoot: NOT IMPLEMENTED
578	controls the allowed factor by which
579	fitted values
580	may exceed the value of input data. The lower limit is defined
581	as min(z) - acceptable_overshoot*delta z and upper limit
582	as max(z) + acceptable_overshoot*delta z
583
584
585	mesh_origin: A geo_reference object or 3-tuples consisting of
586	UTM zone, easting and northing.
587	If specified vertex coordinates are assumed to be
588	relative to their respective origins.
589
590
591	point_attributes: Vector or array of data at the
592	point_coordinates.
593
594	"""
595
596	if mesh is None:
597	# FIXME(DSG): Throw errors if triangles or vertex_coordinates
598	# are None
599
600	#Convert input to Numeric arrays
601	triangles = ensure_numeric(triangles, num.Int)
602	vertex_coordinates = ensure_absolute(vertex_coordinates,
603	geo_reference = mesh_origin)
604
605	if verbose: print 'FitInterpolate: Building mesh'
606	mesh = Mesh(vertex_coordinates, triangles)
607	mesh.check_integrity()
608
609
610	interp = Fit(mesh=mesh,
611	verbose=verbose,
612	alpha=alpha)
613
614	vertex_attributes = interp.fit(point_coordinates,
615	point_attributes,
616	point_origin=data_origin,
617	max_read_lines=max_read_lines,
618	attribute_name=attribute_name,
619	verbose=verbose)
620
621
622	# Add the value checking stuff that's in least squares.
623	# Maybe this stuff should get pushed down into Fit.
624	# at least be a method of Fit.
625	# Or intigrate it into the fit method, saving teh max and min's
626	# as att's.
627
628	return vertex_attributes
629
630
631	#def _fit(args, *kwargs):
632	# """Private function for use with caching. Reason is that classes
633	# may change their byte code between runs which is annoying.
634	# """
635	#
636	# return Fit(args, *kwargs)
637
638
639	def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
640	alpha=DEFAULT_ALPHA, verbose= False,
641	expand_search = False,
642	precrop = False,
643	display_errors = True):
644	"""
645	Given a mesh file (tsh) and a point attribute file, fit
646	point attributes to the mesh and write a mesh file with the
647	results.
648
649	Note: the points file needs titles. If you want anuga to use the tsh file,
650	make sure the title is elevation.
651
652	NOTE: Throws IOErrors, for a variety of file problems.
653
654	"""
655
656	from load_mesh.loadASCII import import_mesh_file, \
657	export_mesh_file, concatinate_attributelist
658
659
660	try:
661	mesh_dict = import_mesh_file(mesh_file)
662	except IOError,e:
663	if display_errors:
664	print "Could not load bad file. ", e
665	raise IOError #Could not load bad mesh file.
666
667	vertex_coordinates = mesh_dict['vertices']
668	triangles = mesh_dict['triangles']
669	if type(mesh_dict['vertex_attributes']) == num.ArrayType:
670	old_point_attributes = mesh_dict['vertex_attributes'].tolist()
671	else:
672	old_point_attributes = mesh_dict['vertex_attributes']
673
674	if type(mesh_dict['vertex_attribute_titles']) == num.ArrayType:
675	old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
676	else:
677	old_title_list = mesh_dict['vertex_attribute_titles']
678
679	if verbose: print 'tsh file %s loaded' %mesh_file
680
681	# load in the points file
682	try:
683	geo = Geospatial_data(point_file, verbose=verbose)
684	except IOError,e:
685	if display_errors:
686	print "Could not load bad file. ", e
687	raise IOError #Re-raise exception
688
689	point_coordinates = geo.get_data_points(absolute=True)
690	title_list,point_attributes = concatinate_attributelist( \
691	geo.get_all_attributes())
692
693	if mesh_dict.has_key('geo_reference') and \
694	not mesh_dict['geo_reference'] is None:
695	mesh_origin = mesh_dict['geo_reference'].get_origin()
696	else:
697	mesh_origin = None
698
699	if verbose: print "points file loaded"
700	if verbose: print "fitting to mesh"
701	f = fit_to_mesh(point_coordinates,
702	vertex_coordinates,
703	triangles,
704	None,
705	point_attributes,
706	alpha = alpha,
707	verbose = verbose,
708	data_origin = None,
709	mesh_origin = mesh_origin)
710	if verbose: print "finished fitting to mesh"
711
712	# convert array to list of lists
713	new_point_attributes = f.tolist()
714	#FIXME have this overwrite attributes with the same title - DSG
715	#Put the newer attributes last
716	if old_title_list <> []:
717	old_title_list.extend(title_list)
718	#FIXME can this be done a faster way? - DSG
719	for i in range(len(old_point_attributes)):
720	old_point_attributes[i].extend(new_point_attributes[i])
721	mesh_dict['vertex_attributes'] = old_point_attributes
722	mesh_dict['vertex_attribute_titles'] = old_title_list
723	else:
724	mesh_dict['vertex_attributes'] = new_point_attributes
725	mesh_dict['vertex_attribute_titles'] = title_list
726
727	if verbose: print "exporting to file ", mesh_output_file
728
729	try:
730	export_mesh_file(mesh_output_file, mesh_dict)
731	except IOError,e:
732	if display_errors:
733	print "Could not write file. ", e
734	raise IOError

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

Download in other formats: