#include "wave.h"

void diffop(Array2_dbl &du,Array2_dbl &u,int NNx,int NNy,
	    int nreq, MPI_Request *reqarr, int *coords,int *dims)
{
  int i,j,Nx,Ny;
  MPI_Status mpistat[8];
  Nx=u.storleken1(); Ny=u.storleken2();

  // Start up all the communication
  MPI_Startall(nreq,reqarr);

  // Differentiate inner points while waiting for the communication 
  for (i=2; i<Nx-2; i++)
    for (j=2; j<Ny-2; j++)
      du(i,j)=(u(i+1,j)-u(i-1,j))/2.0*NNx+(u(i,j+1)-u(i,j-1))/2.0*NNy;

  // Wait for the communication to complete
  MPI_Waitall(nreq,reqarr,mpistat);

  // Differentiate inner boundaries
  if (coords[0]>0)
    {
      for (j=1; j<Ny-1; j++)
	du(1,j)=(u(2,j)-u(0,j))/2.0*NNx+(u(1,j+1)-u(1,j-1))/2.0*NNy;
    }
  if (coords[0]<(dims[0]-1))
    {
      for (j=1; j<Ny-1; j++)
	du(Nx-2,j)=(u(Nx-1,j)-u(Nx-3,j))/2.0*NNx+
	  (u(Nx-2,j+1)-u(Nx-2,j-1))/2.0*NNy;
    }
  if (coords[1]>0)
    {
      for (i=1; i<Nx-1; i++)
	du(i,1)=(u(i+1,1)-u(i-1,1))/2.0*NNx+(u(i,2)-u(i,0))/2.0*NNy;
    }
  if (coords[1]<(dims[1]-1))
    {
      for (i=1; i<Nx-1; i++)
	du(i,Ny-2)=(u(i+1,Ny-2)-u(i-1,Ny-2))/2.0*NNx+
	  (u(i,Ny-1)-u(i,Ny-3))/2.0*NNy;
    }

}