streamFunction.patch

gistfile1.diff

diff -ru streamFunction/Make/files streamFunctionV/Make/files
--- streamFunction/Make/files	2013-06-11 11:33:58.000000000 +0200
+++ streamFunctionV/Make/files	2013-12-17 15:15:36.255538814 +0100
@@ -1,3 +1,3 @@
 streamFunction.C
 
-EXE = $(FOAM_APPBIN)/streamFunction
+EXE = $(FOAM_USER_APPBIN)/streamFunctionV
diff -ru streamFunction/streamFunction.C streamFunctionV/streamFunction.C
--- streamFunction/streamFunction.C	2013-06-11 11:33:58.000000000 +0200
+++ streamFunctionV/streamFunction.C	2013-12-16 16:57:13.830205945 +0100
@@ -100,6 +100,34 @@
                 dimensionedScalar("zero", phi.dimensions(), 0.0)
             );
 
+            volScalarField streamFunctionV
+            (
+                IOobject
+                (
+                    "streamFunctionV",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar("zero", phi.dimensions(), 0.0)
+            );
+
+            volScalarField weight
+            (
+                IOobject
+                (
+                    "weight",
+                    runTime.timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimless
+            );
+
             labelList visitedPoint(mesh.nPoints());
             forAll(visitedPoint, pointI)
             {
@@ -449,6 +477,24 @@
                 Info<< endl;
             } while (!finished);
 
+            streamFunctionV = dimensionedScalar("zero", phi.dimensions(), 0.0);
+            weight = 0.0;
+            const labelListList& cellPoints = streamFunctionV.mesh().cellPoints();
+            forAll(streamFunctionV, cellI)
+            {
+                const labelList& curCellPoints = cellPoints[cellI];
+                forAll(curCellPoints, cellPointI)
+                {
+                    streamFunctionV[cellI] += 
+                        streamFunction[curCellPoints[cellPointI]];
+                    weight[cellI] += 1.0;
+                }
+            }
+
+            streamFunctionV /= (weight + 1.0);
+            streamFunctionV.boundaryField() = 0.0;
+            streamFunctionV.write();
+
             streamFunction.boundaryField() = 0.0;
             streamFunction.write();
         }

stream-function-mpl.py

def plot_stream_function_mpl(data, filename=None, time=0):
    """Plot stream function isolines over velocity field."""
    pp.clf()
    xes, ys, u, v, mag = list2table(data[0])
    pp.streamplot(xes/CM, ys/CM, u, v, density=2, linewidth=1, color='white',
                  arrowstyle='fancy')
    vd = data[0]
    pp.tricontourf(vd[:, 0]/CM, vd[:, 1]/CM,
                   np.sqrt(vd[:, 3]**2 + vd[:, 4]**2), 16)

    x_ = data[2][:, 0]
    y_ = data[2][:, 1]
    fl_ = data[2][:, 3]
    pp.tricontour(x_/CM, y_/CM, fl_, [0.5], colors='#d3d7cf',
                  linewidths=1)

    pp.title(r'Stream function, t = %d s' % time)
    pp.gca().set_aspect('equal', adjustable='box', anchor='C')
    pp.xlim([0, 10])
    pp.ylim([0, 10])
    pp.xlabel('x, cm')
    pp.ylabel('y, cm')

    setup_axes()

    if filename:
        pp.savefig(filename)

stream-function.py

def plot_stream_function(data, filename=None, time=0):
    """Plot stream function isolines over velocity field."""

    pp.clf()

    vd = data[0]
    pp.tricontourf(vd[:, 0]/CM, vd[:, 1]/CM,
                   np.sqrt(vd[:, 3]**2 + vd[:, 4]**2), 16)

    sfd = data[3]
    below = np.linspace(min(sfd[:, 3]), -1e-32, 32)
    above = np.linspace(1e-32, max(sfd[:, 3]), 32)
    pp.tricontour(sfd[:, 0]/CM, sfd[:, 1]/CM, sfd[:, 3], below + above,
                  linewidths=1, colors='white', linestyles='solid')

    x_ = data[2][:, 0]
    y_ = data[2][:, 1]
    fl_ = data[2][:, 3]
    pp.tricontour(x_/CM, y_/CM, fl_, [0.5], colors='#d3d7cf',
                  linewidths=1)

    pp.title(r'Stream function, t = %d s' % time)
    pp.gca().set_aspect('equal', adjustable='box', anchor='C')
    pp.xlim([0, 10])
    pp.ylim([0, 10])
    pp.xlabel('x, cm')
    pp.ylabel('y, cm')

    setup_axes()

    if filename:
        pp.savefig(filename)