Context Navigation

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

Last change on this file since 9679 was 9602, checked in by davies, 10 years ago
Updates
File size: 33.6 KB

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

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

Download in other formats: