Context Navigation

source: trunk/anuga_core/source/anuga/structures/riverwall.py @ 9234

Last change on this file since 9234 was 9234, checked in by steve, 11 years ago
updating version.py
File size: 33.3 KB

Line
1	import os
2	from anuga_parallel import barrier, numprocs, myid
3	import numpy
4
5	class RiverWall:
6	"""Define the elevation of 'riverwalls'.
7
8	These are located along each cell edge, and can have an elevation different
9	from the bed elevation.
10
11	For the DE algorithms, they are used in computing the 'z_half' value [if they
12	are greater than either edge elevation].
13
14	In addition, the DE fluxes at riverwalls are adjusted to agree with a weir relation,
15	so long as the riverwall is not too deeply submerged.
16
17	As of 22/04/2014, they are only implemented for DE algorithms [the shallow water
18	component would be more difficult to implement with other algorithms]
19
20	How the flux over the riverwall is computed:
21	# We have a riverwall, which is treated as a weir.
22	#
23	# Say the headwater-head is the upstream water depth above the weir elevation (min 0), and
24	# the tailwater head is the downstream water depth above the weir elevation (min 0).
25	# By definition headwater-head > tailwater-head.
26	#
27	# Let s = (headwater head) / (tailwater head), and h = (tailwater head)/(weir height)
28	# where 'weir height' is the elevation of the weir crest above the minimum
29	# of the left/right bed elevations
30	#
31	# Denote ID as the 'ideal' weir flow, computed from the hydraulic
32	# formula (including a submergence correction factor from Villemonte,
33	# 1947),
34	# Q1 = 2/3headwater_headsqrt(g2/3headwater_head)*Qfactor
35	# Q2 = 2/3tailwater_headsqrt(g2/3headwater_head)*Qfactor
36	# ID = Q1(1-Q2/Q1)*0.385
37	#
38	# Denote SW as the 'shallow-water' weir flux, computed from the approximate
39	# reimann solver, where the mid-edge-elevation is the weir crest elevation.
40	# This makes clear sense for DE0 and DE1. Cell centroid stage/height/bed_elevation
41	# are used in the flux computation
42	#
43	# Then the flux over the weir is computed from:
44	# w1 = min( max(s-s1, 0.)/(s2-s1), 1.0) # Factor describing relative submergence
45	# w1' = min( max(h-h1,0.)/(h2-h1), 1.0) # Factor describing absolute submergence
46	# flux_over_weir = (w1SW + (1-w1)ID)( 1-w1') + (w1')SW
47	#
48	# The key idea is that if s<s1, h<h1, then the ideal weir solution is
49	# used. Otherwise, we linearly blend with the SW solution,
50	# and the SW solution is used completely if s>s2 or h>h2
51
52	"""
53
54	def __init__(self, domain):
55	"""Riverwall data structure
56
57	Allows reasonably efficient storage of riverwall elevations and hydraulic
58	properties
59
60	NOTE: In domain.edge_flux_type, riverwalls correspond to the value 1
61
62	RiverWall variables are initialised to dummy values, because even
63	if no riverwalls are used, some values have to be passed to the flux computation
64
65	Hydraulic parameters are
66	Qfactor -- Multiplicative factor for ideal weir relation (calibration coef)
67	s1 -- Submergence ratio at which we start blending with the shallow water solution (<s2)
68	s2 -- Submergence ratio at which we entirely use the shallow water solution (>s1)
69	h1 -- Tailwater head / weir height at which we start blending with the shallow water solution (<h2)
70	h2 -- Tailwater head / weir height at which we entirely use the shallow water solution (>h1)
71
72	# Default riverwall hydraulic parameters
73	default_riverwallPar={'Qfactor':1.0,
74	's1': 0.9,
75	's2': 0.95,
76	'h1': 1.0,
77	'h2': 1.5
78	}
79
80	Other variables are:
81
82	riverwall_elevation -- Variable to hold riverwall elevations.
83	len = number of riverwall edges in the domain.
84	First entry corresponds to first riverwall edge in domain.edge_coordinates,
85	second entry to second riverwall edge in domain.edge_coordinates, etc
86
87	hydraulic_properties_rowIndex -- Variable to hold the row-index of the hydraulic properties table
88	len = number of riverwall edges
89	in the domain, ordered like riverwall_elevation
90
91	riverwall_edges -- Holds indices of edges in domain which are riverwalls, ordered like riverwall_elevation
92
93	names -- list with the names of the riverwalls
94	len = number of riverwalls which cover edges in the domain
95
96	hydraulic_variable_names -- Variables to hold the names of variables in columns of the hydraulic
97	properties table. THE ORDER IS IMPORTANT -- C code ASSUMES IT
98
99	ncol_hydraulic_properties -- Number of columns in the hydraulic properties table
100	[ = len(hydraulic_variable_names) ]
101
102	hydraulic_properties -- Array with the hydraulic parameters for each riverwall.
103	number of rows = number of riverwalls which cover edges in the domain
104	number of cols = number of hydraulic variable names
105
106
107	input_riverwall_geo, input_riverwall_par -- holds input information
108
109	"""
110	self.domain=domain
111
112	default_float=-9.0e+20
113	default_int=-9e+20
114	self.riverwall_elevation=numpy.array([default_float])
115
116	self.hydraulic_properties_rowIndex=numpy.array([default_int])
117
118	self.names=[ ]
119
120	# Default riverwall hydraulic parameters
121	self.default_riverwallPar={'Qfactor':1.0,
122	's1': 0.9,
123	's2': 0.95,
124	'h1': 1.0,
125	'h2': 1.5
126	}
127
128	# DO NOT CHANGE THE ORDER OF hydraulic_variable_names
129	# It needs to match hard-coded assumptions in C [compute_fluxes_central]
130	# If you add a variable, append it to the end of hydraulic_variable_names
131	self.hydraulic_variable_names=('Qfactor', 's1', 's2', 'h1', 'h2')
132
133	self.ncol_hydraulic_properties=len(self.hydraulic_variable_names)
134	# Variable to hold the riverwall hydraulic properties in a table
135	# number of rows = number of riverwalls which cover edges in the domain
136	# number of cols = number of hydraulic variable names
137	self.hydraulic_properties=numpy.array([ [default_float] ])
138
139	# Variable to hold the indices of riverwall edges
140	# len = number of riverwall edges in the domain
141	self.riverwall_edges=numpy.array([default_int])
142
143	# Input info
144	self.input_riverwall_geo=None
145	self.input_riverwallPar=None
146
147
148	def create_riverwalls(self, riverwalls, riverwallPar={ },
149	default_riverwallPar={ },
150	tol=1.0e-4, verbose=True,
151	output_dir=None):
152	"""Add riverwalls at chosen locations along the mesh
153
154	As of 22/04/2014, these only work with DE algorithms [for which the concept is natural]
155
156	The walls MUST EXACTLY COINCIDE with edges along the mesh
157
158	You can force the mesh to do this by passing riverwalls.values()
159	to the 'breaklines' argument in the function create_mesh_from_regions. You
160	may also need to set the maximum_triangle_area for regions, if the breaklines
161	split the region. Do this with the regionPtArea argument in
162	create_mesh_from_regions.
163
164	As of 22/04/2014, the computational method used here is very 'brute-force'
165	For each segment on each riverwall, every point in the mesh is checked to see
166	if it is on the segment. A possibly faster/less memory method would be to
167	'walk' through connected edges on the mesh.
168
169	Inputs:
170	riverwalls: Dictionary of '3D polygons', containing xyz locations of named riverwalls.
171
172	exampleRiverWalls = { # Left bank n1 --
173	'n1': [[1.0, 1000., 2.],
174	[1.0, 50., 3.]],
175	# Left bank n2
176	'n2': [[2.0, 30., 1.0],
177	[3.0, 20., 2.],
178	[2.5, 15., 1.5]]
179	}
180
181	riverwallPar: Dictionary containing a dictionary of named hydraulic parameters for each named riverwall
182	If parameters are not provided, default values will be used.
183	See the help for class 'RiverWall' for an explanation of these
184
185	exampleRiverWallPar = {'n2': {'Qfactor':0.5} }
186	This would use a Qfactor of 0.5 for riverwall 'n2', while the other riverwall would have the default values
187
188	default_riverwallPar: Dictionary containing default values of the riverwall parameters, to be used
189	if they are not explicitly set.
190	If not provided, defaults from __init__ are used. See the help for class 'RiverWall' for more info
191
192	example_default_riverwallPar = {'Qfactor':1.5,
193	's1': 0.9,
194	's2': 0.95,
195	'h1': 1.0,
196	'h2': 1.5
197	}
198
199	example_default_riverwallPar = {'Qfactor':1.5,
200	's1': 10000.,
201	's2': 20000.
202	} # Here h1/h2 defaults will come from __init__
203
204
205	tol: Edges will be assigned a riverwall elevation if they are within 'tol' of
206	a segment in riverwalls. Round-off error means this should not be too small.
207
208	verbose: TRUE/FALSE Print lots of information
209
210	output_dir: Text representation of riverwalls is written to output_dir, unless output_dir=None
211
212	Outputs:
213	None, but it sets domain.edge_flux_type = 1 for every edge on a riverwall
214	Also, it adds a RiverWall object domain.riverwallData to the domain
215	The latter contains the riverwall heights, names, and hydraulic properties for each edge, in
216	a way we can pass in/out of C code.
217
218	"""
219
220	if(verbose and myid==0):
221	print ' '
222	print 'WARNING: Riverwall is an experimental feature'
223	print ' At each riverwall edge, we place a thin "wall" between'
224	print ' the 2 edges -- so each sees its neighbour edge as having'
225	print ' bed elevation = max(levee height, neighbour bed elevation)'
226	print ' We also force the discharge with a weir relation, blended with'
227	print ' the shallow water flux solution as the ratio (min_head)/(weir_height) becomes '
228	print ' large, or the ratio (downstream_head)/(upstream_head) becomes large'
229	print ' '
230	print ' It works in parallel, but you must use domain.riverwallData.create_riverwall AFTER distributing the mesh'
231	print ' '
232
233	# NOTE: domain.riverwallData is initialised in shallow_water_domain.py for DE algorithms
234	domain=self.domain
235
236
237	# Check flow algorithm
238	if(not domain.get_using_discontinuous_elevation()):
239	raise Exception, 'Riverwalls are currently only supported for discontinuous elevation flow algorithms'
240
241	if(len(self.names)>0):
242	# Delete any existing riverwall data
243	# The function cannot presently be used to partially edit riverwall data
244	if(verbose):
245	print 'Warning: There seems to be existing riverwall data'
246	print 'It will be deleted and overwritten with this function call'
247	domain.riverwallData.__init__(domain)
248
249	# Store input parameters
250	# FIXME: Is this worth it?
251	self.input_riverwall_geo=riverwalls
252	self.input_riverwallPar=riverwallPar
253
254	# Update self.default_riverwallPar (defined in __init__)
255	for i in self.default_riverwallPar.keys():
256	if(default_riverwallPar.has_key(i)):
257	self.default_riverwallPar[i]=default_riverwallPar[i]
258
259	# Check that all the keys in default_riverwallPar are allowed
260	for i in default_riverwallPar.keys():
261	if(not self.default_riverwallPar.has_key(i)):
262	msg='Key ', i + ' in default_riverwallPar not recognized'
263	raise Exception, msg
264	# Final default river-wall parameters
265	default_riverwallPar=self.default_riverwallPar
266
267	# Check that all named inputs in riverwallPar correspond to names in
268	# riverwall
269	for i in riverwallPar.keys():
270	if not riverwalls.has_key(i):
271	msg= 'Key ', i, ' in riverwallPar has no corresponding key in riverwall'
272	raise Exception, msg
273	#
274	# Check that all hydraulic parameter names in riverwallPar correspond
275	# to names in default_riverwallPar
276	#
277	for j in riverwallPar[i].keys():
278	if not default_riverwallPar.has_key(j):
279	msg = 'Hydraulic parameter named ', j ,\
280	' not recognised in default_riverwallPar'
281	raise Exception, msg
282
283	if(verbose):
284	print 'Setting riverwall elevations (P'+str(myid)+')...'
285
286	# Set up geometry
287	exy=domain.edge_coordinates
288	llx=domain.mesh.geo_reference.get_xllcorner()
289	lly=domain.mesh.geo_reference.get_yllcorner()
290
291	# Temporary variables
292	from anuga.config import max_float
293	riverwall_elevation=exy[:,0]*0. - max_float
294	riverwall_Qfactor=exy[:,0]*0.
295	riverwall_rowIndex=exy[:,0]*0 -1.
296	riverwall_rowIndex.astype(int)
297
298	# Loop over all segments in each riverwall, and set its elevation
299	nw=range(len(riverwalls))
300	nw_names=riverwalls.keys() # Not guarenteed to be in deterministic order
301
302	if(verbose):
303	# Use variable to record messages, allows cleaner parallel printing
304	printInfo=''
305
306	for i in nw:
307	# Name of riverwall
308	riverwalli_name=nw_names[i]
309	# Look at the ith riverwall
310	riverwalli=riverwalls[riverwalli_name]
311
312	ns=len(riverwalli)-1
313
314	if(verbose):
315	printInfo=printInfo + ' Wall ' + str(i) +' ....\n'
316
317	for j in range(ns):
318	if(verbose):
319	printInfo=printInfo + ' Segment ' + str(j) +' ....\n'
320	# Start and end xyz coordinates
321	start=riverwalli[j]
322	end=riverwalli[j+1]
323
324	if(len(start)!=3 \| len(end)!=3):
325	msg='Each riverwall coordinate must have at exactly 3 values [xyz]'
326	raise Exception, msg
327
328	# Find length
329	segLen=( (start[0]-end[0])2+(start[1]-end[1])2)**0.5
330	if(segLen<tol):
331	if(verbose):
332	printInfo=printInfo+' Segment with length < tolerance ' + str(tol) +' ignored\n'
333	continue
334
335	# Find edge indices which are within 'tol' of the segment
336	# We use a simple, inefficient method [but likely ok in practice
337	# except for very complex riverwalls]
338
339	# Unit vector along segment
340	se_0=-(start[0]-end[0])/segLen
341	se_1=-(start[1]-end[1])/segLen
342
343	# Vector from 'start' to every point on mesh
344	# NOTE: We account for georeferencing
345	pv_0 = exy[:,0]-(start[0]-llx)
346	pv_1 = exy[:,1]-(start[1]-lly)
347
348	pvLen=( pv_02 + pv_12)**0.5
349
350	# Dot product of pv and se == along-segment distance of projection
351	# of each point onto segment
352	pv_dot_se = pv_0se_0+pv_1se_1
353	# Normal distance^2 of each point to segment
354	perp_len_sq = pvLen2.-pv_dot_se2.
355
356	# Point is on a levee if the perpendicular distance is < tol,
357	# AND it is between start and end [within a tolerance]
358	onLevee=(perp_len_sq<tol*2)(pv_dot_se > 0.-tol)*(pv_dot_se<segLen+tol)
359	onLevee=onLevee.nonzero()
360	onLevee=onLevee[0]
361	if(len(onLevee)==0):
362	continue
363
364	if(verbose):
365	printInfo=printInfo+' Finding ' + str(len(onLevee)) + ' edges on this segment\n'
366
367	# Levee has Edge_flux_type=1
368	domain.edge_flux_type[onLevee]=1
369
370	# Get edge elevations as weighted averages of start/end elevations
371	w0=pv_dot_se[onLevee]/segLen
372	w0=w0*(w0>=0.0) # Enforce min of 0
373	w0=w0(w0<=1.0) + 1.0(w0>1.0) # Max of 1
374	riverwall_elevation[onLevee]= start[2](1.0-w0)+w0end[2]
375
376	# Record row index
377	riverwall_rowIndex[onLevee] = i
378
379	# Now, condense riverwall_elevation to array with length = number of riverwall edges
380	#
381	# We do this to avoid storing a riverwall_elevation for every edge in the mesh
382	# However, the data structure ends up being quite complex -- maybe there is a better way?
383	#
384	# The zeroth index in domain.edge_flux_type which = 1 will correspond to riverwall_elevation[0]
385	# The first index will correspond to riverwall_elevation[1]
386	# etc
387	#
388	riverwallInds=(domain.edge_flux_type==1).nonzero()[0]
389	# elevation
390	self.riverwall_elevation=\
391	riverwall_elevation[riverwallInds]
392	# corresponding row in the hydraulic properties table
393	self.hydraulic_properties_rowIndex=\
394	riverwall_rowIndex[riverwallInds].astype(int)
395	# index of edges which are riverwalls
396	self.riverwall_edges=riverwallInds
397
398	# Record the names of the riverwalls
399	self.names=nw_names
400
401
402	# Now create the hydraulic properties table
403
404	# Temporary variable to hold hydraulic properties table
405	# This will have as many rows are there are distinct riverwalls,
406	# and as many columns as there are hydraulic variables
407	hydraulicTmp=numpy.zeros((len(riverwalls), len(default_riverwallPar)))*numpy.nan
408
409	if(verbose):
410	print ' '
411	# Loop over every riverwall / hydraulic parameter, and set its value
412	for i in nw:
413	# Get the riverwall's name and specified parameters
414	name_riverwalli=nw_names[i]
415	if(riverwallPar.has_key(name_riverwalli)):
416	riverwalli_Par=riverwallPar[name_riverwalli]
417	else:
418	riverwalli_Par=None
419
420	# Set the ith riverwall's hydraulic properties
421	for j, hydraulicVar in enumerate(self.hydraulic_variable_names):
422	if((riverwalli_Par is not None) and (riverwalli_Par.has_key(hydraulicVar))):
423	if(verbose):
424	printInfo=printInfo+ ' Using provided '+ str(hydraulicVar)+' '+\
425	str(riverwalli_Par[hydraulicVar])+ ' for riverwall '+ str(name_riverwalli)+'\n'
426	hydraulicTmp[i,j]=riverwalli_Par[hydraulicVar]
427	else:
428	if(verbose):
429	printInfo=printInfo+ ' Using default '+ str(hydraulicVar)+' '+\
430	str(default_riverwallPar[hydraulicVar])+' for riverwall '+ str(name_riverwalli)+'\n'
431	hydraulicTmp[i,j]=default_riverwallPar[hydraulicVar]
432
433	if(verbose):
434	print ' '
435
436	# Check that s1 < s2
437	for i in nw:
438	if(hydraulicTmp[i,1]>= hydraulicTmp[i,2]):
439	msg = 's1 >= s2 on riverwall ' + nw_names[i] +'. This is not allowed'
440	raise Exception, msg
441	if( (hydraulicTmp[i,1]<0.) or (hydraulicTmp[i,2] < 0.)):
442	raise Exception, 's1 and s2 must be positive, with s1<s2'
443
444	# Check that h1 < h2
445	for i in nw:
446	if(hydraulicTmp[i,3]>= hydraulicTmp[i,4]):
447	msg = 'h1 >= h2 on riverwall ' + nw_names[i] +'. This is not allowed'
448	raise Exception, msg
449	if((hydraulicTmp[i,3]<0.) or (hydraulicTmp[i,4] < 0.)):
450	raise Exception, 'h1 and h2 must be positive, with h1<h2'
451
452	# Define the hydraulic properties
453	self.hydraulic_properties=hydraulicTmp
454
455	# Check for riverwall 'connectedness' errors (e.g. theoretically possible
456	# to miss an edge due to round-off)
457	connectedness=self.check_riverwall_connectedness(verbose=verbose)
458
459	self.export_riverwalls_to_text(output_dir=output_dir)
460
461	# Pretty printing of riverwall information in parallel
462	if(verbose):
463	if domain.parallel : barrier()
464	for i in range(numprocs):
465	if(myid==i):
466	print 'Processor '+str(myid)
467	print printInfo
468	print connectedness[0]
469	msg='Riverwall discontinuity -- possible round-off error in'+\
470	'finding edges on wall -- try increasing value of tol'
471	if(not connectedness[1]):
472	raise Exception, msg
473	if domain.parallel : barrier()
474	return
475
476	#####################################################################################
477
478	def get_centroids_corresponding_to_edgeInds(self, riverwalledgeInds):
479	"""
480	Get indices of centroids containing edges with indices riverwalledgeInds
481	"""
482	riverwallCentInds=numpy.floor(riverwalledgeInds/3.)
483	riverwallCentInds=riverwallCentInds.astype(int)
484
485	return riverwallCentInds
486
487	#####################################################################################
488
489	def get_vertices_corresponding_to_edgeInds(self, riverwalledgeInds, checkCoords=True):
490	"""
491	Get indices of vertices corresponding to edges at index riverwalledgeInds
492
493	Since each edge has 2 vertices, use V1 and V2
494
495	There is indeed a simple relationship between the vertex and edge indices
496	"""
497
498
499	#riverwallCentInds=self.get_centroids_corresponding_to_edgeInds(riverwalledgeInds)
500
501	rwEI_mod3=riverwalledgeInds%3
502
503	# Figure out the difference between each vertex index and the edge index. Is either
504	# -2, -1, 1, 2
505	rwV1_adjustment= -2(rwEI_mod3==2) -1(rwEI_mod3==1) +1*(rwEI_mod3==0)
506	rwV2_adjustment= -1(rwEI_mod3==2) +1(rwEI_mod3==1) +2*(rwEI_mod3==0)
507	riverwallV1Inds=riverwalledgeInds+rwV1_adjustment
508	riverwallV2Inds=riverwalledgeInds+rwV2_adjustment
509
510	if(checkCoords):
511	####################################################
512	# Check that vertices and edges really do correspond
513	domain=self.domain
514	# X coordinates
515	assert( numpy.allclose(
516	domain.edge_coordinates[riverwalledgeInds,0],
517	0.5*(domain.vertex_coordinates[riverwallV1Inds,0]+domain.vertex_coordinates[riverwallV2Inds,0]))
518	)
519	# Y coordinates
520	assert( numpy.allclose(
521	domain.edge_coordinates[riverwalledgeInds,1],
522	0.5*(domain.vertex_coordinates[riverwallV1Inds,1]+domain.vertex_coordinates[riverwallV2Inds,1]))
523	)
524	####################################################
525
526	return riverwallV1Inds, riverwallV2Inds
527
528	#####################################################################################
529	def is_vertex_on_boundary(self, vertexIndices):
530	"""
531	Determine whether a vertex is on the boundary of the domain
532	(i.e. it's connected with an edge that is a boundary edge)
533
534	INPUTS: self -- riverwallData
535	vertexIndices -- indices of vertices on the domain which are on the riverwall
536
537	OUTPUT:
538	TRUE if the vertex is on a domain boundary, FALSE otherwise
539
540	"""
541	domain=self.domain
542
543	# Get edge/vertex indices for boundaries
544	boundary_index_info=domain.boundary.keys()
545	boundary_edges=[ boundary_index_info[i][0]*3+boundary_index_info[i][1] for i in range(len(boundary_index_info))]
546	boundary_edges=numpy.array(boundary_edges)
547	tmp=self.get_vertices_corresponding_to_edgeInds(boundary_edges, checkCoords=False)
548	boundary_vertices=numpy.hstack([tmp[0], tmp[1]]).tolist()
549
550	# Get 'unique' vertex coordinates on boundary
551	node_complex=domain.vertex_coordinates[boundary_vertices,0]+1j*domain.vertex_coordinates[boundary_vertices,1]
552
553	# Get riverwall vertex coordinates as complex numbers (for equality testing)
554	complex_vertex_coords=domain.vertex_coordinates[vertexIndices.tolist(),0]+\
555	1j*domain.vertex_coordinates[vertexIndices.tolist(),1]
556
557	# Flag telling us if the vertex is on the boundary (1=True, 0=False)
558	isOnBoundary=[ 1 if complex_vertex_coords[i] in node_complex else 0 for i in range(len(complex_vertex_coords))]
559	isOnBoundary=numpy.array(isOnBoundary)
560
561	return isOnBoundary
562
563	#####################################################################################
564	def check_riverwall_connectedness(self, verbose=True):
565	"""
566	We expect riverwalls to be connected
567	(although they can pass through the bounding polygon several times, especially in parallel)
568	Round-off can potentially cause riverwalls to be dis-connected
569	Use this routine to check for that
570
571	Basically, we search for edges which are connected to vertices which
572	themselves are not connected to other edges. We ignore vertices on the domain's bounding-polygon
573
574	For a continuous riverwall, there can be at most 2 endpoints inside the domain
575
576	Otherwise, the riverwall is discontinuous, which is an error
577
578	"""
579
580	domain=self.domain
581
582	# Preliminary definitions
583	isConnected=True
584	printInfo=''
585
586	if(len(self.names)==0):
587	if(verbose):
588	printInfo=printInfo+' There are no riverwalls (P'+str(myid)+')\n'
589	return [printInfo, isConnected]
590
591	# Shorthand notation
592	rwd=self
593
594	for i, name in enumerate(rwd.names):
595	# Get indices of edges on this riverwall
596	riverwalledgeInds=rwd.riverwall_edges[(rwd.hydraulic_properties_rowIndex==i).nonzero()[0]]
597
598	if(len(riverwalledgeInds)==0):
599	printInfo=printInfo+'Riverwall '+name+' was not found on this mesh (if this is wrong, adjust tol in create_riverwalls)\n'
600	continue
601	# Get their corresponding vertices
602	riverwallV1Inds, riverwallV2Inds = rwd.get_vertices_corresponding_to_edgeInds(riverwalledgeInds)
603
604	# Flag telling us if vertex points are on the boundary of the model
605	# Used to help detect disconnected riverwalls (due to round-off)
606	v1_on_boundary=rwd.is_vertex_on_boundary(riverwallV1Inds)
607	v2_on_boundary=rwd.is_vertex_on_boundary(riverwallV2Inds)
608
609	# With discontinuous triangles, we expect edges to occur twice
610	# Let's remove duplicates to simplify the analysis
611	repeat=riverwalledgeInds*0
612	lre=len(riverwalledgeInds)
613	# edge coordinates as a complex number, for easy equality checking
614	complex_edge_coordinates=domain.edge_coordinates[riverwalledgeInds,0]+\
615	1j*domain.edge_coordinates[riverwalledgeInds,1]
616	for j in range(lre-1):
617	# Ignore if already checked
618	if(repeat[j]==1):
619	continue
620	# Check for a dupulicate
621	dups=(complex_edge_coordinates[(j+1):lre]==complex_edge_coordinates[j]).nonzero()[0]
622	if(len(dups)>0):
623	repeat[dups+j+1]=1
624
625	unique_riverwall_edge_indices=(repeat==0).nonzero()[0]
626
627	# Finally, get 'unqiue' edges in the riverwall
628	uEdges=riverwalledgeInds[unique_riverwall_edge_indices]
629	uV1=riverwallV1Inds[unique_riverwall_edge_indices]
630	uV2=riverwallV2Inds[unique_riverwall_edge_indices]
631	uV1_boundary=v1_on_boundary[unique_riverwall_edge_indices]
632	uV2_boundary=v2_on_boundary[unique_riverwall_edge_indices]
633
634	# Next, count how many times each vertex value occurs.
635	# For a 'connected' riverwall, we only want 2 edges where a vertex occurs only once,
636	# unless the vertex is on the boundary of the domain
637	lure=len(uEdges)
638	complex_v1_coordinates=domain.vertex_coordinates[uV1,0]+\
639	1j*domain.vertex_coordinates[uV1,1]
640	complex_v2_coordinates=domain.vertex_coordinates[uV2,0]+\
641	1j*domain.vertex_coordinates[uV2,1]
642	v1Counter=uEdges*0
643	v2Counter=uEdges*0
644	for j in range(lure):
645	v1Counter[j]=(complex_v1_coordinates==complex_v1_coordinates[j]).sum()+\
646	(complex_v2_coordinates==complex_v1_coordinates[j]).sum()
647	v2Counter[j]=(complex_v1_coordinates==complex_v2_coordinates[j]).sum()+\
648	(complex_v2_coordinates==complex_v2_coordinates[j]).sum()
649
650	num_disconnected_edges=((v1Counter==1)*(1-uV1_boundary)).sum()+\
651	((v2Counter==1)*(1-uV2_boundary)).sum()
652
653	if(verbose):
654	printInfo=printInfo+ ' On riverwall '+ str(name) +' there are '+ str(num_disconnected_edges)+\
655	' endpoints inside the domain [ignoring points on the boundary polygon] (P'+str(myid)+')\n'
656
657	if(num_disconnected_edges<=2):
658	if(verbose):
659	pass
660	#printInfo=printInfo+ " This is consistent with a continuous wall \n"
661	else:
662	isConnected=False
663	printInfo=printInfo+' Riverwall ' + name +' appears to be discontinuous. (P'+str(myid)+')\n'+\
664	' This suggests there is a gap in the wall, which should not occur\n'
665
666	return [printInfo, isConnected]
667
668	###################################################################
669
670	def export_riverwalls_to_text(self, output_dir=None):
671	"""
672	Function for dumping riverwall outputs to text file (useful for QC)
673
674	This will overwrite any existing files with the same location/name
675
676	INPUT: output_dir = Directory where the outputs will be written
677
678	OUTPUT:
679	None, but writes files as a side effect
680
681	"""
682	if(output_dir is None):
683	return
684
685	if(myid==0):
686	# Make output directory
687	try:
688	os.mkdir(output_dir)
689	except:
690	pass
691	# Make output files with empty contents
692	for i, riverWallFile in enumerate(self.names):
693	newFile=open(output_dir+'/'+riverWallFile+'.txt','w')
694	# Write hydraulic variable information
695	hydraulicVars=self.hydraulic_properties[i,:]
696	newFile.write('## Hydraulic Variable values below ## \n')
697	newFile.write(str(self.hydraulic_variable_names)+'\n')
698	newFile.write(str(hydraulicVars)+'\n')
699	newFile.write('\n')
700	newFile.write('## xyElevation at edges below. Order may be erratic for parallel runs ##\n')
701	newFile.close()
702	else:
703	pass
704
705	# Files must be created before writing out
706	domain=self.domain
707
708	#if domain.parallel: barrier()
709
710	# Now dump the required info to the files
711	for i in range(numprocs):
712	# Write 1 processor at a time
713	if(myid==i):
714	for j, riverWallname in enumerate(self.names):
715	# Get xyz data for riverwall j
716	riverWallInds=(self.hydraulic_properties_rowIndex==j).nonzero()[0].tolist()
717	riverWallDomainInds=self.riverwall_edges[riverWallInds].tolist()
718	myXCoords=domain.edge_coordinates[riverWallDomainInds,0]+domain.geo_reference.xllcorner
719	myYCoords=domain.edge_coordinates[riverWallDomainInds,1]+domain.geo_reference.yllcorner
720	myElev=self.riverwall_elevation[riverWallInds]
721
722	# Open file for appending data
723	theFile=open(output_dir+'/'+riverWallname+'.txt','a')
724	for k in range(len(myElev)):
725	theFile.write(str(myXCoords[k])+','+str(myYCoords[k])+','+str(myElev[k])+'\n')
726	theFile.close()
727
728	else:
729	pass
730
731	#if domain.parallel : barrier()
732
733	return

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

Download in other formats: