| 1 | """ |
|---|
| 2 | Set elevation operators |
|---|
| 3 | |
|---|
| 4 | |
|---|
| 5 | """ |
|---|
| 6 | |
|---|
| 7 | __author__="steve" |
|---|
| 8 | __date__ ="$09/03/2012 4:46:39 PM$" |
|---|
| 9 | |
|---|
| 10 | |
|---|
| 11 | from anuga import Domain |
|---|
| 12 | from anuga import Quantity |
|---|
| 13 | import numpy as num |
|---|
| 14 | import anuga.utilities.log as log |
|---|
| 15 | |
|---|
| 16 | from anuga.geometry.polygon import inside_polygon |
|---|
| 17 | |
|---|
| 18 | from anuga.operators.set_quantity import Set_quantity |
|---|
| 19 | from anuga.config import indent |
|---|
| 20 | |
|---|
| 21 | class Set_elevation(Set_quantity): |
|---|
| 22 | """ |
|---|
| 23 | Helper class to setup calculation of elevation |
|---|
| 24 | associated with a region (defined by indices, polygon or center/radius |
|---|
| 25 | """ |
|---|
| 26 | |
|---|
| 27 | set_elevation = Set_quantity.set_value |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | def __init__(self, |
|---|
| 31 | domain, |
|---|
| 32 | elevation=None, |
|---|
| 33 | indices=None, |
|---|
| 34 | polygon=None, |
|---|
| 35 | center=None, |
|---|
| 36 | radius=None, |
|---|
| 37 | line=None, |
|---|
| 38 | verbose = False): |
|---|
| 39 | |
|---|
| 40 | Set_quantity.__init__(self, domain, 'elevation', |
|---|
| 41 | value=elevation, |
|---|
| 42 | indices=indices, |
|---|
| 43 | polygon=polygon, |
|---|
| 44 | center=center, |
|---|
| 45 | radius=radius, |
|---|
| 46 | line=line, |
|---|
| 47 | verbose=verbose, |
|---|
| 48 | test_elevation=False) |
|---|
| 49 | |
|---|
| 50 | |
|---|
| 51 | |
|---|
| 52 | #----------------------------------------- |
|---|
| 53 | # Extra structures to support maintaining |
|---|
| 54 | # continuity of elevation |
|---|
| 55 | #----------------------------------------- |
|---|
| 56 | self.setup_node_structures() |
|---|
| 57 | |
|---|
| 58 | |
|---|
| 59 | #------------------------------------------ |
|---|
| 60 | # Extra aliases for changing elevation |
|---|
| 61 | # to ensure conservation of mass and |
|---|
| 62 | # continuity at vertices and edges |
|---|
| 63 | #------------------------------------------ |
|---|
| 64 | self.elev_v = self.domain.quantities['elevation'].vertex_values |
|---|
| 65 | self.elev_e = self.domain.quantities['elevation'].edge_values |
|---|
| 66 | self.stage_c = self.domain.quantities['stage'].centroid_values |
|---|
| 67 | self.elev_c = self.domain.quantities['elevation'].centroid_values |
|---|
| 68 | |
|---|
| 69 | #------------------------------------------ |
|---|
| 70 | # x,y coordinates of vertices of cells that are |
|---|
| 71 | # updated |
|---|
| 72 | #------------------------------------------ |
|---|
| 73 | coord_v = self.domain.vertex_coordinates |
|---|
| 74 | ids = self.indices |
|---|
| 75 | if ids is None: |
|---|
| 76 | self.v_x = coord_v[:,0].reshape((-1,3)) |
|---|
| 77 | self.v_y = coord_v[:,1].reshape((-1,3)) |
|---|
| 78 | else: |
|---|
| 79 | self.v_x = coord_v[:,0].reshape((-1,3))[ids] |
|---|
| 80 | self.v_y = coord_v[:,1].reshape((-1,3))[ids] |
|---|
| 81 | |
|---|
| 82 | #print self.v_x.shape |
|---|
| 83 | #print self.v_y.shape |
|---|
| 84 | |
|---|
| 85 | #------------------------------------------ |
|---|
| 86 | # Need to turn off this optimization as it |
|---|
| 87 | # doesn't fixup the relationship between |
|---|
| 88 | # bed and stage vertex values in dry region |
|---|
| 89 | #------------------------------------------ |
|---|
| 90 | self.domain.optimise_dry_cells = 0 |
|---|
| 91 | |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | def __call__(self): |
|---|
| 96 | """ |
|---|
| 97 | Apply rate to those triangles defined in indices |
|---|
| 98 | |
|---|
| 99 | indices == [], then don't apply anywhere |
|---|
| 100 | indices == None, then apply everywhere |
|---|
| 101 | otherwise apply for the specific indices |
|---|
| 102 | """ |
|---|
| 103 | |
|---|
| 104 | if self.value is None: |
|---|
| 105 | return |
|---|
| 106 | |
|---|
| 107 | if self.indices is []: |
|---|
| 108 | return |
|---|
| 109 | |
|---|
| 110 | #------------------------------------------ |
|---|
| 111 | # Apply changes to elevation vertex values |
|---|
| 112 | # via the update_quantites routine |
|---|
| 113 | # Assume vertex values updated and need to |
|---|
| 114 | # fix up centroid values unless |
|---|
| 115 | # domain.get_discontinuous_elevation is true |
|---|
| 116 | #------------------------------------------ |
|---|
| 117 | if not self.update_quantities(): |
|---|
| 118 | return |
|---|
| 119 | |
|---|
| 120 | |
|---|
| 121 | #------------------------------------------ |
|---|
| 122 | # Cleanup elevation and stage quantity values |
|---|
| 123 | #----------------------------------------- |
|---|
| 124 | if self.indices is None: |
|---|
| 125 | |
|---|
| 126 | #-------------------------------------- |
|---|
| 127 | # Make elevation continuous and clean up |
|---|
| 128 | # stage values to ensure conservation |
|---|
| 129 | #-------------------------------------- |
|---|
| 130 | if self.domain.get_using_discontinuous_elevation(): |
|---|
| 131 | pass |
|---|
| 132 | else: |
|---|
| 133 | height_c = self.stage_c - self.elev_c |
|---|
| 134 | self.domain.quantities['elevation'].smooth_vertex_values() |
|---|
| 135 | self.domain.quantities['elevation'].interpolate() |
|---|
| 136 | self.stage_c[:] = self.elev_c + height_c |
|---|
| 137 | |
|---|
| 138 | |
|---|
| 139 | else: |
|---|
| 140 | |
|---|
| 141 | #-------------------------------------- |
|---|
| 142 | # Make elevation continuous and clean up |
|---|
| 143 | # stage values to ensure conservation |
|---|
| 144 | #-------------------------------------- |
|---|
| 145 | |
|---|
| 146 | if self.domain.get_using_discontinuous_elevation(): |
|---|
| 147 | pass |
|---|
| 148 | else: |
|---|
| 149 | height_c = self.stage_c[self.vols] - self.elev_c[self.vols] |
|---|
| 150 | for nid in self.node_ids: |
|---|
| 151 | non = self.domain.number_of_triangles_per_node[nid] |
|---|
| 152 | |
|---|
| 153 | vid = num.arange(self.node_index[nid], self.node_index[nid+1],dtype=num.int) |
|---|
| 154 | vidd = self.domain.vertex_value_indices[vid] |
|---|
| 155 | |
|---|
| 156 | self.elev_v[vidd/3,vidd%3] = num.sum(self.elev_v[vidd/3,vidd%3])/non |
|---|
| 157 | |
|---|
| 158 | #-------------------------------------- |
|---|
| 159 | # clean up the centroid values and edge values |
|---|
| 160 | #-------------------------------------- |
|---|
| 161 | self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1) |
|---|
| 162 | |
|---|
| 163 | self.elev_e[self.vols,0] = 0.5*(self.elev_v[self.vols,1]+ self.elev_v[self.vols,2]) |
|---|
| 164 | self.elev_e[self.vols,1] = 0.5*(self.elev_v[self.vols,2]+ self.elev_v[self.vols,0]) |
|---|
| 165 | self.elev_e[self.vols,2] = 0.5*(self.elev_v[self.vols,0]+ self.elev_v[self.vols,1]) |
|---|
| 166 | |
|---|
| 167 | self.stage_c[self.vols] = self.elev_c[self.vols] + height_c |
|---|
| 168 | |
|---|
| 169 | |
|---|
| 170 | |
|---|
| 171 | def update_quantities(self): |
|---|
| 172 | """Update the vertex values of the quantities to model erosion |
|---|
| 173 | """ |
|---|
| 174 | |
|---|
| 175 | if self.value is None: |
|---|
| 176 | return False |
|---|
| 177 | |
|---|
| 178 | |
|---|
| 179 | updated = True |
|---|
| 180 | |
|---|
| 181 | if self.indices is None: |
|---|
| 182 | |
|---|
| 183 | if self.domain.get_using_discontinuous_elevation(): |
|---|
| 184 | try: |
|---|
| 185 | height_c = self.stage_c - self.elev_c |
|---|
| 186 | value = self.get_value(x=self.coord_c[:,0], y=self.coord_c[:,1]) |
|---|
| 187 | self.elev_c[:] = value |
|---|
| 188 | self.stage_c[:] = self.elev_c + height_c |
|---|
| 189 | except ValueError: |
|---|
| 190 | updated = False |
|---|
| 191 | pass |
|---|
| 192 | else: |
|---|
| 193 | #-------------------------------------- |
|---|
| 194 | # Update all three vertices for each cell |
|---|
| 195 | #-------------------------------------- |
|---|
| 196 | try: |
|---|
| 197 | value = self.get_value(self.v_x, self.v_y) |
|---|
| 198 | self.elev_v[:] = value |
|---|
| 199 | except ValueError: |
|---|
| 200 | updated = False |
|---|
| 201 | pass |
|---|
| 202 | |
|---|
| 203 | #---------------------------------- |
|---|
| 204 | # Apply just to indices |
|---|
| 205 | #---------------------------------- |
|---|
| 206 | else: |
|---|
| 207 | |
|---|
| 208 | if self.domain.get_using_discontinuous_elevation(): |
|---|
| 209 | ids = self.indices |
|---|
| 210 | x = self.coord_c[ids,0] |
|---|
| 211 | y = self.coord_c[ids,1] |
|---|
| 212 | try: |
|---|
| 213 | height_c = self.stage_c[ids] - self.elev_c[ids] |
|---|
| 214 | value = self.get_value(x=x,y=y) |
|---|
| 215 | self.elev_c[ids] = value |
|---|
| 216 | self.stage_c[ids] = self.elev_c[ids] + height_c |
|---|
| 217 | except ValueError: |
|---|
| 218 | updated = False |
|---|
| 219 | pass |
|---|
| 220 | else: |
|---|
| 221 | ids = self.indices |
|---|
| 222 | try: |
|---|
| 223 | value = self.get_value(self.v_x, self.v_y) |
|---|
| 224 | self.elev_v[ids] = value |
|---|
| 225 | except ValueError: |
|---|
| 226 | updated = False |
|---|
| 227 | pass |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | return updated |
|---|
| 231 | |
|---|
| 232 | |
|---|
| 233 | |
|---|
| 234 | |
|---|
| 235 | def setup_node_structures(self): |
|---|
| 236 | """ For setting elevation we need to |
|---|
| 237 | ensure that the elevation quantity remains |
|---|
| 238 | continuous (at least for version 1.3 of anuga) |
|---|
| 239 | |
|---|
| 240 | So we need to find all the vertices that need to |
|---|
| 241 | update within each timestep. |
|---|
| 242 | """ |
|---|
| 243 | |
|---|
| 244 | |
|---|
| 245 | |
|---|
| 246 | if self.indices is None or self.indices is []: |
|---|
| 247 | self.vol_ids = None |
|---|
| 248 | self.vols = None |
|---|
| 249 | self.vert_ids = None |
|---|
| 250 | return |
|---|
| 251 | |
|---|
| 252 | |
|---|
| 253 | |
|---|
| 254 | node_ids = set() |
|---|
| 255 | |
|---|
| 256 | for ind in self.indices: |
|---|
| 257 | for k in [0,1,2]: |
|---|
| 258 | node_ids.add(self.domain.triangles[ind,k]) |
|---|
| 259 | |
|---|
| 260 | self.node_ids = [ id for id in node_ids ] |
|---|
| 261 | |
|---|
| 262 | |
|---|
| 263 | |
|---|
| 264 | node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int) |
|---|
| 265 | |
|---|
| 266 | # FIXME: SR Don't we calculate this for the domain already! |
|---|
| 267 | k = 0 |
|---|
| 268 | node_index[0] = 0 |
|---|
| 269 | for i in range(self.domain.number_of_nodes): |
|---|
| 270 | node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i] |
|---|
| 271 | |
|---|
| 272 | self.node_index = node_index |
|---|
| 273 | |
|---|
| 274 | vertex_ids =[] |
|---|
| 275 | for nid in self.node_ids: |
|---|
| 276 | #print nid,self.domain.number_of_triangles_per_node[nid] |
|---|
| 277 | for vid in range(node_index[nid], node_index[nid+1]): |
|---|
| 278 | vidd = self.domain.vertex_value_indices[vid] |
|---|
| 279 | vertex_ids.append(vidd) |
|---|
| 280 | #print ' ',nid, vid, vidd, vidd/3, vidd%3 |
|---|
| 281 | |
|---|
| 282 | self.vol_ids = num.array(vertex_ids,dtype=num.int)/3 |
|---|
| 283 | self.vols = num.array(list(set(self.vol_ids)), dtype=num.int) |
|---|
| 284 | self.vert_ids = num.array(vertex_ids,dtype=num.int)%3 |
|---|
| 285 | |
|---|
| 286 | |
|---|
| 287 | |
|---|
| 288 | |
|---|
| 289 | |
|---|
| 290 | |
|---|
