#include "petscsnes.h"
#include "petscda.h"

typedef struct {
   InactiveDouble lidvelocity,prandtl,grashof;  /* physical parameters */
   PetscTruth draw_contours;                /* flag - 1 indicates drawing contours */
} AppCtx;

typedef struct {
  Scalar u,v,omega,temp;
} Field;

 int localfunction2d(Field **x,Field **f,int xs, int xm, int ys, int ym,
		                int mx,int my, void *ptr) {

  AppCtx  *user = (AppCtx*)ptr;
  int     ierr,i,j;
  int     xints,xinte,yints,yinte;
  double  two = 2.0,one = 1.0,p5 = 0.5,hx,hy,dhx,dhy,hxdhy,hydhx;
  double  grashof,prandtl,lid;
  Scalar  u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym;

  grashof = user->grashof;  
  prandtl = user->prandtl;
  lid     = user->lidvelocity;

  /* 
     Define mesh intervals ratios for uniform grid.
     [Note: FD formulae below are normalized by multiplying through by
     local volume element to obtain coefficients O(1) in two dimensions.]
  */
  dhx = (double)(mx-1);     dhy = (double)(my-1);
  hx = one/dhx;             hy = one/dhy;
  hxdhy = hx*dhy;           hydhx = hy*dhx;

  xints = xs; xinte = xs+xm; yints = ys; yinte = ys+ym;

  /* Test whether we are on the bottom edge of the global array */
  if (yints == 0) {
    j = 0;
    yints = yints + 1;
    /* bottom edge */
    for (i=xs; i<xs+xm; i++) {
        f[j][i].u     = x[j][i].u;
        f[j][i].v     = x[j][i].v;
        f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy; 
	f[j][i].temp  = x[j][i].temp-x[j+1][i].temp;
    }
  }

  /* Test whether we are on the top edge of the global array */
  if (yinte == my) {
    j = my - 1;
    yinte = yinte - 1;
    /* top edge */
    for (i=xs; i<xs+xm; i++) {
        f[j][i].u     = x[j][i].u - lid;
        f[j][i].v     = x[j][i].v;
        f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy; 
	f[j][i].temp  = x[j][i].temp-x[j-1][i].temp;
    }
  }

  /* Test whether we are on the left edge of the global array */
  if (xints == 0) {
    i = 0;
    xints = xints + 1;
    /* left edge */
    for (j=ys; j<ys+ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx; 
      f[j][i].temp  = x[j][i].temp;
    }
  }

  /* Test whether we are on the right edge of the global array */
  if (xinte == mx) {
    i = mx - 1;
    xinte = xinte - 1;
    /* right edge */ 
    for (j=ys; j<ys+ym; j++) {
      f[j][i].u     = x[j][i].u;
      f[j][i].v     = x[j][i].v;
      f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx; 
      f[j][i].temp  = x[j][i].temp - (double)(grashof>0);
    }
  }

  /* Compute over the interior points */
  for (j=yints; j<yinte; j++) {
    for (i=xints; i<xinte; i++) {

	/*
	  convective coefficients for upwinding
        */
	vx = x[j][i].u; avx = PetscAbsScalar(vx);
        vxp = p5*(vx+avx); vxm = p5*(vx-avx);
	vy = x[j][i].v; avy = PetscAbsScalar(vy);
        vyp = p5*(vy+avy); vym = p5*(vy-avy);

	/* U velocity */
        u          = x[j][i].u;
        uxx        = (two*u - x[j][i-1].u - x[j][i+1].u)*hydhx;
        uyy        = (two*u - x[j-1][i].u - x[j+1][i].u)*hxdhy;
        f[j][i].u  = uxx + uyy - p5*(x[j+1][i].omega-x[j-1][i].omega)*hx;

	/* V velocity */
        u          = x[j][i].v;
        uxx        = (two*u - x[j][i-1].v - x[j][i+1].v)*hydhx;
        uyy        = (two*u - x[j-1][i].v - x[j+1][i].v)*hxdhy;
        f[j][i].v  = uxx + uyy + p5*(x[j][i+1].omega-x[j][i-1].omega)*hy;

	/* Omega */
        u          = x[j][i].omega;
        uxx        = (two*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx;
        uyy        = (two*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy;
	f[j][i].omega = uxx + uyy + 
			(vxp*(u - x[j][i-1].omega) +
			  vxm*(x[j][i+1].omega - u)) * hy +
			(vyp*(u - x[j-1][i].omega) +
			  vym*(x[j+1][i].omega - u)) * hx -
			p5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy;

        /* Temperature */
        u             = x[j][i].temp;
        uxx           = (two*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx;
        uyy           = (two*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy;
	f[j][i].temp =  uxx + uyy  + prandtl * (
			(vxp*(u - x[j][i-1].temp) +
			  vxm*(x[j][i+1].temp - u)) * hy +
		        (vyp*(u - x[j-1][i].temp) +
		       	  vym*(x[j+1][i].temp - u)) * hx);
    }
  }

  /*
     Flop count (multiply-adds are counted as 2 operations)
  */
  ierr = PetscLogFlops(84*ym*xm);CHKERRQ(ierr);
  return 0; 
}