#include <fstream.h>
#include <math.h>
#include <assert.h>
#include "FlowSubdiv3D.h"

/*
 * FlowSubdiv is where the actual subdivision takes place. The mess of for loops below takes
 * care of the appropriate convolutions, putting the refined vector field into the memory provided
 * by argument VectorField3D& fine.
 */


const float tiny  = 10E-12;


void FlowSubdiv3D::subdivide(VectorField3D& coarse, VectorField3D& fine) const
{
    //Initializes a few things for the masks
    const int m11mx = m11->mx();
    const int m11my = m11->my();
    const int m11mz = m11->mz();
    const int m12mx = m12->mx();
    const int m12my = m12->my();
    const int m12mz = m12->mz();

    const int maskSize = m11mx;

    const int cnx = coarse.nx();
    const int cny = coarse.ny();
    const int cnz = coarse.nz();

    //The fine grid provided is empty right now - nothing more than some space in memory. This resizes that
    //space so that it is big enought to hold the entire refined grid.
    fine.reset(2*(2*cnx + 1) - 1 + 2*maskSize, 2*(2*cny + 1) - 1 + 2*maskSize, 2*(2*cnz + 1) - 1 + 2*maskSize);

    //Just a whole bunch of convolution
    for (int l = -cnx; l <= cnx; l++)
	for (int m = -cny; m <= cny; m++)
	    for (int n = -cnz; n <=cnz; n++)
	    {
		float coarseXVal = coarse.valX(l,m,n);
		float coarseYVal = coarse.valY(l,m,n);
		float coarseZVal = coarse.valZ(l,m,n);

		int i, j, k;

		if (fabs(coarseXVal) > tiny)
		{
		    //m11(x,y,z)
		    for (i = -m11mx; i <= m11mx; i+=m11Stride)
			for (j = -m11my; j <= m11my; j+=m11Stride)
			    for (k = -m11mz; k <= m11mz; k+=m11Stride)
				fine.valXRef(2*l+i, 2*m+j, 2*n+k) += m11->val(i,j,k) * coarseXVal;
		    //m12(y,x,z)
		    for (i = -m12my; i <= m12my; i+=m12Stride)
			for (j = -m12mx; j <= m12mx; j+=m12Stride)
			    for (k = -m12mz; k <= m12mz; k+=m12Stride)
				fine.valYRef(2*l+i, 2*m+j, 2*n+k) += m12->val(j,i,k) * coarseXVal;
		    //m12(z,x,y)
		    for (i = -m12my; i <= m12my; i+=m12Stride)
			for (j = -m12mz; j <= m12mz; j+=m12Stride)
			    for (k = -m12mx; k <= m12mx; k+=m12Stride)
				fine.valZRef(2*l+i, 2*m+j, 2*n+k) += m12->val(k,i,j) * coarseXVal;
		}

		if (fabs(coarseYVal) > tiny)
		{
		    //m12(x,y,z)
		    for (i = -m12mx; i <= m12mx; i+=m12Stride)
			for (j = -m12my; j <=m12my; j+=m12Stride)
			    for (k=-m12mz; k <=m12mz; k+=m12Stride)
				fine.valXRef(2*l+i, 2*m+j, 2*n+k) += m12->val(i,j,k) * coarseYVal;
		    //m11(y,z,x)
		    for (i = -m11my; i <= m11my; i+=m11Stride)
			for (j = -m11mz; j <= m11mz; j+=m11Stride)
			    for (k = -m11mx; k <= m11mx; k+=m11Stride)
				fine.valYRef(2*l+i, 2*m+j, 2*n+k) += m11->val(j,k,i) * coarseYVal;
		    //m12(z,x,y)
		    for (i = -m12mz; i <= m12mz; i+=m12Stride)
			for (j = -m12mx; j <= m12mx; j+=m12Stride)
			    for (k = -m12my; k <= m12my; k+=m12Stride)
				fine.valZRef(2*l+i, 2*m+j, 2*n+k) += m12->val(k,i,j) * coarseYVal;
		}

		if (fabs(coarseZVal) > tiny)
		{
		    //m12(x,y,z)
		    for (i = -m12mx; i <= m12mx; i+=m12Stride)
			for (j = -m12my; j <= m12my; j+=m12Stride)
			    for (k = -m12mz; k <= m12mz; k+=m12Stride)
				fine.valXRef(2*l+i, 2*m+j, 2*n+k) += m12->val(i,j,k) * coarseZVal;
		    //m12(y,z,x)
		    for (i = -m12my; i <= m12my; i+=m12Stride)
			for (j = -m12mz; j <= m12mz; j+=m12Stride)
			    for (k = -m12mx; k <= m12mx; k+=m12Stride)
				fine.valYRef(2*l+i, 2*m+j, 2*n+k) += m12->val(j,k,i) * coarseZVal;
		    //m11(z,x,y)
		    for (i = -m11mz; i <= m11mz; i+=m11Stride)
			for (j = -m11mx; j <= m11mx; j += m11Stride)
			    for (k = -m11my; k <=m11my; k += m11Stride)
				fine.valZRef(2*l+i, 2*m+j, 2*n+k) += m11->val(k,i,j) * coarseZVal;
		}
	    }
}