Mercurial > hg > cfcfd3

changeset 2251:27ff0e3d15f6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Eilmer4: flow derivative code for 3D. It's in place but untested in this revision. Also, made the 2D derivative code work from the flow points stored in FVVertex objects. Once the flow points are stored in the vertex objects, the flow-property-derivative function has been greatly simplified because every vertex can be treated the same.
author Peter Jacobs <p.jacobs@uq.edu.au>
date Sat, 06 Jun 2015 21:15:20 +1000
parents 314b46321471
children c9cebf43e8d5
files dlang/eilmer4/block.d dlang/eilmer4/fvvertex.d dlang/eilmer4/sblock.d dlang/eilmer4/simcore.d dlang/eilmer4/viscousflux.d
diffstat 5 files changed, 682 insertions(+), 123 deletions(-) [+]
line wrap: on
line diff
--- a/dlang/eilmer4/block.d	Sat Jun 06 10:15:41 2015 +1000
+++ b/dlang/eilmer4/block.d	Sat Jun 06 21:15:20 2015 +1000
@@ -69,6 +69,7 @@
     abstract int count_invalid_cells(int gtl);
     abstract void compute_primary_cell_geometric_data(int gtl);
     abstract void compute_distance_to_nearest_wall_for_all_cells(int gtl);
+    abstract void assign_flow_locations_for_derivative_calc(size_t gtl);
     abstract void read_grid(string filename, size_t gtl=0);
     abstract void write_grid(string filename, double sim_time, size_t gtl=0);
     abstract double read_solution(string filename);
--- a/dlang/eilmer4/fvvertex.d	Sat Jun 06 10:15:41 2015 +1000
+++ b/dlang/eilmer4/fvvertex.d	Sat Jun 06 21:15:20 2015 +1000
@@ -13,6 +13,7 @@
 import fvcore;
 import geom;
 import gas;
+import flowstate;
 
 class FVVertex {
 public:
@@ -30,6 +31,8 @@
     Vector3 grad_omega;  // pseudo vorticity for k-omega turbulence
     Vector3[] grad_f;    // mass fraction derivatives
     Vector3 grad_pe;     // electron pressure derivatives
+    Vector3[] cloud_pos; // Positions of flow points for derivative calculation.
+    FlowState[] cloud_fs; // References to flow states at those points.
 
     this(in GasModel gm, size_t id_init=0)
     {
--- a/dlang/eilmer4/sblock.d	Sat Jun 06 10:15:41 2015 +1000
+++ b/dlang/eilmer4/sblock.d	Sat Jun 06 21:15:20 2015 +1000
@@ -1016,6 +1016,497 @@
 	}
     } // end calc_ghost_cell_geom_2D()
 
+    override void assign_flow_locations_for_derivative_calc(size_t gtl)
+    {
+	size_t i, j, k;
+	if (myConfig.dimensions == 2) {
+	    // First, do all of the internal secondary cells.
+	    // i.e. Those not on a boundary.
+	    for ( i = imin+1; i <= imax; ++i ) {
+		for ( j = jmin+1; j <= jmax; ++j ) {
+		    // Secondary-cell centre is a primary-cell vertex.
+		    FVVertex vtx = get_vtx(i,j);
+		    // These are the corners of the secondary cell.
+		    FVCell A = get_cell(i,j-1);
+		    FVCell B = get_cell(i,j);
+		    FVCell C = get_cell(i-1,j);
+		    FVCell D = get_cell(i-1,j-1);
+		    // Retain locations and references to flow states for later.
+		    vtx.cloud_pos = [A.pos[gtl], B.pos[gtl], C.pos[gtl], D.pos[gtl]];
+		    vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+		} // j loop
+	    } // i loop
+	    // Half-cells along the edges of the block.
+	    // East boundary
+	    i = imax+1;
+	    for (j = jmin+1; j <= jmax; ++j) {
+		FVVertex vtx = get_vtx(i,j);
+		FVInterface A = get_ifi(i,j-1);
+		FVInterface B = get_ifi(i,j);
+		FVCell C = get_cell(i-1,j);
+		FVCell D = get_cell(i-1,j-1);
+		vtx.cloud_pos = [A.pos, B.pos, C.pos[gtl], D.pos[gtl]];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    } // j loop
+	    // West boundary
+	    i = imin;
+	    for (j = jmin+1; j <= jmax; ++j) {
+		FVVertex vtx = get_vtx(i,j);
+		// These are the corners of the secondary cell.
+		FVCell A = get_cell(i,j-1);
+		FVCell B = get_cell(i,j);
+		FVInterface C = get_ifi(i,j);
+		FVInterface D = get_ifi(i,j-1);
+		vtx.cloud_pos = [A.pos[gtl], B.pos[gtl], C.pos, D.pos];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    } // j loop
+	    // North boundary
+	    j = jmax+1;
+	    for (i = imin+1; i <= imax; ++i) {
+		FVVertex vtx = get_vtx(i,j);
+		FVCell A = get_cell(i,j-1);
+		FVInterface B = get_ifj(i,j);
+		FVInterface C = get_ifj(i-1,j);
+		FVCell D = get_cell(i-1,j-1);
+		vtx.cloud_pos = [A.pos[gtl], B.pos, C.pos, D.pos[gtl]];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    } // i loop
+	    // South boundary
+	    j = jmin;
+	    for (i = imin+1; i <= imax; ++i) {
+		FVVertex vtx = get_vtx(i,j);
+		FVInterface A = get_ifj(i,j);
+		FVCell B = get_cell(i,j);
+		FVCell C = get_cell(i-1,j);
+		FVInterface D = get_ifj(i-1,j);
+		vtx.cloud_pos = [A.pos, B.pos[gtl], C.pos[gtl], D.pos];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    } // i loop
+	    // For the corners, we are going to use the same divergence-theorem-based
+	    // gradient calculator and let one edge collapse to a point, thus giving
+	    // it a triangle to compute over.  This should be fine. 
+	    // North-east corner
+	    {
+		i = imax+1; j = jmax+1;
+		FVVertex vtx = get_vtx(i,j);
+		FVInterface A = get_ifi(i,j-1);
+		FVInterface B = get_ifj(i-1,j);
+		FVInterface C = get_ifj(i-1,j);
+		FVCell D = get_cell(i-1,j-1);
+		vtx.cloud_pos = [A.pos, B.pos, C.pos, D.pos[gtl]];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    }
+	    // South-east corner
+	    {
+		i = imax+1; j = jmin;
+		FVVertex vtx = get_vtx(i,j);
+		FVInterface A = get_ifi(i,j);
+		FVInterface B = get_ifi(i,j);
+		FVCell C = get_cell(i-1,j);
+		FVInterface D = get_ifj(i-1,j);
+		vtx.cloud_pos = [A.pos, B.pos, C.pos[gtl], D.pos];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    }
+	    // South-west corner
+	    {
+		i = imin; j = jmin;
+		FVVertex vtx = get_vtx(i,j);
+		FVInterface A = get_ifj(i,j);
+		FVCell B = get_cell(i,j);
+		FVInterface C = get_ifi(i,j);
+		FVInterface D = A;
+		vtx.cloud_pos = [A.pos, B.pos[gtl], C.pos, D.pos];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    }
+	    // North-west corner
+	    {
+		i = imin; j = jmax+1;
+		FVVertex vtx = get_vtx(i,j);
+		FVCell A = get_cell(i,j-1);
+		FVInterface B = get_ifj(i,j);
+		FVInterface C = B;
+		FVInterface D = get_ifi(i,j-1);
+		vtx.cloud_pos = [A.pos[gtl], B.pos, C.pos, D.pos];
+		vtx.cloud_fs = [A.fs, B.fs, C.fs, D.fs];
+	    }
+	} else { // Flow quantity derivatives for 3D.
+	    // Internal secondary cell geometry information
+	    for ( i = imin; i <= imax-1; ++i ) {
+		for ( j = jmin; j <= jmax-1; ++j ) {
+		    for ( k = kmin; k <= kmax-1; ++k ) {
+			FVVertex vtx = get_vtx(i+1,j+1,k+1);
+			FVCell c0 = get_cell(i,j,k);
+			FVCell c1 = get_cell(i+1,j,k);
+			FVCell c2 = get_cell(i+1,j+1,k);
+			FVCell c3 = get_cell(i,j+1,k);
+			FVCell c4 = get_cell(i,j,k+1);
+			FVCell c5 = get_cell(i+1,j,k+1);
+			FVCell c6 = get_cell(i+1,j+1,k+1);
+			FVCell c7 = get_cell(i,j+1,k+1);
+			vtx.cloud_pos = [c0.pos[gtl], c1.pos[gtl], c2.pos[gtl], c3.pos[gtl],
+					 c4.pos[gtl], c5.pos[gtl], c6.pos[gtl], c7.pos[gtl]];
+			vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		    }
+		}
+	    }
+	    // East boundary secondary cell geometry information
+	    i = imax;
+	    for ( j = jmin; j <= jmax-1; ++j ) {
+		for ( k = kmin; k <= kmax-1; ++k ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVCell c0 = get_cell(i,j,k);
+		    FVInterface c1 = get_ifi(i+1,j,k);
+		    FVInterface c2 = get_ifi(i+1,j+1,k);
+		    FVCell c3 = get_cell(i,j+1,k);
+		    FVCell c4 = get_cell(i,j,k+1);
+		    FVInterface c5 = get_ifi(i+1,j,k+1);
+		    FVInterface c6 = get_ifi(i+1,j+1,k+1);
+		    FVCell c7 = get_cell(i,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos[gtl], c1.pos, c2.pos, c3.pos[gtl],
+				     c4.pos[gtl], c5.pos, c6.pos, c7.pos[gtl]];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // West boundary secondary cell geometry information
+	    i = imin - 1;
+	    for ( j = jmin; j <= jmax-1; ++j ) {
+		for ( k = kmin; k <= kmax-1; ++k ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVInterface c0 = get_ifi(i+1,j,k);
+		    FVCell c1 = get_cell(i+1,j,k);
+		    FVCell c2 = get_cell(i+1,j+1,k);
+		    FVInterface c3 = get_ifi(i+1,j+1,k);
+		    FVInterface c4 = get_ifi(i+1,j,k+1);
+		    FVCell c5 = get_cell(i+1,j,k+1);
+		    FVCell c6 = get_cell(i+1,j+1,k+1);
+		    FVInterface c7 = get_ifi(i+1,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos, c1.pos[gtl], c2.pos[gtl], c3.pos,
+				     c4.pos, c5.pos[gtl], c6.pos[gtl], c7.pos];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // North boundary secondary cell geometry information
+	    j = jmax;
+	    for ( i = imin; i <= imax-1; ++i ) {
+		for ( k = kmin; k <= kmax-1; ++k ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVCell c0 = get_cell(i,j,k);
+		    FVCell c1 = get_cell(i+1,j,k);
+		    FVInterface c2 = get_ifj(i+1,j+1,k);
+		    FVInterface c3 = get_ifj(i,j+1,k);
+		    FVCell c4 = get_cell(i,j,k+1);
+		    FVCell c5 = get_cell(i+1,j,k+1);
+		    FVInterface c6 = get_ifj(i+1,j+1,k+1);
+		    FVInterface c7 = get_ifj(i,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos[gtl], c1.pos[gtl], c2.pos, c3.pos,
+				     c4.pos[gtl], c5.pos[gtl], c6.pos, c7.pos];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // South boundary secondary cell geometry information
+	    j = jmin - 1;
+	    for ( i = imin; i <= imax-1; ++i ) {
+		for ( k = kmin; k <= kmax-1; ++k ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVInterface c0 = get_ifj(i,j+1,k);
+		    FVInterface c1 = get_ifj(i+1,j+1,k);
+		    FVCell c2 = get_cell(i+1,j+1,k);
+		    FVCell c3 = get_cell(i,j+1,k);
+		    FVInterface c4 = get_ifj(i,j+1,k+1);
+		    FVInterface c5 = get_ifj(i+1,j+1,k+1);
+		    FVCell c6 = get_cell(i+1,j+1,k+1);
+		    FVCell c7 = get_cell(i,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos, c1.pos, c2.pos[gtl], c3.pos[gtl],
+				     c4.pos, c5.pos, c6.pos[gtl], c7.pos[gtl]];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // Top boundary secondary cell geometry information
+	    k = kmax;
+	    for ( i = imin; i <= imax-1; ++i ) {
+		for ( j = jmin; j <= jmax-1; ++j ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVCell c0 = get_cell(i,j,k);
+		    FVCell c1 = get_cell(i+1,j,k);
+		    FVCell c2 = get_cell(i+1,j+1,k);
+		    FVCell c3 = get_cell(i,j+1,k);
+		    FVInterface c4 = get_ifk(i,j,k+1);
+		    FVInterface c5 = get_ifk(i+1,j,k+1);
+		    FVInterface c6 = get_ifk(i+1,j+1,k+1);
+		    FVInterface c7 = get_ifk(i,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos[gtl], c1.pos[gtl], c2.pos[gtl], c3.pos[gtl],
+				     c4.pos, c5.pos, c6.pos, c7.pos];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // Bottom boundary secondary cell geometry information
+	    k = kmin - 1;
+	    for ( i = imin; i <= imax-1; ++i ) {
+		for ( j = jmin; j <= jmax-1; ++j ) {
+		    FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		    FVInterface c0 = get_ifk(i,j,k+1);
+		    FVInterface c1 = get_ifk(i+1,j,k+1);
+		    FVInterface c2 = get_ifk(i+1,j+1,k+1);
+		    FVInterface c3 = get_ifk(i,j+1,k+1);
+		    FVCell c4 = get_cell(i,j,k+1);
+		    FVCell c5 = get_cell(i+1,j,k+1);
+		    FVCell c6 = get_cell(i+1,j+1,k+1);
+		    FVCell c7 = get_cell(i,j+1,k+1);
+		    vtx.cloud_pos = [c0.pos, c1.pos, c2.pos, c3.pos,
+				     c4.pos[gtl], c5.pos[gtl], c6.pos[gtl], c7.pos[gtl]];
+		    vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs, c6.fs, c7.fs];
+		}
+	    }
+	    // Now, do the 4 edges around the bottom face.
+	    // Bottom-South edge [0]-->[1]
+	    j = jmin; k = kmin;    	
+	    for ( i = imin+1; i <= imax; ++i ) {
+		FVVertex vtx = get_vtx(i,j,k);
+		FVCell c0 = get_cell(i-1,j,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifj(i-1,j,k);
+		FVInterface c3 = get_ifk(i-1,j,k);
+		FVInterface c4 = get_ifj(i,j,k);
+		FVInterface c5 = get_ifk(i,j,k);
+		vtx.cloud_pos = [c0.pos[gtl], c1.pos[gtl], c2.pos, c3.pos, c4.pos, c5.pos];
+		vtx.cloud_fs = [c0.fs, c1.fs, c2.fs, c3.fs, c4.fs, c5.fs];
+	    }
+	    // Bottom-North edge [3]-->[2]
+	    j = jmax; k = kmin;
+	    for ( i = imin+1; i <= imax; ++i ) {
+		FVVertex vtx = get_vtx(i,j+1,k);
+		FVCell c0 = get_cell(i-1,j,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifj(i-1,j+1,k);
+		FVInterface c3 = get_ifk(i-1,j,k);
+		FVInterface c4 = get_ifj(i,j+1,k);
+		FVInterface c5 = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Bottom-West edge [0]-->[3]
+	    i = imin; k = kmin;
+	    for ( j = jmin+1; j <= jmax; ++j ) {
+		FVVertex vtx = get_vtx(i,j,k);
+		FVCell c0 = get_cell(i,j-1,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i,j-1,k);
+		FVInterface c3 = get_ifk(i,j-1,k);
+		FVInterface c4 = get_ifi(i,j,k);
+		FVInterface c5 = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Bottom-East edge [1]-->[2]
+	    i = imax; k = kmin;
+	    for ( j = jmin+1; j <= jmax; ++j ) {
+		FVVertex vtx = get_vtx(i+1,j,k);
+		FVCell c0 = get_cell(i,j-1,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i+1,j-1,k);
+		FVInterface c3 = get_ifk(i,j-1,k);
+		FVInterface c4 = get_ifi(i+1,j,k);
+		FVInterface c5 = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // 4 edges around the top face.
+	    // Top-South edge [4]-->[5]
+	    j = jmin; k = kmax;
+	    for ( i = imin+1; i <= imax; ++i ) {
+		FVVertex vtx = get_vtx(i,j,k+1);
+		FVCell c0 = get_cell(i-1,j,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifj(i-1,j,k);
+		FVInterface c3 = get_ifk(i-1,j,k+1);
+		FVInterface c4 = get_ifj(i,j,k);
+		FVInterface c5 = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Top-North edge [7]-->[6]
+	    j = jmax; k = kmax;
+	    for ( i = imin+1; i <= imax; ++i ) {
+		FVVertex vtx = get_vtx(i,j+1,k+1);
+		FVCell c0 = get_cell(i-1,j,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifj(i-1,j+1,k);
+		FVInterface c3 = get_ifk(i-1,j,k+1);
+		FVInterface c4 = get_ifj(i,j+1,k);
+		FVInterface c5 = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Top-West edge [4]-->[7]
+	    i = imin; k = kmax;
+	    for ( j = jmin+1; j <= jmax; ++j ) {
+		FVVertex vtx = get_vtx(i,j,k+1);
+		FVCell c0 = get_cell(i,j-1,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i,j-1,k);
+		FVInterface c3 = get_ifk(i,j-1,k+1);
+		FVInterface c4 = get_ifi(i,j,k);
+		FVInterface c5 = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Top-East edge [5]-->[6]
+	    i = imax; k = kmax;
+	    for ( j = jmin+1; j <= jmax; ++j ) {
+		FVVertex vtx = get_vtx(i+1,j,k+1);
+		FVCell c0 = get_cell(i,j-1,k);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i+1,j-1,k);
+		FVInterface c3 = get_ifk(i,j-1,k+1);
+		FVInterface c4 = get_ifi(i+1,j,k);
+		FVInterface c5 = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // 4 edges running from bottom to top.
+	    // South-West edge [0]-->[4]
+	    i = imin; j = jmin;
+	    for ( k = kmin+1; k <= kmax; ++k ) {
+		FVVertex vtx = get_vtx(i,j,k);
+		FVCell c0 = get_cell(i,j,k-1);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i,j,k-1);
+		FVInterface c3 = get_ifj(i,j,k-1);
+		FVInterface c4 = get_ifi(i,j,k);
+		FVInterface c5 = get_ifj(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // South-East edge [1]-->[5]
+	    i = imax; j = jmin;
+	    for ( k = kmin+1; k <= kmax; ++k ) {
+		FVVertex vtx = get_vtx(i+1,j,k);
+		FVCell c0 = get_cell(i,j,k-1);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i+1,j,k-1);
+		FVInterface c3 = get_ifj(i,j,k-1);
+		FVInterface c4 = get_ifi(i+1,j,k);
+		FVInterface c5 = get_ifj(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-East edge [2]-->[6]
+	    i = imax; j = jmax;
+	    for ( k = kmin+1; k <= kmax; ++k ) {
+		FVVertex vtx = get_vtx(i+1,j+1,k);
+		FVCell c0 = get_cell(i,j,k-1);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i+1,j,k-1);
+		FVInterface c3 = get_ifj(i,j+1,k-1);
+		FVInterface c4 = get_ifi(i+1,j,k);
+		FVInterface c5 = get_ifj(i,j+1,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-West edge [3]-->[7]
+	    i = imin; j = jmax;
+	    for ( k = kmin+1; k <= kmax; ++k ) {
+		FVVertex vtx = get_vtx(i,j+1,k);
+		FVCell c0 = get_cell(i,j,k-1);
+		FVCell c1 = get_cell(i,j,k);
+		FVInterface c2 = get_ifi(i,j,k-1);
+		FVInterface c3 = get_ifj(i,j+1,k-1);
+		FVInterface c4 = get_ifi(i,j,k);
+		FVInterface c5 = get_ifj(i,j+1,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // Finally, the 8 corners.
+	    // South-West-Bottom corner [0]
+	    i = imin; j = jmin; k = kmin;
+	    {
+		FVCell c0 = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i,j,k);
+		FVInterface a = get_ifi(i,j,k);
+		FVInterface b = get_ifj(i,j,k);
+		FVInterface d = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // South-East-Bottom corner [1]
+	    i = imax; j = jmin; k = kmin;
+	    {
+		FVCell c0 = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i+1,j,k);
+		FVInterface a = get_ifi(i+1,j,k);
+		FVInterface b = get_ifj(i,j,k);
+		FVInterface d = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-East-Bottom corner [2]
+	    i = imax; j = jmax; k = kmin;
+	    {
+		FVCell c0 = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i+1,j+1,k);
+		FVInterface a = get_ifi(i+1,j,k);
+		FVInterface b = get_ifj(i,j+1,k);
+		FVInterface d = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-West-Bottom corner [3]
+	    i = imin; j = jmax; k = kmin;
+	    {
+		FVCell c = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i,j+1,k);
+		FVInterface a = get_ifi(i,j,k);
+		FVInterface b = get_ifj(i,j+1,k);
+		FVInterface d = get_ifk(i,j,k);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // South-West-Top corner [4]
+	    i = imin; j = jmin; k = kmax;
+	    {
+		FVCell c = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i,j,k+1);
+		FVInterface a = get_ifi(i,j,k);
+		FVInterface b = get_ifj(i,j,k);
+		FVInterface d = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // South-East-Top corner [5]
+	    i = imax; j = jmin; k = kmax;
+	    {
+		FVCell c = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i+1,j,k+1);
+		FVInterface a = get_ifi(i+1,j,k);
+		FVInterface b = get_ifj(i,j,k);
+		FVInterface d = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-East-Top corner [6]
+	    i = imax; j = jmax; k = kmax;
+	    {
+		FVCell c = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i+1,j+1,k+1);
+		FVInterface a = get_ifi(i+1,j,k);
+		FVInterface b = get_ifj(i,j+1,k);
+		FVInterface d = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	    // North-West-Top corner [7]
+	    i = imin; j = jmax; k = kmax;
+	    {
+		FVCell c = get_cell(i,j,k);
+		FVVertex vtx = get_vtx(i,j+1,k+1);
+		FVInterface a = get_ifi(i,j,k);
+		FVInterface b = get_ifj(i,j+1,k);
+		FVInterface d = get_ifk(i,j,k+1);
+		vtx.cloud_pos = [];
+		vtx.cloud_fs = [];
+	    }
+	} // end if (myConfig.dimensions
+    } // end assign_flow_locations_for_derivative_calc()
 
     override void read_grid(string filename, size_t gtl=0)
     // Read the grid vertices from a gzip file.
@@ -1284,127 +1775,30 @@
     @nogc
     override void flow_property_derivatives(int gtl)
     {
-	size_t i, j, k;
 	if (myConfig.dimensions == 2) {
-	    // First, do all of the internal secondary cells.
-	    // i.e. Those not on a boundary.
-	    for ( i = imin+1; i <= imax; ++i ) {
-		for ( j = jmin+1; j <= jmax; ++j ) {
-		    // Secondary-cell centre is a primary-cell vertex.
+	    size_t k = 0;
+	    for ( size_t i = imin; i <= imax+1; ++i ) {
+		for ( size_t j = jmin; j <= jmax+1; ++j ) {
 		    FVVertex vtx = get_vtx(i,j);
-		    // These are the corners of the secondary cell.
-		    FVCell A = get_cell(i,j-1);
-		    FVCell B = get_cell(i,j);
-		    FVCell C = get_cell(i-1,j);
-		    FVCell D = get_cell(i-1,j-1);
-		    // Delegate the work of computing the actual gradients.
-		    gradients_xy(vtx, A.pos[gtl], B.pos[gtl], C.pos[gtl], D.pos[gtl],
-				 A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-		} // j loop
-	    } // i loop
-	    // Half-cells along the edges of the block.
-	    // [TODO] Check that the secondary-cell geometry functions actually
-	    // compute the correct xy-plane area for these "half-cells" or, maybe,
-	    // the divergence function needs to compute its own copy of area.
-	    // It will be a bit of extra work but then we can pass it all sorts of 
-	    // shapes and it will look after itself, even if the grid moves.
-	    // East boundary
-	    i = imax+1;
-	    for (j = jmin+1; j <= jmax; ++j) {
-		FVVertex vtx = get_vtx(i,j);
-		FVInterface A = get_ifi(i,j-1);
-		FVInterface B = get_ifi(i,j);
-		FVCell C = get_cell(i-1,j);
-		FVCell D = get_cell(i-1,j-1);
-		gradients_xy(vtx, A.pos, B.pos, C.pos[gtl], D.pos[gtl],
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    } // j loop
-	    // West boundary
-	    i = imin;
-	    for (j = jmin+1; j <= jmax; ++j) {
-		FVVertex vtx = get_vtx(i,j);
-		// These are the corners of the secondary cell.
-		FVCell A = get_cell(i,j-1);
-		FVCell B = get_cell(i,j);
-		FVInterface C = get_ifi(i,j);
-		FVInterface D = get_ifi(i,j-1);
-		// Delegate the work of computing the actual gradients.
-		gradients_xy(vtx, A.pos[gtl], B.pos[gtl], C.pos, D.pos,
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    } // j loop
-	    // North boundary
-	    j = jmax+1;
-	    for (i = imin+1; i <= imax; ++i) {
-		FVVertex vtx = get_vtx(i,j);
-		FVCell A = get_cell(i,j-1);
-		FVInterface B = get_ifj(i,j);
-		FVInterface C = get_ifj(i-1,j);
-		FVCell D = get_cell(i-1,j-1);
-		gradients_xy(vtx, A.pos[gtl], B.pos, C.pos, D.pos[gtl],
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    } // i loop
-	    // South boundary
-	    j = jmin;
-	    for (i = imin+1; i <= imax; ++i) {
-		FVVertex vtx = get_vtx(i,j);
-		FVInterface A = get_ifj(i,j);
-		FVCell B = get_cell(i,j);
-		FVCell C = get_cell(i-1,j);
-		FVInterface D = get_ifj(i-1,j);
-		gradients_xy(vtx, A.pos, B.pos[gtl], C.pos[gtl], D.pos,
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    } // i loop
-	    // For the corners, we are going to use the same divergence-theorem-based
-	    // gradient calculator and let one edge collapse to a point, thus giving
-	    // it a triangle to compute over.  This should be fine. 
-	    // North-east corner
-	    {
-		i = imax+1; j = jmax+1;
-		FVVertex vtx = get_vtx(i,j);
-		FVInterface A = get_ifi(i,j-1);
-		FVInterface B = get_ifj(i-1,j);
-		FVInterface C = get_ifj(i-1,j);
-		FVCell D = get_cell(i-1,j-1);
-		gradients_xy(vtx, A.pos, B.pos, C.pos, D.pos[gtl],
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    }
-	    // South-east corner
-	    {
-		i = imax+1; j = jmin;
-		FVVertex vtx = get_vtx(i,j);
-		FVInterface A = get_ifi(i,j);
-		FVInterface B = get_ifi(i,j);
-		FVCell C = get_cell(i-1,j);
-		FVInterface D = get_ifj(i-1,j);
-		gradients_xy(vtx, A.pos, B.pos, C.pos[gtl], D.pos,
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    }
-	    // South-west corner
-	    {
-		i = imin; j = jmin;
-		FVVertex vtx = get_vtx(i,j);
-		FVInterface A = get_ifj(i,j);
-		FVCell B = get_cell(i,j);
-		FVInterface C = get_ifi(i,j);
-		FVInterface D = A;
-		gradients_xy(vtx, A.pos, B.pos[gtl], C.pos, D.pos,
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
-	    }
-	    // North-west corner
-	    {
-		i = imin; j = jmax+1;
-		FVVertex vtx = get_vtx(i,j);
-		FVCell A = get_cell(i,j-1);
-		FVInterface B = get_ifj(i,j);
-		FVInterface C = B;
-		FVInterface D = get_ifi(i,j-1);
-		gradients_xy(vtx, A.pos[gtl], B.pos, C.pos, D.pos,
-			     A.fs, B.fs, C.fs, D.fs, myConfig.diffusion);
+		    // For the moment, don't change the signature of the gradients fn.
+		    gradients_xy(vtx, vtx.cloud_pos[0], vtx.cloud_pos[1],
+				 vtx.cloud_pos[2], vtx.cloud_pos[3],
+				 vtx.cloud_fs[0], vtx.cloud_fs[1],
+				 vtx.cloud_fs[2], vtx.cloud_fs[3], myConfig.diffusion);
+		}
 	    }
 	} else {
-	    assert(false, "[TODO] flow_property_derivatives() in 3D not implemented yet.");
-	}
-    }
+	    // Flow quantity derivatives for 3D.
+	    for ( size_t i = imin; i <= imax+1; ++i ) {
+		for ( size_t j = jmin; j <= jmax+1; ++j ) {
+		    for ( size_t k = kmin; k <= kmax+1; ++k ) {
+			FVVertex vtx = get_vtx(i,j,k);
+			gradients_xyz(vtx, myConfig.diffusion);
+		    }
+		}
+	    }
+	} // end if (myConfig.dimensions
+    } // end flow_property_derivatives()
 
     override void applyPreReconAction(double t, int gtl, int ftl)
     {
--- a/dlang/eilmer4/simcore.d	Sat Jun 06 10:15:41 2015 +1000
+++ b/dlang/eilmer4/simcore.d	Sat Jun 06 21:15:20 2015 +1000
@@ -75,6 +75,7 @@
     foreach (myblk; parallel(gasBlocks,1)) {
 	myblk.compute_primary_cell_geometric_data(0);
 	myblk.compute_distance_to_nearest_wall_for_all_cells(0);
+	myblk.assign_flow_locations_for_derivative_calc(0);
 	myblk.identify_reaction_zones(0);
 	myblk.identify_turbulent_zones(0);
 	myblk.set_grid_velocities(sim_time);
--- a/dlang/eilmer4/viscousflux.d	Sat Jun 06 10:15:41 2015 +1000
+++ b/dlang/eilmer4/viscousflux.d	Sat Jun 06 21:15:20 2015 +1000
@@ -297,7 +297,7 @@
 //
 // Since we are embedding this function in a nominally 3D code,
 // the z-coordinate derivatives are set to zero.
-// [TODO] Check that this contour-integral approach works when the 
+// This contour-integral approach works when the 
 // quadrilateral degenerates to a triangle in the plane. 
 // We can then double-up on one of the points and still get an estimate
 // for the flow gradients.
@@ -376,12 +376,172 @@
     vtx.grad_omega.refz = 0.0;
 } // end gradients_xy()
 
+@nogc
+void computeInverse(ref double[6][3] c, double very_small_value=1.0e-16)
+// Perform Gauss-Jordan elimination on an augmented matrix.
+// c = [A|b] such that the mutated matrix becomes [I|x]
+// where x is the solution vector(s) to A.x = b
+{
+    foreach(j; 0 .. 3) {
+	// Select pivot.
+	size_t p = j;
+	foreach(i; j+1 .. 3) {
+	    if ( abs(c[i][j]) > abs(c[p][j]) ) p = i;
+	}
+	assert(abs(c[p][j]) > very_small_value, "matrix is essentially singular");
+	if ( p != j ) { // Swap rows
+	    foreach(col; 0 .. 6) {
+		double tmp = c[p][col]; c[p][col] = c[j][col]; c[j][col] = tmp;
+	    }
+	}
+	// Scale row j to get unity on the diagonal.
+	double cjj = c[j][j];
+	foreach(col; 0 .. 6) c[j][col] /= cjj;
+	// Do the elimination to get zeros in all off diagonal values in column j.
+	foreach(i; 0 .. 3) {
+	    if ( i == j ) continue;
+	    double cij = c[i][j];
+	    foreach(col; 0 .. 6) c[i][col] -= cij * c[j][col]; 
+	}
+    } // end foreach j
+} // end gaussJordanElimination()
+
+@nogc
+void solveGradients(ref double[6][3] c, ref double[3] rhs, ref double[3] qgrad)
+// Multiply right-hand-side by the inverse part of the augmented matrix.
+{
+    foreach(i; 0 .. 3) {
+	qgrad[i] = 0.0;
+	foreach(j; 0 .. 3) {
+	    qgrad[i] += c[i][3+j] * rhs[j];
+	}
+    }
+} // end solveGradients()
 
 @nogc
-void gradients_xyz(ref FVVertex vtx,
-		   const ref Vector3[] pos, const ref FlowState[] fs,
-		   bool diffusion)
+void gradients_xyz(ref FVVertex vtx, bool diffusion)
+// Fit a linear model to the cloud of flow-quantity points
+// in order to extract approximations to the flow-field gradients.
 {
-    // pos, fs represent the cloud of surrounding points and 
-    // corresponding flow states.
+    double[6][3] xTx; // normal matrix, augmented to give 6 entries per row
+    double[3] rhs, gradients;
+    // Assemble and invert the normal matrix.
+    // We'll reuse the resulting inverse for each flow-field quantity.
+    size_t N = vtx.cloud_pos.length;
+    double xx = 0.0;
+    double xy = 0.0;
+    double xz = 0.0;
+    double yy = 0.0;
+    double yz = 0.0;
+    double zz = 0.0;
+    foreach (i; 1 .. N) {
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	xx += dx*dx; xy += dx*dy; xz += dx*dz;
+	yy += dy*dy; yz += dy*dz; zz += dz*dz;
+    }
+    xTx[0][0] = xx; xTx[0][1] = xy; xTx[0][2] = xz;
+    xTx[1][0] = xy; xTx[1][1] = yy; xTx[1][2] = yz;
+    xTx[2][0] = xz; xTx[2][1] = yz; xTx[2][2] = zz;
+    xTx[0][3] = 1.0; xTx[0][4] = 0.0; xTx[0][5] = 0.0;
+    xTx[1][3] = 0.0; xTx[1][4] = 1.0; xTx[1][5] = 0.0;
+    xTx[2][3] = 0.0; xTx[2][4] = 0.0; xTx[2][5] = 1.0;
+    computeInverse(xTx);
+    // x-velocity
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double dvx = vtx.cloud_fs[i].vel.x - vtx.cloud_fs[i].vel.x;
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*dvx; rhs[1] += dy*dvx; rhs[2] += dz*dvx;
+    }
+    solveGradients(xTx, rhs, gradients);
+    foreach (j; 0 .. 3) { vtx.grad_vel[0][j] = gradients[j]; }
+    // y-velocity
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double dvy = vtx.cloud_fs[i].vel.y - vtx.cloud_fs[i].vel.y;
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*dvy; rhs[1] += dy*dvy; rhs[2] += dz*dvy;
+    }
+    solveGradients(xTx, rhs, gradients);
+    foreach (j; 0 .. 3) { vtx.grad_vel[1][j] = gradients[j]; }
+    // z-velocity
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double dvz = vtx.cloud_fs[i].vel.z - vtx.cloud_fs[i].vel.z;
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*dvz; rhs[1] += dy*dvz; rhs[2] += dz*dvz;
+    }
+    solveGradients(xTx, rhs, gradients);
+    foreach (j; 0 .. 3) { vtx.grad_vel[2][j] = gradients[j]; }
+    // T[0]
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double dT = vtx.cloud_fs[i].gas.T[0] - vtx.cloud_fs[i].gas.T[0];
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*dT; rhs[1] += dy*dT; rhs[2] += dz*dT;
+    }
+    solveGradients(xTx, rhs, gradients);
+    vtx.grad_T.refx = gradients[0];
+    vtx.grad_T.refy = gradients[1];
+    vtx.grad_T.refz = gradients[2];
+    // massf
+    size_t nsp = vtx.cloud_fs[0].gas.massf.length;
+    if (diffusion) {
+	foreach(isp; 0 .. nsp) {
+	    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+	    foreach (i; 1 .. N) {
+		double df = vtx.cloud_fs[i].gas.massf[isp] - vtx.cloud_fs[i].gas.massf[isp];
+		double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+		double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+		double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+		rhs[0] += dx*df; rhs[1] += dy*df; rhs[2] += dz*df;
+	    }
+	    solveGradients(xTx, rhs, gradients);
+	    vtx.grad_f[isp].refx = 0.0;
+	    vtx.grad_f[isp].refy = 0.0;
+	    vtx.grad_f[isp].refz = 0.0;
+	} // foreach isp
+    } else {
+	foreach(isp; 0 .. nsp) {
+	    vtx.grad_f[isp].refx = 0.0;
+	    vtx.grad_f[isp].refy = 0.0;
+	    vtx.grad_f[isp].refz = 0.0;
+	} // foreach isp
+    }
+    // tke
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double dtke = vtx.cloud_fs[i].tke - vtx.cloud_fs[i].tke;
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*dtke; rhs[1] += dy*dtke; rhs[2] += dz*dtke;
+    }
+    solveGradients(xTx, rhs, gradients);
+    vtx.grad_tke.refx = gradients[0];
+    vtx.grad_tke.refy = gradients[1];
+    vtx.grad_tke.refz = gradients[2];
+    // omega
+    foreach (j; 0 .. 3) { rhs[j] = 0.0; }
+    foreach (i; 1 .. N) {
+	double domega = vtx.cloud_fs[i].omega - vtx.cloud_fs[i].omega;
+	double dx = vtx.cloud_pos[i].x - vtx.cloud_pos[0].x;
+	double dy = vtx.cloud_pos[i].y - vtx.cloud_pos[0].y;
+	double dz = vtx.cloud_pos[i].z - vtx.cloud_pos[0].z;
+	rhs[0] += dx*domega; rhs[1] += dy*domega; rhs[2] += dz*domega;
+    }
+    solveGradients(xTx, rhs, gradients);
+    vtx.grad_omega.refx = gradients[0];
+    vtx.grad_omega.refy = gradients[1];
+    vtx.grad_omega.refz = gradients[2];
 } // end gradients_xyz()