# Changeset 1387

Ignore:
Timestamp:
May 13, 2005, 6:15:08 PM (19 years ago)
Message:

Need to look at uint test for inscribed circle

Location:
inundation/ga/storm_surge
Files:
17 edited

### Legend:

Unmodified
Removed

 r1363 domain.visualise = True #domain.visualise_range_z = 1.0 domain.visualise_color_stage = True domain.visualise_range_z = 0.5 #domain.visualise_color_stage = True

 r1363 from config import g, epsilon from Numeric import allclose, array, zeros, ones, Float from advection import * from parallel_advection import * from Numeric import array from mesh_factory import rectangular points, vertices, boundary = rectangular(60, 60) points, vertices, boundary = rectangular(30, 30) #Create advection domain with direction (1,-1) domain = Domain(points, vertices, boundary, velocity=[1.0, 0.0]) domain = Parallel_Domain(points, vertices, boundary, velocity=[1.0, 0.0]) # Initial condition is zero by default domain.set_boundary( {'left': D, 'right': T, 'bottom': T, 'top': T} ) domain.set_boundary( {'left': T, 'right': T, 'bottom': T, 'top': T} ) domain.check_integrity() class Set_Stage: """Set an initial condition with constant water height, for xself.x0)&(x
• ## inundation/ga/storm_surge/parallel/small.tsh

 r1363 31 0 # <# of verts> <# of vert attributes>, next lines [attributes] ...Triangulation Vertices... 0 554.0 -242.0 1 443.0 -190.0 2 340.0 -184.0 3 276.0 -237.0 4 281.0 -347.0 5 282.0 -380.0 6 312.0 -425.0 7 387.0 -472.0 8 465.0 -487.0 9 549.0 -469.0 10 566.0 -413.0 11 589.0 -359.0 12 472.0 -258.0 13 377.0 -411.0 14 428.0 -257.0 15 386.0 -278.0 16 372.0 -325.0 17 369.0 -371.0 18 433.0 -422.0 19 476.0 -361.0 20 502.0 -304.0 21 333.0 -225.0 22 305.0 -244.0 23 333.0 -270.0 24 369.0 -237.0 25 278.5 -292.0 26 499.90712817 -423.516449623 27 532.621862616 -366.885237781 28 498.5 -216.0 29 322.359287054 -318.872505543 30 571.5 -300.5 0 554.0 -242.0 1 443.0 -190.0 2 340.0 -184.0 3 276.0 -237.0 4 281.0 -347.0 5 282.0 -380.0 6 312.0 -425.0 7 387.0 -472.0 8 465.0 -487.0 9 549.0 -469.0 10 566.0 -413.0 11 589.0 -359.0 12 472.0 -258.0 13 377.0 -411.0 14 428.0 -257.0 15 386.0 -278.0 16 372.0 -325.0 17 369.0 -371.0 18 433.0 -422.0 19 476.0 -361.0 20 502.0 -304.0 21 333.0 -225.0 22 305.0 -244.0 23 333.0 -270.0 24 369.0 -237.0 25 278.5 -292.0 26 499.90712817 -423.516449623 27 532.621862616 -366.885237781 28 498.5 -216.0 29 322.359287054 -318.872505543 30 571.5 -300.5 # attribute column titles ...Triangulation Vertex Titles... 38 # <# of triangles>, next lines [] [] [attribute of region] ...Triangulation Triangles... 0 29 5 17 1 5 26 1 17 5 6 -1 3 0 2 10 11 27 28 21 -1 3 13 17 6 1 6 -1 4 15 23 16 24 -1 16 5 16 29 17 0 -1 24 6 6 7 13 37 3 -1 7 29 23 25 9 8 24 8 25 4 29 26 7 -1 9 25 23 22 10 18 7 10 22 23 21 13 11 9 11 3 22 21 10 15 18 12 21 24 2 14 15 13 13 24 21 23 10 16 12 14 2 24 1 27 -1 12 15 21 2 3 -1 11 12 16 24 23 15 4 30 13 17 27 30 20 31 36 28 18 22 3 25 -1 9 11 19 7 8 18 35 37 -1 20 8 9 26 23 35 -1 21 26 10 27 2 25 23 22 26 19 18 -1 35 25 23 10 26 9 20 -1 21 24 16 23 29 7 5 4 25 26 27 19 36 22 21 26 29 4 5 -1 0 8 27 14 1 24 14 30 29 28 27 11 30 -1 17 2 29 1 14 12 -1 33 27 30 14 24 15 16 -1 27 31 20 30 0 -1 34 17 32 12 0 28 -1 33 34 33 12 28 1 -1 29 32 34 0 12 20 -1 31 32 35 8 26 18 22 19 20 36 27 20 19 -1 25 17 37 18 13 7 6 19 -1 0 29 5 17 1 5 26 building 1 17 5 6 -1 3 0 2 10 11 27 28 21 -1 3 13 17 6 1 6 -1 4 15 23 16 24 -1 16 5 16 29 17 0 -1 24 6 6 7 13 37 3 -1 7 29 23 25 9 8 24 8 25 4 29 26 7 -1 9 25 23 22 10 18 7 10 22 23 21 13 11 9 building 11 3 22 21 10 15 18 12 21 24 2 14 15 13 13 24 21 23 10 16 12 14 2 24 1 27 -1 12 15 21 2 3 -1 11 12 16 24 23 15 4 30 13 17 27 30 20 31 36 28 house 18 22 3 25 -1 9 11 19 7 8 18 35 37 -1 20 8 9 26 23 35 -1 21 26 10 27 2 25 23 22 26 19 18 -1 35 25 building 23 10 26 9 20 -1 21 24 16 23 29 7 5 4 25 26 27 19 36 22 21 26 29 4 5 -1 0 8 27 14 1 24 14 30 29 house 28 27 11 30 -1 17 2 29 1 14 12 -1 33 27 oval 30 14 24 15 16 -1 27 31 20 30 0 -1 34 17 32 12 0 28 -1 33 34 33 12 28 1 -1 29 32 34 0 12 20 -1 31 32 35 8 26 18 22 19 20 36 27 20 19 -1 25 17 37 18 13 7 6 19 -1 28 # <# of segments>, next lines   [boundary tag] ...Triangulation Segments... 0 14 12 exterior 7 19 20 exterior 8 12 20 exterior 9 21 22 10 23 22 11 23 24 12 21 24 9 21 22 10 23 22 11 23 24 12 21 24 13 4 25 exterior 14 7 6 exterior 27 30 11 exterior 25 0 # <# of verts> <# of vert attributes>, next lines [attributes] ...Mesh Vertices... 0 554.0 -242.0 1 443.0 -190.0 2 340.0 -184.0 3 276.0 -237.0 4 281.0 -347.0 5 282.0 -380.0 6 312.0 -425.0 7 387.0 -472.0 8 465.0 -487.0 9 549.0 -469.0 10 566.0 -413.0 11 589.0 -359.0 12 472.0 -258.0 13 377.0 -411.0 14 428.0 -257.0 15 386.0 -278.0 16 372.0 -325.0 17 369.0 -371.0 18 433.0 -422.0 19 476.0 -361.0 20 502.0 -304.0 21 333.0 -225.0 22 305.0 -244.0 23 333.0 -270.0 24 369.0 -237.0 0 554.0 -242.0 1 443.0 -190.0 2 340.0 -184.0 3 276.0 -237.0 4 281.0 -347.0 5 282.0 -380.0 6 312.0 -425.0 7 387.0 -472.0 8 465.0 -487.0 9 549.0 -469.0 10 566.0 -413.0 11 589.0 -359.0 12 472.0 -258.0 13 377.0 -411.0 14 428.0 -257.0 15 386.0 -278.0 16 372.0 -325.0 17 369.0 -371.0 18 433.0 -422.0 19 476.0 -361.0 20 502.0 -304.0 21 333.0 -225.0 22 305.0 -244.0 23 333.0 -270.0 24 369.0 -237.0 25 # <# of segments>, next lines   [boundary tag] ...Mesh Segments... 0 14 12 1 15 14 2 16 15 3 17 16 4 13 17 5 18 13 6 19 18 7 20 19 8 12 20 9 22 21 10 23 22 11 24 23 12 21 24 13 3 4 14 6 7 15 5 6 16 2 3 17 4 5 18 7 8 19 9 10 20 8 9 21 0 1 22 1 2 23 11 0 24 10 11 0 14 12 1 15 14 2 16 15 3 17 16 4 13 17 5 18 13 6 19 18 7 20 19 8 12 20 9 22 21 10 23 22 11 24 23 12 21 24 13 3 4 14 6 7 15 5 6 16 2 3 17 4 5 18 7 8 19 9 10 20 8 9 21 0 1 22 1 2 23 11 0 24 10 11 1 # <# of holes>, next lines ...Mesh Holes... 0 418.0 -332.0 1 # <# of regions>, next lines ...Mesh Regions... 0 334.0 -243.0 0 334.0 -243.0 1 # <# of regions>, next lines [Max Area] ...Mesh Regions... 0 0 # Geo reference 56
• ## inundation/ga/storm_surge/pyvolution/domain.py

 r1360 x = self.boundary.keys() x.sort() for k, (vol_id, edge_id) in enumerate(x): tag = self.boundary[ (vol_id, edge_id) ]
• ## inundation/ga/storm_surge/pyvolution/general_mesh.py

 r1178 A triangular element is defined in terms of three vertex ids, ordered counter clock-wise, ordered counter clock-wise, each corresponding to a given coordinate set. Vertices from different elements can point to the same Coordinate sets are implemented as an N x 2 Numeric array containing x and y coordinates. x and y coordinates. To instantiate: c = [2.0,0.0] e = [2.0, 2.0] points = [a, b, c, e] triangles = [ [1,0,2], [1,2,3] ]   #bac, bce triangles = [ [1,0,2], [1,2,3] ]   #bac, bce mesh = Mesh(points, triangles) #creates two triangles: bac and bce Other: In addition mesh computes an Nx6 array called vertex_coordinates. This structure is derived from coordinates and contains for each This structure is derived from coordinates and contains for each triangle the three x,y coordinates at the vertices. This is a cut down version of mesh from pyvolution\mesh.py This is a cut down version of mesh from pyvolution mesh.py """ origin is a 3-tuple consisting of UTM zone, easting and northing. If specified coordinates are assumed to be relative to this origin. If specified coordinates are assumed to be relative to this origin. """ if geo_reference is None: self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) else: self.geo_reference = geo_reference #Input checks msg = 'Triangles must an Nx2 Numeric array or a sequence of 2-tuples' msg = 'Triangles must an Nx3 Numeric array or a sequence of 3-tuples' assert len(self.triangles.shape) == 2, msg msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples' assert len(self.coordinates.shape) == 2, msg assert len(self.coordinates.shape) == 2, msg msg = 'Vertex indices reference non-existing coordinate sets' #Register number of elements (N) self.number_of_elements = N = self.triangles.shape[0] #Allocate space for geometric quantities #Allocate space for geometric quantities self.normals = zeros((N, 6), Float) self.areas = zeros(N, Float) #Get x,y coordinates for all triangles and store self.vertex_coordinates = V = self.compute_vertex_coordinates() #Initialise each triangle for i in range(N): #if i % (N/10) == 0: print '(%d/%d)' %(i, N) x0 = V[i, 0]; y0 = V[i, 1] x1 = V[i, 2]; y1 = V[i, 3] x2 = V[i, 4]; y2 = V[i, 5] x2 = V[i, 4]; y2 = V[i, 5] #Area msg += ' is degenerate:  area == %f' %self.areas[i] assert self.areas[i] > 0.0, msg #Normals n0 = array([x2 - x1, y2 - y1]) l0 = sqrt(sum(n0**2)) l0 = sqrt(sum(n0**2)) n1 = array([x0 - x2, y0 - y2]) n1[1], -n1[0], n2[1], -n2[0]] #Edgelengths self.edgelengths[i, :] = [l0, l1, l2] #Build vertex list self.build_vertexlist() self.build_vertexlist() def __len__(self): return self.number_of_elements """Return all normal vectors. Return normal vectors for all triangles as an Nx6 array (ordered as x0, y0, x1, y1, x2, y2 for each triangle) (ordered as x0, y0, x1, y1, x2, y2 for each triangle) """ return self.normals Return value is the numeric array slice [x, y] """ return self.normals[i, 2*j:2*j+2] return self.normals[i, 2*j:2*j+2] def get_vertex_coordinates(self): """Return all vertex coordinates. Return all vertex coordinates for all triangles as an Nx6 array (ordered as x0, y0, x1, y1, x2, y2 for each triangle) (ordered as x0, y0, x1, y1, x2, y2 for each triangle) """ return self.vertex_coordinates Return value is the numeric array slice [x, y] """ return self.vertex_coordinates[i, 2*j:2*j+2] return self.vertex_coordinates[i, 2*j:2*j+2] #FIXME: Perhaps they should be ordered as in obj files?? #See quantity.get_vertex_values from Numeric import zeros, Float N = self.number_of_elements vertex_coordinates = zeros((N, 6), Float) vertex_coordinates = zeros((N, 6), Float) for i in range(N): return vertex_coordinates def get_vertices(self,  indexes=None): """Get connectivity indexes is the set of element ids of interest """ from Numeric import take if (indexes ==  None): indexes = range(len(self))  #len(self)=number of elements return  take(self.triangles, indexes) unique_verts = {} for triangle in triangles: unique_verts[triangle[0]] = 0 unique_verts[triangle[1]] = 0 unique_verts[triangle[0]] = 0 unique_verts[triangle[1]] = 0 unique_verts[triangle[2]] = 0 return unique_verts.keys() return unique_verts.keys() def build_vertexlist(self): """Build vertexlist index by vertex ids and for each entry (point id) Preconditions: self.coordinates and self.triangles are defined self.coordinates and self.triangles are defined Postcondition: self.vertexlist is built a = self.triangles[i, 0] b = self.triangles[i, 1] c = self.triangles[i, 2] c = self.triangles[i, 2] #Register the vertices v as lists of vertexlist[v] = [] vertexlist[v].append( (i, vertex_id) ) vertexlist[v].append( (i, vertex_id) ) self.vertexlist = vertexlist def get_extent(self): X = C[:,0:6:2].copy() Y = C[:,1:6:2].copy() xmin = min(X.flat) xmax = max(X.flat) ymin = min(Y.flat) ymax = max(Y.flat) ymax = max(Y.flat) return xmin, xmax, ymin, ymax """ return sum(self.areas) return sum(self.areas)

• ## inundation/ga/storm_surge/pyvolution/netherlands.py

 r1360 # from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ Transmissive_boundary, Constant_height Transmissive_boundary, Constant_height, Constant_stage from mesh_factory import rectangular N = 130 #size = 33800 N = 600 #Size = 720000 N = 60 N = 50 M = 40 #N = 15 print 'Creating domain' #Create basic mesh points, vertices, boundary = rectangular(N, N) points, vertices, boundary = rectangular(N, M) #Create shallow water domain #Set bed-slope and friction inflow_stage = 0.1 manning = 0.02 inflow_stage = 0.2 manning = 0.00 Z = Weir(inflow_stage) # print 'Initial condition' domain.set_quantity('stage', Constant_height(Z, 0.)) domain.set_quantity('stage', Constant_height(Z, 0.0)) #domain.set_quantity('stage', Constant_stage(-1.0)) #Evolve t0 = time.time() for t in domain.evolve(yieldstep = 0.05, finaltime = 10.0): for t in domain.evolve(yieldstep = 0.01, finaltime = 5.0): domain.write_time() #domain.visualiser.scale_z = 1.0 domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0)
• ## inundation/ga/storm_surge/pyvolution/realtime_visualisation_new.py

 r1363 elevation_color = (0.3,0.4,0.5) stage_color = (0.1,0.4,0.99) other_color = (0.99,0.4,0.1) class Visualiser: def update_quantity(self,qname,qcolor=None): def update_quantity(self,qname,qcolor=None,scale_z=None): #print 'update '+qname+' arrays' if qcolor is None: qcolor = other_color if qname=='elevation': qcolor = elevation_color qcolor = stage_color if scale_z is None: scale_z = self.scale_z try: self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor) self.update_arrays(self.domain.quantities[qname].vertex_values, qcolor, scale_z) #print 'update bed image' pass def update_quantity_color(self,qname,qcolor=None): def update_quantity_color(self,qname,qcolor=None,scale_z=None): #print 'update '+qname+' arrays' qcolor = stage_color if scale_z is None: scale_z = self.scale_z #try: self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor) self.update_arrays_color(self.domain.quantities[qname].vertex_values, qcolor, scale_z) #print 'update bed image' def update_arrays(self,quantity,qcolor): def update_arrays(self,quantity,qcolor,scale_z): col = asarray(qcolor) N = len(self.domain) scale_z = self.scale_z min_x = self.min_x min_y = self.min_y def update_arrays_color(self,quantity,qcolor): def update_arrays_color(self,quantity,qcolor,scale_z): col = asarray(qcolor) N = len(self.domain) scale_z = self.scale_z min_x = self.min_x min_y = self.min_y normal(1) = normal(1)/norm; normal(2) = normal(2)/norm; for (int j=0; j<3; j++) {
• ## inundation/ga/storm_surge/pyvolution/shallow_water.py

 r1363 class Constant_stage: """Set an initial condition with constant stage """ def __init__(self, s): self.s = s def __call__(self, x, y): return self.s class Constant_height: """Set an initial condition with constant water height, e.g stage s = z+h """ def __init__(self, W, h): self.W = W
• ## inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

 r1017 // // To compile (Python2.3): //  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O //  gcc -shared domain_ext.o  -o domain_ext.so //  gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O //  gcc -shared domain_ext.o  -o domain_ext.so // // or use python compile.py // // Ole Nielsen, GA 2004 #include "Python.h" #include "Numeric/arrayobject.h" #include "math.h" #include //Shared code snippets /*Rotate the momentum component q (q[1], q[2]) from x,y coordinates to coordinates based on normal vector (n1, n2). Result is returned in array 3x1 r Result is returned in array 3x1 r To rotate in opposite direction, call rotate with (q, n1, -n2) Contents of q are changed by this function */ Contents of q are changed by this function */ double q1, q2; //Shorthands q1 = q[1];  //uh momentum //Rotate q[1] =  n1*q1 + n2*q2; q[2] = -n2*q1 + n1*q2; q[2] = -n2*q1 + n1*q2; return 0; } } // Computational function for flux computation (using stage w=z+h) int flux_function(double *q_left, double *q_right, double z_left, double z_right, double n1, double n2, double epsilon, double g, int flux_function(double *q_left, double *q_right, double z_left, double z_right, double n1, double n2, double epsilon, double g, double *edgeflux, double *max_speed) { /*Compute fluxes between volumes for the shallow water wave equation cast in terms of the 'stage', w = h+z using cast in terms of the 'stage', w = h+z using the 'central scheme' as described in Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'. Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740. The implemented formula is given in equation (3.15) on page 714 The implemented formula is given in equation (3.15) on page 714 */ int i; double w_left, h_left, uh_left, vh_left, u_left; double w_right, h_right, uh_right, vh_right, u_right; double s_min, s_max, soundspeed_left, soundspeed_right; double denom, z; double q_left_copy[3], q_right_copy[3]; double q_left_copy[3], q_right_copy[3]; double flux_right[3], flux_left[3]; //Copy conserved quantities to protect from modification for (i=0; i<3; i++) { q_left_copy[i] = q_left[i]; q_right_copy[i] = q_right[i]; } } //Align x- and y-momentum with x-axis _rotate(q_left_copy, n1, n2); _rotate(q_right_copy, n1, n2); _rotate(q_right_copy, n1, n2); z = (z_left+z_right)/2; //Take average of field values h_left = 0.0;  //Could have been negative u_left = 0.0; } else { } else { u_left = uh_left/h_left; } w_right = q_right_copy[0]; h_right = w_right-z; h_right = 0.0; //Could have been negative u_right = 0.0; } else { } else { u_right = uh_right/h_right; } //Momentum in y-direction //Momentum in y-direction vh_left  = q_left_copy[2]; vh_right = q_right_copy[2]; vh_right = q_right_copy[2]; //Maximal and minimal wave speeds soundspeed_left  = sqrt(g*h_left); soundspeed_left  = sqrt(g*h_left); soundspeed_right = sqrt(g*h_right); s_max = max(u_left+soundspeed_left, u_right+soundspeed_right); if (s_max < 0.0) s_max = 0.0; if (s_max < 0.0) s_max = 0.0; s_min = min(u_left-soundspeed_left, u_right-soundspeed_right); if (s_min > 0.0) s_min = 0.0; //Flux formulas if (s_min > 0.0) s_min = 0.0; //Flux formulas flux_left[0] = u_left*h_left; flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left; flux_left[2] = u_left*vh_left; flux_right[0] = u_right*h_right; flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right; flux_right[2] = u_right*vh_right; //Flux computation //Flux computation denom = s_max-s_min; if (denom == 0.0) { for (i=0; i<3; i++) edgeflux[i] = 0.0; *max_speed = 0.0; } else { } else { for (i=0; i<3; i++) { edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i]; edgeflux[i] += s_max*s_min*(q_right_copy[i]-q_left_copy[i]); edgeflux[i] /= denom; } } //Maximal wavespeed *max_speed = max(fabs(s_max), fabs(s_min)); //Rotate back //Rotate back _rotate(edgeflux, n1, -n2); } return 0; } void _manning_friction(double g, double eps, int N, double* w, double* z, double* uh, double* vh, double* eta, double* xmom, double* ymom) { void _manning_friction(double g, double eps, int N, double* w, double* z, double* uh, double* vh, double* eta, double* xmom, double* ymom) { int k; double S, h; for (k=0; k eps) { } } int _balance_deep_and_shallow(int N, double* wc, double* zc, double* hc, double* wv, double* zv, double* zc, double* hc, double* wv, double* zv, double* hv, double* hvbar, double* xmomc, double* ymomc, double* xmomv, double* ymomv) { double* hvbar, double* xmomc, double* ymomc, double* xmomv, double* ymomv) { int k, k3, i; double dz, hmin, alpha; //Compute linear combination between w-limited stages and //h-limited stages close to the bed elevation. //h-limited stages close to the bed elevation. for (k=0; k dz/2 then alpha = 1 and the bed will have no effect //If hmin < 0 then alpha = 0 reverting to constant height above bed. if (dz > 0.0) //if (hmin<0.0) //if (hmin<0.0) //        alpha = 0.0; //else //  alpha = max( min( hc[k]/dz, 1.0), 0.0 ); alpha = max( min( 2.0*hmin/dz, 1.0), 0.0 ); else alpha = 1.0;  //Flat bed //alpha = 1.0; //printf("dz = %.3f, alpha = %.8f\n", dz, alpha); //  Let // //    wvi be the w-limited stage (wvi = zvi + hvi) //  Let // //    wvi be the w-limited stage (wvi = zvi + hvi) //    wvi- be the h-limited state (wvi- = zvi + hvi-) // // // // //  where i=0,1,2 denotes the vertex ids // //  Weighted balance between w-limited and h-limited stage is // //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi) // // //  Weighted balance between w-limited and h-limited stage is // //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi) // //  It follows that the updated wvi is //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi // // //   Momentum is balanced between constant and limited if (alpha < 1) { if (alpha < 1) { for (i=0; i<3; i++) { wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i]; //Update momentum as a linear combination of //Update momentum as a linear combination of //xmomc and ymomc (shallow) and momentum //from extrapolator xmomv and ymomv (deep). xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i]; xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i]; ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i]; } } } } } return 0; } int _protect(int N, double minimum_allowed_height, int _protect(int N, double minimum_allowed_height, double* wc, double* zc, double* xmomc, double* zc, double* xmomc, double* ymomc) { int k; int k; double hc; //Protect against initesimal and negative heights for (k=0; k dimensions[0]; for (k=0; k data)[k3+i]; } } avg_h /= 3; //Compute bed slope x0 = ((double*) x -> data)[k6 + 0]; y0 = ((double*) x -> data)[k6 + 1]; y0 = ((double*) x -> data)[k6 + 1]; x1 = ((double*) x -> data)[k6 + 2]; y1 = ((double*) x -> data)[k6 + 3]; y1 = ((double*) x -> data)[k6 + 3]; x2 = ((double*) x -> data)[k6 + 4]; y2 = ((double*) x -> data)[k6 + 5]; y2 = ((double*) x -> data)[k6 + 5]; z0 = ((double*) v -> data)[k3 + 0]; z1 = ((double*) v -> data)[k3 + 1]; z2 = ((double*) v -> data)[k3 + 2]; z2 = ((double*) v -> data)[k3 + 2]; _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy); //Update momentum ((double*) xmom -> data)[k] += -g*zx*avg_h; ((double*) ymom -> data)[k] += -g*zy*avg_h; } ((double*) ymom -> data)[k] += -g*zy*avg_h; } return Py_BuildValue(""); } // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update) // PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom; int N; double g, eps; if (!PyArg_ParseTuple(args, "ddOOOOOOO", &g, &eps, &w, &z, &uh, &vh, &eta, &xmom, &ymom)) return NULL; N = w -> dimensions[0]; &g, &eps, &w, &z, &uh, &vh, &eta, &xmom, &ymom)) return NULL; N = w -> dimensions[0]; _manning_friction(g, eps, N, (double*) w -> data, (double*) z -> data, (double*) uh -> data, (double*) vh -> data, (double*) z -> data, (double*) uh -> data, (double*) vh -> data, (double*) eta -> data, (double*) xmom -> data, (double*) xmom -> data, (double*) ymom -> data); return Py_BuildValue(""); } } PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) { // normal a Float numeric array of length 2. PyObject *Q, *Normal; PyArrayObject *q, *r, *normal; static char *argnames[] = {"q", "normal", "direction", NULL}; int dimensions[1], i, direction=1; double n1, n2; // Convert Python arguments to C if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, &Q, &Normal, &direction)) return NULL; // Convert Python arguments to C if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames, &Q, &Normal, &direction)) return NULL; //Input checks (convert sequences into numeric arrays) q = (PyArrayObject *) q = (PyArrayObject *) PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0); normal = (PyArrayObject *) normal = (PyArrayObject *) PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0); if (normal -> dimensions[0] != 2) { PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components"); } //Allocate space for return vector r (don't DECREF) //Allocate space for return vector r (don't DECREF) dimensions[0] = 3; r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); //Copy for (i=0; i<3; i++) { ((double *) (r -> data))[i] = ((double *) (q -> data))[i]; } ((double *) (r -> data))[i] = ((double *) (q -> data))[i]; } //Get normal and direction n1 = ((double *) normal -> data)[0]; n2 = ((double *) normal -> data)[1]; n1 = ((double *) normal -> data)[0]; n2 = ((double *) normal -> data)[1]; if (direction == -1) n2 = -n2; //Release numeric arrays Py_DECREF(q); Py_DECREF(q); Py_DECREF(normal); //return result using PyArray to avoid memory leak return PyArray_Return(r); } } Fluxes across each edge are scaled by edgelengths and summed up Resulting flux is then scaled by area and stored in explicit_update for each of the three conserved quantities explicit_update for each of the three conserved quantities stage, xmomentum and ymomentum is converted to a timestep that must not be exceeded. The minimum of those is computed as the next overall timestep. Python call: domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g, domain.neighbours, domain.neighbour_edges, domain.normals, domain.edgelengths, domain.edgelengths, domain.radii, domain.areas, Stage.edge_values, Xmom.edge_values, Ymom.edge_values, Bed.edge_values, Stage.edge_values, Xmom.edge_values, Ymom.edge_values, Bed.edge_values, Stage.boundary_values, Xmom.boundary_values, Xmom.explicit_update, Ymom.explicit_update) Post conditions: domain.explicit_update is reset to computed flux values domain.timestep is set to the largest step satisfying all volumes. domain.timestep is set to the largest step satisfying all volumes. */ PyArrayObject *neighbours, *neighbour_edges, *normals, *edgelengths, *radii, *areas, *stage_edge_values, *xmom_edge_values, *ymom_edge_values, *bed_edge_values, *stage_edge_values, *xmom_edge_values, *ymom_edge_values, *bed_edge_values, *stage_boundary_values, *xmom_boundary_values, *ymom_explicit_update; //Local variables //Local variables double timestep, max_speed, epsilon, g; double normal[2], ql[3], qr[3], zl, zr; int number_of_elements, k, i, j, m, n; int ki, nm, ki2; //Index shorthands // Convert Python arguments to C // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO", ×tep, &epsilon, &g, &neighbours, &neighbours, &neighbour_edges, &normals, &normals, &edgelengths, &radii, &areas, &stage_edge_values, &xmom_edge_values, &ymom_edge_values, &bed_edge_values, &stage_edge_values, &xmom_edge_values, &ymom_edge_values, &bed_edge_values, &stage_boundary_values, &xmom_boundary_values, number_of_elements = stage_edge_values -> dimensions[0]; for (k=0; k data)[ki]; ql[1] = ((double *) xmom_edge_values -> data)[ki]; ql[2] = ((double *) ymom_edge_values -> data)[ki]; zl =    ((double *) bed_edge_values -> data)[ki]; ql[2] = ((double *) ymom_edge_values -> data)[ki]; zl =    ((double *) bed_edge_values -> data)[ki]; //Quantities at neighbour on nearest face if (n < 0) { m = -n-1; //Convert negative flag to index qr[0] = ((double *) stage_boundary_values -> data)[m]; qr[1] = ((double *) xmom_boundary_values -> data)[m]; qr[2] = ((double *) ymom_boundary_values -> data)[m]; qr[0] = ((double *) stage_boundary_values -> data)[m]; qr[1] = ((double *) xmom_boundary_values -> data)[m]; qr[2] = ((double *) ymom_boundary_values -> data)[m]; zr = zl; //Extend bed elevation to boundary } else { } else { m = ((long *) neighbour_edges -> data)[ki]; nm = n*3+m; nm = n*3+m; qr[0] = ((double *) stage_edge_values -> data)[nm]; qr[1] = ((double *) xmom_edge_values -> data)[nm]; qr[2] = ((double *) ymom_edge_values -> data)[nm]; zr =    ((double *) bed_edge_values -> data)[nm]; qr[2] = ((double *) ymom_edge_values -> data)[nm]; zr =    ((double *) bed_edge_values -> data)[nm]; } //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m, //             ql[0], ql[1], ql[2]); // Outward pointing normal vector //             ql[0], ql[1], ql[2]); // Outward pointing normal vector // normal = domain.normals[k, 2*i:2*i+2] ki2 = 2*ki; //k*6 + i*2 normal[0] = ((double *) normals -> data)[ki2]; normal[1] = ((double *) normals -> data)[ki2+1]; normal[1] = ((double *) normals -> data)[ki2+1]; //Edge flux computation flux_function(ql, qr, zl, zr, flux_function(ql, qr, zl, zr, normal[0], normal[1], epsilon, g, epsilon, g, edgeflux, &max_speed); //flux -= edgeflux * edgelengths[k,i] for (j=0; j<3; j++) { for (j=0; j<3; j++) { flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki]; } //Update timestep //timestep = min(timestep, domain.radii[k]/max_speed) if (max_speed > epsilon) { timestep = min(timestep, ((double *) radii -> data)[k]/max_speed); } } } // end for i //Normalise by area and store for when all conserved //quantities get updated // flux /= areas[k] for (j=0; j<3; j++) { for (j=0; j<3; j++) { flux[j] /= ((double *) areas -> data)[k]; } ((double *) stage_explicit_update -> data)[k] = flux[0]; ((double *) xmom_explicit_update -> data)[k] = flux[1]; ((double *) ymom_explicit_update -> data)[k] = flux[2]; ((double *) ymom_explicit_update -> data)[k] = flux[2]; } //end for k return Py_BuildValue("d", timestep); } } // //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc) PyArrayObject PyArrayObject *wc,            //Stage at centroids *zc,            //Elevation at centroids *zc,            //Elevation at centroids *xmomc,         //Momentums at centroids *ymomc; int N; *ymomc; int N; double minimum_allowed_height; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dOOOO", // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dOOOO", &minimum_allowed_height, &wc, &zc, &xmomc, &ymomc)) N = wc -> dimensions[0]; _protect(N, minimum_allowed_height, (double*) wc -> data, (double*) zc -> data, (double*) xmomc -> data, (double*) zc -> data, (double*) xmomc -> data, (double*) ymomc -> data); return Py_BuildValue(""); return Py_BuildValue(""); } //    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, //                             xmomc, ymomc, xmomv, ymomv) PyArrayObject PyArrayObject *wc,            //Stage at centroids *zc,            //Elevation at centroids *hc,            //Height at centroids *zc,            //Elevation at centroids *hc,            //Height at centroids *wv,            //Stage at vertices *zv,            //Elevation at vertices *hv,            //Depths at vertices *hvbar,         //h-Limited depths at vertices *hv,            //Depths at vertices *hvbar,         //h-Limited depths at vertices *xmomc,         //Momentums at centroids and vertices *ymomc, *xmomv, *ymomv; *ymomc, *xmomv, *ymomv; int N; //, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOOOOOOOOOO", &wc, &zc, &hc, // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOOOOOOOOOO", &wc, &zc, &hc, &wv, &zv, &hv, &hvbar, &xmomc, &ymomc, &xmomv, &ymomv)) N = wc -> dimensions[0]; _balance_deep_and_shallow(N, (double*) wc -> data, (double*) zc -> data, (double*) hc -> data, (double*) wv -> data, (double*) zv -> data, (double*) zc -> data, (double*) hc -> data, (double*) wv -> data, (double*) zv -> data, (double*) hv -> data, (double*) hvbar -> data, (double*) xmomc -> data, (double*) ymomc -> data, (double*) xmomv -> data, (double*) ymomv -> data); return Py_BuildValue(""); (double*) hvbar -> data, (double*) xmomc -> data, (double*) ymomc -> data, (double*) xmomv -> data, (double*) ymomv -> data); return Py_BuildValue(""); } PyObject *h_limiter(PyObject *self, PyObject *args) { PyObject *domain, *Tmp; PyArrayObject *hv, *hc, //Depth at vertices and centroids PyArrayObject *hv, *hc, //Depth at vertices and centroids *hvbar,   //Limited depth at vertices (return values) *neighbours; int k, i, n, N, k3; int dimensions[2]; double beta_h; //Safety factor (see config.py) double *hmin, *hmax, hn; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) return NULL; //Get safety factor beta_h Tmp = PyObject_GetAttrString(domain, "beta_h"); if (!Tmp) return NULL; beta_h = PyFloat_AsDouble(Tmp); Py_DECREF(Tmp); Tmp = PyObject_GetAttrString(domain, "beta_h"); if (!Tmp) return NULL; beta_h = PyFloat_AsDouble(Tmp); Py_DECREF(Tmp); N = hc -> dimensions[0]; //Create hvbar dimensions[0] = N; dimensions[1] = 3; dimensions[1] = 3; hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE); //Find min and max of this and neighbour's centroid values hmin = malloc(N * sizeof(double)); hmax = malloc(N * sizeof(double)); hmax = malloc(N * sizeof(double)); for (k=0; k data)[k]; hmax[k] = hmin[k]; for (i=0; i<3; i++) { n = ((long*) neighbours -> data)[k3+i]; //Initialise hvbar with values from hv ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i]; ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i]; if (n >= 0) { hn = ((double*) hc -> data)[n]; //Neighbour's centroid value hmin[k] = min(hmin[k], hn); hmax[k] = max(hmax[k], hn); } } // Call underlying standard routine _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax); // // //Py_DECREF(domain); //FIXME: NEcessary? // // //Py_DECREF(domain); //FIXME: NEcessary? free(hmin); free(hmax); //return result using PyArray to avoid memory leak free(hmax); //return result using PyArray to avoid memory leak return PyArray_Return(hvbar); //return Py_BuildValue(""); //return Py_BuildValue(""); } //                              s_vec, phi_vec, self.const) PyArrayObject   //(one element per triangle) *s_vec,         //Speeds *phi_vec,       //Bearings *s_vec,         //Speeds *phi_vec,       //Bearings *xmom_update,   //Momentum updates *ymom_update; int N; *ymom_update; int N; double cw; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOOOd", // Convert Python arguments to C if (!PyArg_ParseTuple(args, "OOOOd", &xmom_update, &ymom_update, &s_vec, &phi_vec, &ymom_update, &s_vec, &phi_vec, &cw)) return NULL; N = xmom_update -> dimensions[0]; _assign_wind_field_values(N, (double*) xmom_update -> data, (double*) ymom_update -> data, (double*) s_vec -> data, (double*) ymom_update -> data, (double*) s_vec -> data, (double*) phi_vec -> data, cw); return Py_BuildValue(""); } ////////////////////////////////////////// return Py_BuildValue(""); } ////////////////////////////////////////// // Method table for python module static struct PyMethodDef MethodTable[] = { * three. */ {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"}, {"gravity", gravity, METH_VARARGS, "Print out"}, {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, {"balance_deep_and_shallow", balance_deep_and_shallow, METH_VARARGS, "Print out"}, {"h_limiter", h_limiter, METH_VARARGS, "Print out"}, {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"}, {"assign_windfield_values", assign_windfield_values, METH_VARARGS | METH_KEYWORDS, "Print out"}, //{"distribute_to_vertices_and_edges", // distribute_to_vertices_and_edges, METH_VARARGS}, //{"update_conserved_quantities", // update_conserved_quantities, METH_VARARGS}, //{"set_initialcondition", // set_initialcondition, METH_VARARGS}, {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"}, {"gravity", gravity, METH_VARARGS, "Print out"}, {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, {"balance_deep_and_shallow", balance_deep_and_shallow, METH_VARARGS, "Print out"}, {"h_limiter", h_limiter, METH_VARARGS, "Print out"}, {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"}, {"assign_windfield_values", assign_windfield_values, METH_VARARGS | METH_KEYWORDS, "Print out"}, //{"distribute_to_vertices_and_edges", // distribute_to_vertices_and_edges, METH_VARARGS}, //{"update_conserved_quantities", // update_conserved_quantities, METH_VARARGS}, //{"set_initialcondition", // set_initialcondition, METH_VARARGS}, {NULL, NULL} }; // Module initialisation // Module initialisation void initshallow_water_ext(void){ Py_InitModule("shallow_water_ext", MethodTable); import_array();     //Necessary for handling of NumPY structures } import_array();     //Necessary for handling of NumPY structures }