| 1 | import numpy |
|---|
| 2 | |
|---|
| 3 | class RiverWall: |
|---|
| 4 | """Define the elevation of 'riverwalls'. |
|---|
| 5 | |
|---|
| 6 | These are located along each cell edge, and can have an elevation different |
|---|
| 7 | from the bed elevation. |
|---|
| 8 | |
|---|
| 9 | For the DE algorithms, they are used in computing the 'z_half' value [if they |
|---|
| 10 | are greater than either edge elevation]. |
|---|
| 11 | |
|---|
| 12 | In addition, the DE fluxes at riverwalls are adjusted to agree with a weir relation, |
|---|
| 13 | so long as the riverwall is not too deeply submerged. |
|---|
| 14 | |
|---|
| 15 | As of 22/04/2014, they are only implemented for DE0 and DE1 [the shallow water |
|---|
| 16 | component would be more difficult to implement with other algorithms] |
|---|
| 17 | |
|---|
| 18 | How the flux over the riverwall is computed: |
|---|
| 19 | # We have a riverwall, which is treated as a weir. |
|---|
| 20 | # |
|---|
| 21 | # Say the headwater-head is the upstream water depth above the weir elevation (min 0), and |
|---|
| 22 | # the tailwater head is the downstream water depth above the weir elevation (min 0). |
|---|
| 23 | # By definition headwater-head > tailwater-head. |
|---|
| 24 | # |
|---|
| 25 | # Let s = (headwater head) / (tailwater head), and h = (tailwater head)/(weir height) |
|---|
| 26 | # where 'weir height' is the elevation of the weir crest above the minimum |
|---|
| 27 | # of the left/right bed elevations |
|---|
| 28 | # |
|---|
| 29 | # Denote ID as the 'ideal' weir flow, computed from the hydraulic |
|---|
| 30 | # formula (including a submergence correction factor from Villemonte, |
|---|
| 31 | # 1947), |
|---|
| 32 | # Q1 = 2/3*headwater_head*sqrt(g*2/3*headwater_head)*Qfactor |
|---|
| 33 | # Q2 = 2/3*tailwater_head*sqrt(g*2/3*headwater_head)*Qfactor |
|---|
| 34 | # ID = Q1*(1-Q2/Q1)**0.385 |
|---|
| 35 | # |
|---|
| 36 | # Denote SW as the 'shallow-water' weir flux, computed from the approximate |
|---|
| 37 | # reimann solver, where the mid-edge-elevation is the weir crest elevation. |
|---|
| 38 | # This makes clear sense for DE0 and DE1. Cell centroid stage/height/bed_elevation |
|---|
| 39 | # are used in the flux computation |
|---|
| 40 | # |
|---|
| 41 | # Then the flux over the weir is computed from: |
|---|
| 42 | # w1 = min( max(s-s1, 0.)/(s2-s1), 1.0) # Factor describing relative submergence |
|---|
| 43 | # w1' = min( max(h-h1,0.)/(h2-h1), 1.0) # Factor describing absolute submergence |
|---|
| 44 | # flux_over_weir = (w1*SW + (1-w1)*ID)*( 1-w1') + (w1')*SW |
|---|
| 45 | # |
|---|
| 46 | # The key idea is that if s<s1, h<h1, then the ideal weir solution is |
|---|
| 47 | # used. Otherwise, we linearly blend with the SW solution, |
|---|
| 48 | # and the SW solution is used completely if s>s2 or h>h2 |
|---|
| 49 | |
|---|
| 50 | """ |
|---|
| 51 | |
|---|
| 52 | def __init__(self, domain): |
|---|
| 53 | """Riverwall data structure |
|---|
| 54 | |
|---|
| 55 | Allows reasonably efficient storage of riverwall elevations and hydraulic |
|---|
| 56 | properties |
|---|
| 57 | |
|---|
| 58 | NOTE: In domain.edge_flux_type, riverwalls correspond to the value 1 |
|---|
| 59 | |
|---|
| 60 | RiverWall variables are initialised to dummy values, because even |
|---|
| 61 | if no riverwalls are used, some values have to be passed to the flux computation |
|---|
| 62 | |
|---|
| 63 | Hydraulic parameters are |
|---|
| 64 | Qfactor -- Multiplicative factor for ideal weir relation (calibration coef) |
|---|
| 65 | s1 -- Submergence ratio at which we start blending with the shallow water solution (<s2) |
|---|
| 66 | s2 -- Submergence ratio at which we entirely use the shallow water solution (>s1) |
|---|
| 67 | h1 -- Tailwater head / weir height at which we start blending with the shallow water solution (<h2) |
|---|
| 68 | h2 -- Tailwater head / weir height at which we entirely use the shallow water solution (>h1) |
|---|
| 69 | |
|---|
| 70 | # Default riverwall hydraulic parameters |
|---|
| 71 | default_riverwallPar={'Qfactor':1.0, |
|---|
| 72 | 's1': 0.9, |
|---|
| 73 | 's2': 0.95, |
|---|
| 74 | 'h1': 1.0, |
|---|
| 75 | 'h2': 1.5 |
|---|
| 76 | } |
|---|
| 77 | |
|---|
| 78 | Other variables are: |
|---|
| 79 | |
|---|
| 80 | riverwall_elevation -- Variable to hold riverwall elevations. |
|---|
| 81 | len = number of riverwall edges in the domain. |
|---|
| 82 | First entry corresponds to first riverwall edge in domain.edge_coordinates, |
|---|
| 83 | second entry to second riverwall edge in domain.edge_coordinates, etc |
|---|
| 84 | |
|---|
| 85 | hydraulic_properties_rowIndex -- Variable to hold the row-index of the hydraulic properties table |
|---|
| 86 | len = number of riverwall edges |
|---|
| 87 | in the domain, ordered like riverwall_elevation |
|---|
| 88 | |
|---|
| 89 | riverwall_edges -- Holds indices of edges in domain which are riverwalls, ordered like riverwall_elevation |
|---|
| 90 | |
|---|
| 91 | names -- list with the names of the riverwalls |
|---|
| 92 | len = number of riverwalls which cover edges in the domain |
|---|
| 93 | |
|---|
| 94 | hydraulic_variable_names -- Variables to hold the names of variables in columns of the hydraulic |
|---|
| 95 | properties table. THE ORDER IS IMPORTANT -- C code ASSUMES IT |
|---|
| 96 | |
|---|
| 97 | ncol_hydraulic_properties -- Number of columns in the hydraulic properties table |
|---|
| 98 | [ = len(hydraulic_variable_names) ] |
|---|
| 99 | |
|---|
| 100 | hydraulic_properties -- Array with the hydraulic parameters for each riverwall. |
|---|
| 101 | number of rows = number of riverwalls which cover edges in the domain |
|---|
| 102 | number of cols = number of hydraulic variable names |
|---|
| 103 | |
|---|
| 104 | |
|---|
| 105 | input_riverwall_geo, input_riverwall_par -- holds input information |
|---|
| 106 | |
|---|
| 107 | """ |
|---|
| 108 | self.domain=domain |
|---|
| 109 | |
|---|
| 110 | default_float=-9.0e+20 |
|---|
| 111 | default_int=-9e+20 |
|---|
| 112 | self.riverwall_elevation=numpy.array([default_float]) |
|---|
| 113 | |
|---|
| 114 | self.hydraulic_properties_rowIndex=numpy.array([default_int]) |
|---|
| 115 | |
|---|
| 116 | self.names=[ ] |
|---|
| 117 | |
|---|
| 118 | # Default riverwall hydraulic parameters |
|---|
| 119 | self.default_riverwallPar={'Qfactor':1.0, |
|---|
| 120 | 's1': 0.9, |
|---|
| 121 | 's2': 0.95, |
|---|
| 122 | 'h1': 1.0, |
|---|
| 123 | 'h2': 1.5 |
|---|
| 124 | } |
|---|
| 125 | |
|---|
| 126 | # DO NOT CHANGE THE ORDER OF hydraulic_variable_names |
|---|
| 127 | # It needs to match hard-coded assumptions in C [compute_fluxes_central] |
|---|
| 128 | # If you add a variable, append it to the end of hydraulic_variable_names |
|---|
| 129 | self.hydraulic_variable_names=('Qfactor', 's1', 's2', 'h1', 'h2') |
|---|
| 130 | |
|---|
| 131 | self.ncol_hydraulic_properties=len(self.hydraulic_variable_names) |
|---|
| 132 | # Variable to hold the riverwall hydraulic properties in a table |
|---|
| 133 | # number of rows = number of riverwalls which cover edges in the domain |
|---|
| 134 | # number of cols = number of hydraulic variable names |
|---|
| 135 | self.hydraulic_properties=numpy.array([ [default_float] ]) |
|---|
| 136 | |
|---|
| 137 | # Variable to hold the indices of riverwall edges |
|---|
| 138 | # len = number of riverwall edges in the domain |
|---|
| 139 | self.riverwall_edges=numpy.array([default_int]) |
|---|
| 140 | |
|---|
| 141 | # Input info |
|---|
| 142 | self.input_riverwall_geo=None |
|---|
| 143 | self.input_riverwallPar=None |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | def create_riverwalls(self, riverwalls, riverwallPar={ }, |
|---|
| 147 | default_riverwallPar={ }, |
|---|
| 148 | tol=1.0e-04, verbose=True): |
|---|
| 149 | """Add riverwalls at chosen locations along the mesh |
|---|
| 150 | |
|---|
| 151 | As of 22/04/2014, these only work with DE0 and DE1 [for which the concept is natural] |
|---|
| 152 | |
|---|
| 153 | The walls MUST EXACTLY COINCIDE with edges along the mesh |
|---|
| 154 | |
|---|
| 155 | You can force the mesh to do this by passing riverwalls.values() |
|---|
| 156 | to the 'breaklines' argument in the function create_mesh_from_regions. You |
|---|
| 157 | may also need to set the maximum_triangle_area for regions, if the breaklines |
|---|
| 158 | split the region. Do this with the regionPtArea argument in |
|---|
| 159 | create_mesh_from_regions. |
|---|
| 160 | |
|---|
| 161 | As of 22/04/2014, the computational method used here is very 'brute-force' |
|---|
| 162 | For each segment on each riverwall, every point in the mesh is checked to see |
|---|
| 163 | if it is on the segment. A possibly faster/less memory method would be to |
|---|
| 164 | 'walk' through connected edges on the mesh. |
|---|
| 165 | |
|---|
| 166 | Inputs: |
|---|
| 167 | riverwalls: Dictionary of '3D polygons', containing xyz locations of named riverwalls. |
|---|
| 168 | |
|---|
| 169 | exampleRiverWalls = { # Left bank n1 -- |
|---|
| 170 | 'n1': [[1.0, 1000., 2.], |
|---|
| 171 | [1.0, 50., 3.]], |
|---|
| 172 | # Left bank n2 |
|---|
| 173 | 'n2': [[2.0, 30., 1.0], |
|---|
| 174 | [3.0, 20., 2.], |
|---|
| 175 | [2.5, 15., 1.5]] |
|---|
| 176 | } |
|---|
| 177 | |
|---|
| 178 | riverwallPar: Dictionary containing a dictionary of named hydraulic parameters for each named riverwall |
|---|
| 179 | If parameters are not provided, default values will be used. |
|---|
| 180 | See the help for class 'RiverWall' for an explanation of these |
|---|
| 181 | |
|---|
| 182 | exampleRiverWallPar = {'n2': {'Qfactor':0.5} } |
|---|
| 183 | This would use a Qfactor of 0.5 for riverwall 'n2', while the other riverwall would have the default values |
|---|
| 184 | |
|---|
| 185 | default_riverwallPar: Dictionary containing default values of the riverwall parameters, to be used |
|---|
| 186 | if they are not explicitly set. |
|---|
| 187 | If not provided, defaults from __init__ are used. See the help for class 'RiverWall' for more info |
|---|
| 188 | |
|---|
| 189 | example_default_riverwallPar = {'Qfactor':1.5, |
|---|
| 190 | 's1': 0.9, |
|---|
| 191 | 's2': 0.95, |
|---|
| 192 | 'h1': 1.0, |
|---|
| 193 | 'h2': 1.5 |
|---|
| 194 | } |
|---|
| 195 | |
|---|
| 196 | example_default_riverwallPar = {'Qfactor':1.5, |
|---|
| 197 | 's1': 10000., |
|---|
| 198 | 's2': 20000. |
|---|
| 199 | } # Here h1/h2 defaults will come from __init__ |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | tol: Edges will be assigned a riverwall elevation if they are within 'tol' of |
|---|
| 203 | a segment in riverwalls. Round-off error means this should not be too small. |
|---|
| 204 | |
|---|
| 205 | verbose: TRUE/FALSE Print lots of information |
|---|
| 206 | |
|---|
| 207 | Outputs: |
|---|
| 208 | None, but it sets domain.edge_flux_type = 1 for every edge on a riverwall |
|---|
| 209 | Also, it adds a RiverWall object domain.riverwallData to the domain |
|---|
| 210 | The latter contains the riverwall heights, names, and hydraulic properties for each edge, in |
|---|
| 211 | a way we can pass in/out of C code. |
|---|
| 212 | |
|---|
| 213 | """ |
|---|
| 214 | |
|---|
| 215 | if(verbose): |
|---|
| 216 | print ' ' |
|---|
| 217 | print 'WARNING: Riverwall is an experimental feature' |
|---|
| 218 | print ' At each riverwall edge, we place a thin "wall" between' |
|---|
| 219 | print ' the 2 edges -- so each sees its neighbour edge as having' |
|---|
| 220 | print ' bed elevation = max(levee height, neighbour bed elevation)' |
|---|
| 221 | print ' We also force the discharge with a weir relation, blended with' |
|---|
| 222 | print ' the shallow water flux solution as the ratio (min_head)/(weir_height) becomes ' |
|---|
| 223 | print ' large, or the ratio (downstream_head)/(upstream_head) becomes large' |
|---|
| 224 | print ' ' |
|---|
| 225 | print ' It works in parallel, but you must use domain.riverwallData.create_riverwall AFTER distributing the mesh' |
|---|
| 226 | print ' ' |
|---|
| 227 | |
|---|
| 228 | # NOTE: domain.riverwallData is initialised in shallow_water_domain.py for DE0 and DE1 |
|---|
| 229 | |
|---|
| 230 | domain=self.domain |
|---|
| 231 | |
|---|
| 232 | |
|---|
| 233 | # Check flow algorithm |
|---|
| 234 | if(not domain.get_using_discontinuous_elevation()): |
|---|
| 235 | raise Exception, 'Riverwalls are currently only supported when domain.flow_algorithm="DE0" and "DE1"' |
|---|
| 236 | |
|---|
| 237 | if(len(self.names)>0): |
|---|
| 238 | # Delete any existing riverwall data |
|---|
| 239 | # The function cannot presently be used to partially edit riverwall data |
|---|
| 240 | if(verbose): |
|---|
| 241 | print 'Warning: There seems to be existing riverwall data' |
|---|
| 242 | print 'It will be deleted and overwritten with this function call' |
|---|
| 243 | domain.riverwallData.__init__(domain) |
|---|
| 244 | |
|---|
| 245 | # Store input parameters |
|---|
| 246 | # FIXME: Is this worth it? |
|---|
| 247 | self.input_riverwall_geo=riverwalls |
|---|
| 248 | self.input_riverwallPar=riverwallPar |
|---|
| 249 | |
|---|
| 250 | # Update self.default_riverwallPar |
|---|
| 251 | for i in self.default_riverwallPar.keys(): |
|---|
| 252 | if(default_riverwallPar.has_key(i)): |
|---|
| 253 | self.default_riverwallPar[i]=default_riverwallPar[i] |
|---|
| 254 | |
|---|
| 255 | # Check that all the keys in default_riverwallPar are allowed |
|---|
| 256 | for i in default_riverwallPar.keys(): |
|---|
| 257 | if(not self.default_riverwallPar.has_key(i)): |
|---|
| 258 | msg='Key ', i + ' in default_riverwallPar not recognized' |
|---|
| 259 | raise Exception, msg |
|---|
| 260 | # Final default river-wall parameters |
|---|
| 261 | default_riverwallPar=self.default_riverwallPar |
|---|
| 262 | |
|---|
| 263 | # Check that all named inputs in riverwallPar correspond to names in |
|---|
| 264 | # riverwall |
|---|
| 265 | for i in riverwallPar.keys(): |
|---|
| 266 | if not riverwalls.has_key(i): |
|---|
| 267 | msg= 'Key ', i, ' in riverwallPar has no corresponding key in riverwall' |
|---|
| 268 | raise Exception, msg |
|---|
| 269 | # |
|---|
| 270 | # Check that all hydraulic parameter names in riverwallPar correspond |
|---|
| 271 | # to names in default_riverwallPar |
|---|
| 272 | # |
|---|
| 273 | for j in riverwallPar[i].keys(): |
|---|
| 274 | if not default_riverwallPar.has_key(j): |
|---|
| 275 | msg = 'Hydraulic parameter named ', j ,\ |
|---|
| 276 | ' not recognised in default_riverwallPar' |
|---|
| 277 | raise Exception, msg |
|---|
| 278 | |
|---|
| 279 | if(verbose): |
|---|
| 280 | print 'Setting riverwall elevations ...' |
|---|
| 281 | |
|---|
| 282 | # Set up geometry |
|---|
| 283 | exy=domain.edge_coordinates |
|---|
| 284 | llx=domain.mesh.geo_reference.get_xllcorner() |
|---|
| 285 | lly=domain.mesh.geo_reference.get_yllcorner() |
|---|
| 286 | |
|---|
| 287 | # Temporary variables |
|---|
| 288 | from anuga.config import max_float |
|---|
| 289 | riverwall_elevation=exy[:,0]*0. - max_float |
|---|
| 290 | riverwall_Qfactor=exy[:,0]*0. |
|---|
| 291 | riverwall_rowIndex=exy[:,0]*0 -1. |
|---|
| 292 | riverwall_rowIndex.astype(int) |
|---|
| 293 | |
|---|
| 294 | # Loop over all segments in each riverwall, and set its elevation |
|---|
| 295 | nw=range(len(riverwalls)) |
|---|
| 296 | nw_names=riverwalls.keys() # Not guarenteed to be in deterministic order |
|---|
| 297 | for i in nw: |
|---|
| 298 | # Name of riverwall |
|---|
| 299 | riverwalli_name=nw_names[i] |
|---|
| 300 | # Look at the ith riverwall |
|---|
| 301 | riverwalli=riverwalls[riverwalli_name] |
|---|
| 302 | |
|---|
| 303 | ns=len(riverwalli)-1 |
|---|
| 304 | |
|---|
| 305 | if(verbose): |
|---|
| 306 | print 'Wall ' + str(i) +' ....' |
|---|
| 307 | |
|---|
| 308 | for j in range(ns): |
|---|
| 309 | if(verbose): |
|---|
| 310 | print ' Segment ' + str(j) +' ....' |
|---|
| 311 | # Start and end xyz coordinates |
|---|
| 312 | start=riverwalli[j] |
|---|
| 313 | end=riverwalli[j+1] |
|---|
| 314 | |
|---|
| 315 | if(len(start)!=3 | len(end)!=3): |
|---|
| 316 | msg='Each riverwall coordinate must have at exactly 3 values [xyz]' |
|---|
| 317 | raise Exception, msg |
|---|
| 318 | |
|---|
| 319 | # Find length |
|---|
| 320 | segLen=( (start[0]-end[0])**2+(start[1]-end[1])**2)**0.5 |
|---|
| 321 | if(segLen<tol): |
|---|
| 322 | print 'Segment with length < tolerance ' + str(tol) +' ignored' |
|---|
| 323 | continue |
|---|
| 324 | |
|---|
| 325 | # Find edge indices which are within 'tol' of the segment |
|---|
| 326 | # We use a simple, inefficient method [but likely ok in practice |
|---|
| 327 | # except for very complex riverwalls] |
|---|
| 328 | |
|---|
| 329 | # Unit vector along segment |
|---|
| 330 | se_0=-(start[0]-end[0])/segLen |
|---|
| 331 | se_1=-(start[1]-end[1])/segLen |
|---|
| 332 | |
|---|
| 333 | # Vector from 'start' to every point on mesh |
|---|
| 334 | # NOTE: We account for georeferencing |
|---|
| 335 | pv_0 = exy[:,0]-(start[0]-llx) |
|---|
| 336 | pv_1 = exy[:,1]-(start[1]-lly) |
|---|
| 337 | |
|---|
| 338 | pvLen=( pv_0**2 + pv_1**2)**0.5 |
|---|
| 339 | |
|---|
| 340 | # Dot product of pv and se == along-segment distance of projection |
|---|
| 341 | # of each point onto segment |
|---|
| 342 | pv_dot_se = pv_0*se_0+pv_1*se_1 |
|---|
| 343 | # Normal distance^2 of each point to segment |
|---|
| 344 | perp_len_sq = pvLen**2.-pv_dot_se**2. |
|---|
| 345 | |
|---|
| 346 | # Point is on a levee if the perpendicular distance is < tol, |
|---|
| 347 | # AND it is between start and end [within a tolerance] |
|---|
| 348 | onLevee=(perp_len_sq<tol**2)*(pv_dot_se > 0.-tol)*(pv_dot_se<segLen+tol) |
|---|
| 349 | onLevee=onLevee.nonzero() |
|---|
| 350 | onLevee=onLevee[0] |
|---|
| 351 | if(len(onLevee)==0): |
|---|
| 352 | continue |
|---|
| 353 | |
|---|
| 354 | if(verbose): |
|---|
| 355 | print ' Finding ' + str(len(onLevee)) + ' edges on this segment ' |
|---|
| 356 | |
|---|
| 357 | # Levee has Edge_flux_type=1 |
|---|
| 358 | domain.edge_flux_type[onLevee]=1 |
|---|
| 359 | |
|---|
| 360 | # Get edge elevations as weighted averages of start/end elevations |
|---|
| 361 | w0=pv_dot_se[onLevee]/segLen |
|---|
| 362 | w0=w0*(w0>=0.0) # Enforce min of 0 |
|---|
| 363 | w0=w0*(w0<=1.0) + 1.0*(w0>1.0) # Max of 1 |
|---|
| 364 | riverwall_elevation[onLevee]= start[2]*(1.0-w0)+w0*end[2] |
|---|
| 365 | |
|---|
| 366 | # Record row index |
|---|
| 367 | riverwall_rowIndex[onLevee] = i |
|---|
| 368 | |
|---|
| 369 | # Now, condense riverwall_elevation to array with length = number of riverwall edges |
|---|
| 370 | # |
|---|
| 371 | # FIXME: We do this to avoid storing a riverwall_elevation for every edge in the mesh |
|---|
| 372 | # However, the data structure ends up being quite complex -- maybe there is a better way? |
|---|
| 373 | # |
|---|
| 374 | # The zeroth index in domain.edge_flux_type which = 1 will correspond to riverwall_elevation[0] |
|---|
| 375 | # The first index will correspond to riverwall_elevation[1] |
|---|
| 376 | # etc |
|---|
| 377 | # |
|---|
| 378 | riverwallInds=(domain.edge_flux_type==1).nonzero()[0] |
|---|
| 379 | # elevation |
|---|
| 380 | self.riverwall_elevation=\ |
|---|
| 381 | riverwall_elevation[riverwallInds] |
|---|
| 382 | # corresponding row in the hydraulic properties table |
|---|
| 383 | self.hydraulic_properties_rowIndex=\ |
|---|
| 384 | riverwall_rowIndex[riverwallInds].astype(int) |
|---|
| 385 | # index of edges which are riverwalls [FIXME: Not really needed, can easily back-calculate from edge_flux_type] |
|---|
| 386 | self.riverwall_edges=riverwallInds |
|---|
| 387 | |
|---|
| 388 | # Record the names of the riverwalls |
|---|
| 389 | self.names=nw_names |
|---|
| 390 | |
|---|
| 391 | |
|---|
| 392 | # Now create the hydraulic properties table |
|---|
| 393 | |
|---|
| 394 | # Temporary variable to hold hydraulic properties table |
|---|
| 395 | # This will have as many rows are there are distinct riverwalls, |
|---|
| 396 | # and as many columns as there are hydraulic variables |
|---|
| 397 | hydraulicTmp=numpy.zeros((len(riverwalls), len(default_riverwallPar)))*numpy.nan |
|---|
| 398 | |
|---|
| 399 | if(verbose): |
|---|
| 400 | print ' ' |
|---|
| 401 | # Loop over every riverwall / hydraulic parameter, and set its value |
|---|
| 402 | for i in nw: |
|---|
| 403 | # Get the riverwall's name and specified parameters |
|---|
| 404 | name_riverwalli=nw_names[i] |
|---|
| 405 | if(riverwallPar.has_key(name_riverwalli)): |
|---|
| 406 | riverwalli_Par=riverwallPar[name_riverwalli] |
|---|
| 407 | else: |
|---|
| 408 | riverwalli_Par=None |
|---|
| 409 | |
|---|
| 410 | # Set the ith riverwall's hydraulic properties |
|---|
| 411 | for j, hydraulicVar in enumerate(self.hydraulic_variable_names): |
|---|
| 412 | if((riverwalli_Par is not None) and (riverwalli_Par.has_key(hydraulicVar))): |
|---|
| 413 | if(verbose): |
|---|
| 414 | print 'Using provided ', hydraulicVar,\ |
|---|
| 415 | riverwalli_Par[hydraulicVar], ' for riverwall ', name_riverwalli |
|---|
| 416 | hydraulicTmp[i,j]=riverwalli_Par[hydraulicVar] |
|---|
| 417 | else: |
|---|
| 418 | if(verbose): |
|---|
| 419 | print 'Using default ', hydraulicVar,\ |
|---|
| 420 | default_riverwallPar[hydraulicVar], ' for riverwall ', name_riverwalli |
|---|
| 421 | hydraulicTmp[i,j]=default_riverwallPar[hydraulicVar] |
|---|
| 422 | |
|---|
| 423 | if(verbose): |
|---|
| 424 | print ' ' |
|---|
| 425 | |
|---|
| 426 | # Check that s1 < s2 |
|---|
| 427 | for i in nw: |
|---|
| 428 | if(hydraulicTmp[i,1]>= hydraulicTmp[i,2]): |
|---|
| 429 | msg = 's1 >= s2 on riverwall ' + nw_names[i] +'. This is not allowed' |
|---|
| 430 | raise Exception, msg |
|---|
| 431 | if( (hydraulicTmp[i,1]<0.) or (hydraulicTmp[i,2] < 0.)): |
|---|
| 432 | raise Exception, 's1 and s2 must be positive, with s1<s2' |
|---|
| 433 | |
|---|
| 434 | # Check that h1 < h2 |
|---|
| 435 | for i in nw: |
|---|
| 436 | if(hydraulicTmp[i,3]>= hydraulicTmp[i,4]): |
|---|
| 437 | msg = 'h1 >= h2 on riverwall ' + nw_names[i] +'. This is not allowed' |
|---|
| 438 | raise Exception, msg |
|---|
| 439 | if((hydraulicTmp[i,3]<0.) or (hydraulicTmp[i,4] < 0.)): |
|---|
| 440 | raise Exception, 'h1 and h2 must be positive, with h1<h2' |
|---|
| 441 | |
|---|
| 442 | # Define the hydraulic properties |
|---|
| 443 | self.hydraulic_properties=hydraulicTmp |
|---|
| 444 | |
|---|
| 445 | # Check for riverwall 'connectedness' errors (e.g. theoretically possible |
|---|
| 446 | # to miss an edge due to round-off) |
|---|
| 447 | self.check_riverwall_connectedness(domain, verbose,tol=tol) |
|---|
| 448 | |
|---|
| 449 | return() |
|---|
| 450 | |
|---|
| 451 | ##################################################################################### |
|---|
| 452 | |
|---|
| 453 | def get_centroids_corresponding_to_edgeInds(self, riverwalledgeInds): |
|---|
| 454 | """ |
|---|
| 455 | Get indices of centroids containing edges with indices riverwalledgeInds |
|---|
| 456 | """ |
|---|
| 457 | riverwallCentInds=numpy.floor(riverwalledgeInds/3.) |
|---|
| 458 | riverwallCentInds=riverwallCentInds.astype(int) |
|---|
| 459 | |
|---|
| 460 | return riverwallCentInds |
|---|
| 461 | |
|---|
| 462 | ##################################################################################### |
|---|
| 463 | |
|---|
| 464 | def get_vertices_corresponding_to_edgeInds(self, riverwalledgeInds, tol=1.0e-04): |
|---|
| 465 | """ |
|---|
| 466 | Get indices of vertices corresponding to edges at index riverwalledgeInds |
|---|
| 467 | Since each edge has 2 vertices, use V1 and V2 |
|---|
| 468 | We have to work to compute these, I can't see a simple formula |
|---|
| 469 | """ |
|---|
| 470 | |
|---|
| 471 | domain=self.domain |
|---|
| 472 | |
|---|
| 473 | riverwallCentInds=self.get_centroids_corresponding_to_edgeInds(riverwalledgeInds) |
|---|
| 474 | |
|---|
| 475 | rwV0Inds=3*riverwallCentInds #+((resid+1)%3) |
|---|
| 476 | rwV1Inds=3*riverwallCentInds+1 #((resid+2)%3) |
|---|
| 477 | rwV2Inds=3*riverwallCentInds+2 #((resid+2)%3) |
|---|
| 478 | |
|---|
| 479 | # Find which 2 vertices have the edge value, by testing the |
|---|
| 480 | # x coordinate. |
|---|
| 481 | X01=0.5*(domain.vertex_coordinates[rwV0Inds,0]+domain.vertex_coordinates[rwV1Inds,0]) |
|---|
| 482 | k01=(abs(X01- domain.edge_coordinates[riverwalledgeInds,0])<tol) |
|---|
| 483 | X12=0.5*(domain.vertex_coordinates[rwV2Inds,0]+domain.vertex_coordinates[rwV1Inds,0]) |
|---|
| 484 | k12=(abs(X12 - domain.edge_coordinates[riverwalledgeInds,0])<tol) |
|---|
| 485 | X02=0.5*(domain.vertex_coordinates[rwV2Inds,0]+domain.vertex_coordinates[rwV0Inds,0]) |
|---|
| 486 | k02=(abs(X02 - domain.edge_coordinates[riverwalledgeInds,0])<tol) |
|---|
| 487 | |
|---|
| 488 | riverwallV2Inds=riverwallCentInds*0 |
|---|
| 489 | riverwallV1Inds=riverwallCentInds*0 |
|---|
| 490 | |
|---|
| 491 | t1=(k01==1).nonzero()[0] |
|---|
| 492 | riverwallV2Inds[t1]=rwV0Inds[t1] |
|---|
| 493 | riverwallV1Inds[t1]=rwV1Inds[t1] |
|---|
| 494 | |
|---|
| 495 | t1=(k12==1).nonzero()[0] |
|---|
| 496 | riverwallV2Inds[t1]=rwV2Inds[t1] |
|---|
| 497 | riverwallV1Inds[t1]=rwV1Inds[t1] |
|---|
| 498 | |
|---|
| 499 | t1=(k02==1).nonzero()[0] |
|---|
| 500 | riverwallV2Inds[t1]=rwV0Inds[t1] |
|---|
| 501 | riverwallV1Inds[t1]=rwV2Inds[t1] |
|---|
| 502 | |
|---|
| 503 | return riverwallV1Inds, riverwallV2Inds |
|---|
| 504 | |
|---|
| 505 | ##################################################################################### |
|---|
| 506 | |
|---|
| 507 | def check_riverwall_connectedness(self, warn_only=False, verbose=True, tol=1.0e-04): |
|---|
| 508 | """ |
|---|
| 509 | We expect riverwalls to be connected [unless they pass over a hole] |
|---|
| 510 | However, round-off could potentially cause riverwalls to be dis-connected |
|---|
| 511 | Use this routine to check for that |
|---|
| 512 | |
|---|
| 513 | Basically, we search for edges which are connected to vertices which |
|---|
| 514 | themselves are not connected to other edges |
|---|
| 515 | |
|---|
| 516 | For a continuous riverwall, this should only occur twice (at the end points) |
|---|
| 517 | |
|---|
| 518 | Otherwise, the riverwall appears discontinous |
|---|
| 519 | """ |
|---|
| 520 | |
|---|
| 521 | domain=self.domain |
|---|
| 522 | |
|---|
| 523 | if(len(self.names)==0): |
|---|
| 524 | print 'There are no riverwalls' |
|---|
| 525 | return |
|---|
| 526 | |
|---|
| 527 | # Shorthand notation |
|---|
| 528 | rwd=self |
|---|
| 529 | |
|---|
| 530 | for i, name in enumerate(rwd.names): |
|---|
| 531 | # Get indices of edges on this riverwall |
|---|
| 532 | riverwalledgeInds=rwd.riverwall_edges[(rwd.hydraulic_properties_rowIndex==i).nonzero()[0]] |
|---|
| 533 | # Get their corresponding centroids |
|---|
| 534 | riverwallCentInds=rwd.get_centroids_corresponding_to_edgeInds(riverwalledgeInds) |
|---|
| 535 | # Get their corresponding vertices |
|---|
| 536 | riverwallV1Inds, riverwallV2Inds = rwd.get_vertices_corresponding_to_edgeInds(riverwalledgeInds, tol=tol) |
|---|
| 537 | |
|---|
| 538 | # With discontinuous triangles, we expect edges to occur twice |
|---|
| 539 | # Let's remove duplicates to simplify the analysis |
|---|
| 540 | repeat=riverwalledgeInds*0 |
|---|
| 541 | lre=len(riverwalledgeInds) |
|---|
| 542 | # edge coordinates as a complex number, for easy equality checking |
|---|
| 543 | complex_edge_coordinates=domain.edge_coordinates[riverwalledgeInds,0]+\ |
|---|
| 544 | 1j*domain.edge_coordinates[riverwalledgeInds,1] |
|---|
| 545 | for j in range(lre-1): |
|---|
| 546 | # Ignore if already checked |
|---|
| 547 | if(repeat[j]==1): |
|---|
| 548 | continue |
|---|
| 549 | # Check for a dupulicate |
|---|
| 550 | dups=(complex_edge_coordinates[(j+1):lre]==complex_edge_coordinates[j]).nonzero()[0] |
|---|
| 551 | if(len(dups)>0): |
|---|
| 552 | repeat[dups+j+1]=1 |
|---|
| 553 | |
|---|
| 554 | unique_riverwall_edge_indices=(repeat==0).nonzero()[0] |
|---|
| 555 | |
|---|
| 556 | # Finally, get 'unqiue' edges in the riverwall |
|---|
| 557 | uEdges=riverwalledgeInds[unique_riverwall_edge_indices] |
|---|
| 558 | #uCentroids=riverwallCentInds[unique_riverwall_edge_indices] |
|---|
| 559 | uV1=riverwallV1Inds[unique_riverwall_edge_indices] |
|---|
| 560 | uV2=riverwallV2Inds[unique_riverwall_edge_indices] |
|---|
| 561 | |
|---|
| 562 | # Next, count how many times each vertex value occurs. |
|---|
| 563 | # For a 'connected' riverwall, we only want 2 edges where a vertex occurs only once. |
|---|
| 564 | lure=len(uEdges) |
|---|
| 565 | complex_v1_coordinates=domain.vertex_coordinates[uV1,0]+\ |
|---|
| 566 | 1j*domain.vertex_coordinates[uV1,1] |
|---|
| 567 | complex_v2_coordinates=domain.vertex_coordinates[uV2,0]+\ |
|---|
| 568 | 1j*domain.vertex_coordinates[uV2,1] |
|---|
| 569 | v1Counter=uEdges*0 |
|---|
| 570 | v2Counter=uEdges*0 |
|---|
| 571 | for j in range(lure): |
|---|
| 572 | v1Counter[j]=(complex_v1_coordinates==complex_v1_coordinates[j]).sum()+\ |
|---|
| 573 | (complex_v2_coordinates==complex_v1_coordinates[j]).sum() |
|---|
| 574 | v2Counter[j]=(complex_v1_coordinates==complex_v2_coordinates[j]).sum()+\ |
|---|
| 575 | (complex_v2_coordinates==complex_v2_coordinates[j]).sum() |
|---|
| 576 | |
|---|
| 577 | num_disconnected_edges=(v1Counter==1).sum()+(v2Counter==1).sum() |
|---|
| 578 | |
|---|
| 579 | if(verbose): |
|---|
| 580 | print 'There are ', num_disconnected_edges, ' vertices with only 1 connection on riverwall ', name |
|---|
| 581 | |
|---|
| 582 | if(num_disconnected_edges==0): |
|---|
| 583 | print 'Riverwall '+name+' is not on this mesh [typical for a sub-mesh in parallel]' |
|---|
| 584 | elif(num_disconnected_edges==2): |
|---|
| 585 | if(verbose): |
|---|
| 586 | print " That's Good, it means the wall is continuous" |
|---|
| 587 | print ' ' |
|---|
| 588 | else: |
|---|
| 589 | msg='Riverwall ' + name +', appears to be discontinuous, possibly due to round-off issues,'+\ |
|---|
| 590 | ' be careful [try adjusting tol, or use warn_only=True to avoid error] ' |
|---|
| 591 | if(warn_only): |
|---|
| 592 | print msg |
|---|
| 593 | else: |
|---|
| 594 | raise Exception, msg |
|---|
| 595 | |
|---|
| 596 | return |
|---|
| 597 | |
|---|
