Changeset 8009

Timestamp:

Sep 8, 2010, 5:12:10 PM (14 years ago)

Author:

steve

Message:

Slight conflict with init.py due to addition of hole_tags

Location:

trunk/anuga_core/source/anuga

Files:

• __init__.py (modified) (4 diffs)
• abstract_2d_finite_volumes/mesh_factory.py (modified) (1 diff)
• geometry/polygon.py (modified) (1 diff)
• shallow_water/shallow_water_ext.c (modified) (3 diffs)

trunk/anuga_core/source/anuga/__init__.py

@@ -166 +166 @@
                                interior_regions=None,
                                interior_holes=None,
-                                                           hole_tags = None,
+                               hole_tags=None,
                                poly_geo_reference=None,
                                mesh_geo_reference=None,
@@ -235 +235 @@
               'interior_regions': interior_regions,
               'interior_holes': interior_holes,
-                          'hole_tags': hole_tags,
+              'hole_tags': hole_tags,
               'poly_geo_reference': poly_geo_reference,
               'mesh_geo_reference': mesh_geo_reference,
@@ -269 +269 @@
                                 interior_regions=None,
                                 interior_holes=None,
-                                                                hole_tags=None,
+                                hole_tags=None,
                                 poly_geo_reference=None,
                                 mesh_geo_reference=None,
@@ -289 +289 @@
                              filename=mesh_filename,
                              interior_holes=interior_holes,
-                                                         hole_tags=hole_tags,
+                             hole_tags=hole_tags,
                              poly_geo_reference=poly_geo_reference,
                              mesh_geo_reference=mesh_geo_reference,

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

@@ -202 +202 @@
 
     return points, elements, boundary
+
+
+
+def rectangular_cross_slit(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+    """Setup a rectangular grid of triangles
+    with m+1 by n+1 grid points
+    and side lengths len1, len2. If side lengths are omitted
+    the mesh defaults to the unit square.
+
+    len1: x direction (left to right)
+    len2: y direction (bottom to top)
+
+    Return to lists: points and elements suitable for creating a Mesh or
+    FVMesh object, e.g. Mesh(points, elements)
+    """
+
+    from anuga.config import epsilon
+
+    delta1 = float(len1)/m
+    delta2 = float(len2)/n
+
+    #Dictionary of vertex objects
+    vertices = {}
+    points = []
+
+    for i in range(m+1):
+        for j in range(n+1):
+            vertices[i,j] = len(points)
+            points.append([delta1*i + origin[0], delta2*j  + origin[1]])
+
+    # Construct 4 triangles per element
+    elements = []
+    boundary = {}
+    for i in range(m):
+        for j in range(n):
+            v1 = vertices[i,j+1]
+            v2 = vertices[i,j]
+            v3 = vertices[i+1,j+1]
+            v4 = vertices[i+1,j]
+            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
+            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
+
+            # Create centre point
+            v5 = len(points)
+            points.append([x, y])
+
+            #Create left triangle
+            if i == 0:
+                boundary[(len(elements), 1)] = 'left'
+            elements.append([v2,v5,v1])
+
+            #Create bottom triangle
+            if j == 0:
+                boundary[(len(elements), 1)] = 'bottom'
+            elements.append([v4,v5,v2])
+
+            #Create right triangle
+            if i == m-1:
+                boundary[(len(elements), 1)] = 'right'
+            elements.append([v3,v5,v4])
+
+            #Create top triangle
+            if j == n-1:
+                boundary[(len(elements), 1)] = 'top'
+            elements.append([v1,v5,v3])
+
+
+    return points, elements, boundary
+
 
 

trunk/anuga_core/source/anuga/geometry/polygon.py

@@ -599 +599 @@
         assert polygon.shape[1] == 2, msg
 
+
         msg = ('Points array must be 1 or 2 dimensional. '
                'I got %d dimensions' % len(points.shape))

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -587 +587 @@
         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
 
+
+        //Update momentum
+        xmom_update[k] += S*uh[k];
+        ymom_update[k] += S*vh[k];
+      }
+    }
+  }
+}
+
+
+
+void _chezy_friction(double g, double eps, int N,
+                           double* x, double* w, double* zv,
+                           double* uh, double* vh,
+                           double* chezy, double* xmom_update, double* ymom_update) {
+
+  int k, k3, k6;
+  double S, h, z, z0, z1, z2, zs, zx, zy;
+  double x0,y0,x1,y1,x2,y2;
+
+  for (k=0; k<N; k++) {
+    if (chezy[k] > eps) {
+      k3 = 3*k;
+      // Get bathymetry
+      z0 = zv[k3 + 0];
+      z1 = zv[k3 + 1];
+      z2 = zv[k3 + 2];
+
+      // Compute bed slope
+      k6 = 6*k;  // base index
+
+      x0 = x[k6 + 0];
+      y0 = x[k6 + 1];
+      x1 = x[k6 + 2];
+      y1 = x[k6 + 3];
+      x2 = x[k6 + 4];
+      y2 = x[k6 + 5];
+
+      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
+
+      zs = sqrt(1.0 + zx*zx + zy*zy);
+      z = (z0+z1+z2)/3.0;
+      h = w[k]-z;
+      if (h >= eps) {
+        S = -g * chezy[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
+        S /= pow(h, 2.0);
 
         //Update momentum
@@ -1131 +1177 @@
                         (double*) vh -> data,
                         (double*) eta -> data,
+                        (double*) xmom -> data,
+                        (double*) ymom -> data);
+
+  return Py_BuildValue("");
+}
+
+
+PyObject *chezy_friction(PyObject *self, PyObject *args) {
+  //
+  // chezy_friction(g, eps, x, h, uh, vh, z, chezy, xmom_update, ymom_update)
+  //
+
+
+  PyArrayObject *x, *w, *z, *uh, *vh, *chezy, *xmom, *ymom;
+  int N;
+  double g, eps;
+
+  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
+                        &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: chezy_friction could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(x);
+  CHECK_C_CONTIG(w);
+  CHECK_C_CONTIG(z);
+  CHECK_C_CONTIG(uh);
+  CHECK_C_CONTIG(vh);
+  CHECK_C_CONTIG(chezy);
+  CHECK_C_CONTIG(xmom);
+  CHECK_C_CONTIG(ymom);
+
+  N = w -> dimensions[0];
+
+  _chezy_friction(g, eps, N,
+                        (double*) x -> data,
+                        (double*) w -> data,
+                        (double*) z -> data,
+                        (double*) uh -> data,
+                        (double*) vh -> data,
+                        (double*) chezy -> data,
                         (double*) xmom -> data,
                         (double*) ymom -> data);
@@ -2818 +2906 @@
   {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"},
   {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"},
+  {"chezy_friction", chezy_friction, METH_VARARGS, "Print out"},
   {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
   {"balance_deep_and_shallow", balance_deep_and_shallow,