Changeset 8009


Ignore:
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:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/__init__.py

    r8006 r8009  
    166166                               interior_regions=None,
    167167                               interior_holes=None,
    168                                                            hole_tags = None,
     168                               hole_tags=None,
    169169                               poly_geo_reference=None,
    170170                               mesh_geo_reference=None,
     
    235235              'interior_regions': interior_regions,
    236236              'interior_holes': interior_holes,
    237                           'hole_tags': hole_tags,
     237              'hole_tags': hole_tags,
    238238              'poly_geo_reference': poly_geo_reference,
    239239              'mesh_geo_reference': mesh_geo_reference,
     
    269269                                interior_regions=None,
    270270                                interior_holes=None,
    271                                                                 hole_tags=None,
     271                                hole_tags=None,
    272272                                poly_geo_reference=None,
    273273                                mesh_geo_reference=None,
     
    289289                             filename=mesh_filename,
    290290                             interior_holes=interior_holes,
    291                                                         hole_tags=hole_tags,
     291                            hole_tags=hole_tags,
    292292                             poly_geo_reference=poly_geo_reference,
    293293                             mesh_geo_reference=mesh_geo_reference,
  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r7317 r8009  
    202202
    203203    return points, elements, boundary
     204
     205
     206
     207def rectangular_cross_slit(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
     208
     209    """Setup a rectangular grid of triangles
     210    with m+1 by n+1 grid points
     211    and side lengths len1, len2. If side lengths are omitted
     212    the mesh defaults to the unit square.
     213
     214    len1: x direction (left to right)
     215    len2: y direction (bottom to top)
     216
     217    Return to lists: points and elements suitable for creating a Mesh or
     218    FVMesh object, e.g. Mesh(points, elements)
     219    """
     220
     221    from anuga.config import epsilon
     222
     223    delta1 = float(len1)/m
     224    delta2 = float(len2)/n
     225
     226    #Dictionary of vertex objects
     227    vertices = {}
     228    points = []
     229
     230    for i in range(m+1):
     231        for j in range(n+1):
     232            vertices[i,j] = len(points)
     233            points.append([delta1*i + origin[0], delta2*j  + origin[1]])
     234
     235    # Construct 4 triangles per element
     236    elements = []
     237    boundary = {}
     238    for i in range(m):
     239        for j in range(n):
     240            v1 = vertices[i,j+1]
     241            v2 = vertices[i,j]
     242            v3 = vertices[i+1,j+1]
     243            v4 = vertices[i+1,j]
     244            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
     245            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
     246
     247            # Create centre point
     248            v5 = len(points)
     249            points.append([x, y])
     250
     251            #Create left triangle
     252            if i == 0:
     253                boundary[(len(elements), 1)] = 'left'
     254            elements.append([v2,v5,v1])
     255
     256            #Create bottom triangle
     257            if j == 0:
     258                boundary[(len(elements), 1)] = 'bottom'
     259            elements.append([v4,v5,v2])
     260
     261            #Create right triangle
     262            if i == m-1:
     263                boundary[(len(elements), 1)] = 'right'
     264            elements.append([v3,v5,v4])
     265
     266            #Create top triangle
     267            if j == n-1:
     268                boundary[(len(elements), 1)] = 'top'
     269            elements.append([v1,v5,v3])
     270
     271
     272    return points, elements, boundary
     273
    204274
    205275
  • trunk/anuga_core/source/anuga/geometry/polygon.py

    r7858 r8009  
    599599        assert polygon.shape[1] == 2, msg
    600600
     601
    601602        msg = ('Points array must be 1 or 2 dimensional. '
    602603               'I got %d dimensions' % len(points.shape))
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r7952 r8009  
    587587        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    588588
     589
     590        //Update momentum
     591        xmom_update[k] += S*uh[k];
     592        ymom_update[k] += S*vh[k];
     593      }
     594    }
     595  }
     596}
     597
     598
     599
     600void _chezy_friction(double g, double eps, int N,
     601                           double* x, double* w, double* zv,
     602                           double* uh, double* vh,
     603                           double* chezy, double* xmom_update, double* ymom_update) {
     604
     605  int k, k3, k6;
     606  double S, h, z, z0, z1, z2, zs, zx, zy;
     607  double x0,y0,x1,y1,x2,y2;
     608
     609  for (k=0; k<N; k++) {
     610    if (chezy[k] > eps) {
     611      k3 = 3*k;
     612      // Get bathymetry
     613      z0 = zv[k3 + 0];
     614      z1 = zv[k3 + 1];
     615      z2 = zv[k3 + 2];
     616
     617      // Compute bed slope
     618      k6 = 6*k;  // base index
     619
     620      x0 = x[k6 + 0];
     621      y0 = x[k6 + 1];
     622      x1 = x[k6 + 2];
     623      y1 = x[k6 + 3];
     624      x2 = x[k6 + 4];
     625      y2 = x[k6 + 5];
     626
     627      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
     628
     629      zs = sqrt(1.0 + zx*zx + zy*zy);
     630      z = (z0+z1+z2)/3.0;
     631      h = w[k]-z;
     632      if (h >= eps) {
     633        S = -g * chezy[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     634        S /= pow(h, 2.0);
    589635
    590636        //Update momentum
     
    11311177                        (double*) vh -> data,
    11321178                        (double*) eta -> data,
     1179                        (double*) xmom -> data,
     1180                        (double*) ymom -> data);
     1181
     1182  return Py_BuildValue("");
     1183}
     1184
     1185
     1186PyObject *chezy_friction(PyObject *self, PyObject *args) {
     1187  //
     1188  // chezy_friction(g, eps, x, h, uh, vh, z, chezy, xmom_update, ymom_update)
     1189  //
     1190
     1191
     1192  PyArrayObject *x, *w, *z, *uh, *vh, *chezy, *xmom, *ymom;
     1193  int N;
     1194  double g, eps;
     1195
     1196  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
     1197                        &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
     1198    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: chezy_friction could not parse input arguments");
     1199    return NULL;
     1200  }
     1201
     1202  // check that numpy array objects arrays are C contiguous memory
     1203  CHECK_C_CONTIG(x);
     1204  CHECK_C_CONTIG(w);
     1205  CHECK_C_CONTIG(z);
     1206  CHECK_C_CONTIG(uh);
     1207  CHECK_C_CONTIG(vh);
     1208  CHECK_C_CONTIG(chezy);
     1209  CHECK_C_CONTIG(xmom);
     1210  CHECK_C_CONTIG(ymom);
     1211
     1212  N = w -> dimensions[0];
     1213
     1214  _chezy_friction(g, eps, N,
     1215                        (double*) x -> data,
     1216                        (double*) w -> data,
     1217                        (double*) z -> data,
     1218                        (double*) uh -> data,
     1219                        (double*) vh -> data,
     1220                        (double*) chezy -> data,
    11331221                        (double*) xmom -> data,
    11341222                        (double*) ymom -> data);
     
    28182906  {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"},
    28192907  {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"},
     2908  {"chezy_friction", chezy_friction, METH_VARARGS, "Print out"},
    28202909  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
    28212910  {"balance_deep_and_shallow", balance_deep_and_shallow,
Note: See TracChangeset for help on using the changeset viewer.