// pde_nl22.c  solve non linear second order, second degree nonlinear
// solve  uxx*uxx + 2*uxx*uxy + uyy*uyy = f(x,y)
//        f(x,y) = 4*sin(x-y)^2 
// with boundary  ub(x,y) = sin(x-y)  0 < x < 1.0 < y < 1.0 
// set up system of equations, use specialized Jacobian, simeq_newton3
// solution u(x,y) = sin(x-y)
// x at i, y at ii rderiv or nuderiv cxx, cyy  with nx=5, ny=5
// xg is X grid coordinates, yg is Y grid coordinates
// subscripts 0 and 4 are known on boundary, thus there are 9 unknowns
// u(xg[1],yg[1]), u(xg[2],yg[1]), u(xg[3],yg[1]), 
// u(xg[1],yg[2]), u(xg[2],yg[2]), u(xg[3],yg[2]),
// u(xg[1],yg[3]), u(xg[2],yg[3]), u(xg[3],yg[3]) 
// thus 1 <= i <= 3, 1 <= ii <= 3 for inside boundary, boundary at 0 and 4
//
// discrete uxx(xg[i],yg[ii) = cxx[0]*ub(xg[0],yg[ii]) + 
//   cxx[1]* u(xg[1],yg[ii]) + cxx[2]* u(xg[2],yg[ii]) +
//   cxx[3]* u(xg[3],yg[ii]) + cxx[4]*ub(xg[4],yg[ii]);
//
// ub(xg[0]) and ub(xg[4]) are known boundary conditions, at xg 1,2,3 u unknown
// similarly uyy(xg[i],yg[ii) = cyy[0]*ub(xg[i],yg[0]) +
//     cyy[1]* u(xg[i],yg[1]) + cyy[2]* u(xg[i],yg[2]) +
//     cyy[3]* u(xg[i],yg[3]) + cyy[4]*ub(xg[i],yg[4]);
//
// the product uxx(xg[i],yg[ii])*uyy(xg[i],yg[ii]) has 25 terms
//  4 terms constant cxx[0]*ub(xg[0],yg[ii]) * cyy[0]*ub(xg[i],yg[0])
//                   cxx[4]*ub(xg[4],yg[ii]) * cyy[0]*ub(xg[i],yg[0])
//                   cxx[0]*ub(xg[0],yg[ii]) * cyy[4]*ub(xg[i],yg[4])
//                   cxx[4]*ub(xg[4],yg[ii]) * cyy[4]*ub(xg[i],yg[4])
// 12 terms one unknown
//                   cxx[0]*ub(xg[0],yg[ii]) * cyy[1]* u(xg[i],yg[1])
//          ... 
//                   cxx[4]*ub(xg[4],yg[ii]) * cyy[3]* u(xg[i],yg[3])
//  9 terms product of two unknonws
//                   cxx[1]* u(xg[1],yg[ii]) * cyy[1]* u(xg[i],yg[1])
//          ...
//                   cxx[3]* u(xg[3],yg[ii]) * cyy[3]* u(xg[i],yg[3])
//
// similarly for uxx*uxx and uyy*uxx
 

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "deriv.h"
#include "simeq_newton3.h" // solve system of non linear equations
// also needs invert.h and invert.c 

#undef  abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))
#undef  max
#define max(a,b) ((a)>(b)?(a):(b))

#define nx 7
#define ny 6
// number of equations and number of linear unknowns 
#define neqn (nx-2)*(ny-2)
// number of nonlinear terms, all pairs of unknowns 
#define nlin neqn*neqn
// total number of terms needed for nonlinear solution
#define ntot neqn+nlin

static double A[neqn][ntot]; // nonlinear system of equations 
static double F[neqn];       // RHS of nonlinear equations 
static double xg[nx];        // xgrid 
static double yg[ny];        // ygrid 
static double U[ntot];       // solution being computed and initial 
static double Ua[neqn];      // known solution for testing
static double cxx[nx][nx];   // numerical second derivative i,j 
static double cyy[ny][ny];   // numerical second derivative ii,jj 
static int var1[ntot];       // unknown variable index in U
static int var2[ntot];       // unknown variable index in U second power
static int var3[ntot];       // unknown variable index in U third power
static double xmin = 0.0;
static double xmax = 0.6;
static double ymin = 0.0;
static double ymax = 0.4;
static double hx, hy;
static int check = 0; // 1 for much printout

static double f(double x, double y) // RHS, forcing function 
{
  return 4.0*sin(x-y)*sin(x-y);
} // end f 

// boundary values, in this case, exact solution for checking 
static double ub(double x, double y)
{
  return sin(x-y);
} // end ub 

static int S(int i, int ii) // 0 <= i <= nx, 0 <= ii <= ny
{
  return (i-1)*(ny-2) + (ii-1); // unknown zero based inside boundary
} // end S index into neqn 

static int SS(int unk, int nlk) // unknown i * ii nonlinear term
{
  return neqn+unk*neqn+nlk; // zero based
} // end SS index into ntot 

static void nlterms() // build var1, var2, var3 for simeq_newton3
{
  int i, j, ii, jj, unk, nlk, k;
  // linear, then nonlinear terms neqn, nlin
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      unk = S(i,ii);   // unknown zero based index
      var1[unk] = unk; // linear term
      var2[unk] = -1;
      var3[unk] = -1;
      if(check)
        printf("unk=S(%2d,%2d)=%d, u(xg[%2d],yg[%2d])\n", i, ii, unk, i, ii);
    } // end ii
  } // end i
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      unk = S(i,ii);
      for(j=1; j<nx-1; j++)
      {
        for(jj=1; jj<ny-1; jj++)
	{
          nlk = S(j,jj);
          var1[SS(unk,nlk)] = unk;
          var2[SS(unk,nlk)] = nlk; // product of two terms
          var3[SS(unk,nlk)] = -1;  // no third power
          if(check) printf("unk=%d, nlk=%d, SS(unk,nlk)=%d, u(xg[%2d],yg[%2d])*u(xg[%2d],yg[%2d]) \n",
                            unk, nlk, SS(unk,nlk), i, ii, j, jj);
	} // end jj
      } // end j
    } // end ii
  } // end i
  if(check)
  {
    printf("     var1 var2 var3 \n");
    for(i=0; i<ntot; i++)
    {
      printf("%4d %4d %4d %4d \n", i, var1[i], var2[i], var3[i]);
    }
  } // end check
} // end nlterms

static void initU() // initial guess for non linear solution
{
  int i, j, ii, jj, unk;
  // linear, then nonlinear terms neqn, nlin
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      unk = S(i,ii);
      //U[unk] = ub(xg[i], yg[ii]); // should be fit of boundary
      U[unk] = 0.5; // should be avarage of boundary values
    } // end ii
  } // end i
} // end initU

static void check_row(int row) // needs initU and nlterms
{
  double err = 0.0;
  double val1, val2, val = 0.0;
  int i, j, ii, jj, unk, nlk, nlx;
  // linear, then nonlinear terms neqn, nlin
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      unk = S(i,ii);   // evaluate as though known
      val += A[row][unk]*Ua[unk];
      if(check && (row==0 || row==19))
	printf("unk=%2d, A=%f, val=%f \n", unk, A[row][unk], val);
      // for each linear term, accumulate for nonlinear terms
      for(j=1; j<nx-1; j++)
      {
        for(jj=1; jj<ny-1; jj++)
	{
          nlk = S(j,jj);
          nlx = SS(unk,nlk); // non linear term index in A row
          if(var1[nlx]<0 || var1[nlx]>neqn)
            printf("error var1=%d at nlx=%d \n", var1[nlx], nlx);
          if(var2[nlx]<0 || var2[nlx]>neqn)
            printf("error var2=%d at nlx=%d \n", var2[nlx], nlx);
          if(var3[nlx] != -1)
            printf("error var3=%d at nlx=%d \n", var3[nlx], nlx);
	  val += A[row][nlx]*Ua[var1[nlx]]*Ua[var2[nlx]];
	  if(check && (row==0 || row==19))
            printf("nlx=%d, A[][]=%f, val=%f \n",
		   nlx, A[row][nlx], val);
	} // end jj
      } // end j
    } // end ii
  } // end i
  err = val - F[row];
  printf("check row=%d, F[row]=%f, err=%e \n", row, F[row], err);
  fflush(stdout);
} // end check_row

static void buildA() // A * U = F  U unknown
{
  int i, j, k, ii, jj, kk, row, unk, unk1, unk2, nlx;
  double val, valj, valk, valjj, valkk, txx, txxc, tyy, tyyc;
  // initialize A to zero, F was initialized 
  for(i=0; i<neqn; i++)
  {
    for(j=0; j<ntot; j++)
    {
      A[i][j] = 0.0;
    } // end j
  } // end i

  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      row = S(i,ii); // building a row of A
      printf("buildA working on row %d \n", row);
      F[row] = f(xg[i],yg[ii]);
      val = 0.0;
      // loop through terms in row, 25 txx*txx + 2.0*txx*tyy + tyy*tyy
      for(j=0; j<nx; j++)
      {
	for(k=0; k<nx; k++) // x
        {
          if((j==0 || j==nx-1) && (k==0 || k==nx-1)) // both constant        
	  {
	    valj = cxx[i][j]*ub(xg[j],yg[ii]);
	    valk = cxx[i][k]*ub(xg[k],yg[ii]);
	    F[row] -= valj*valk;
	    if(check && row==0)
	      printf("j=%2d, k=%2d, txx 2con \n", j, k);
	  } // end constant
          else if(j==0 || j==nx-1) // j constant, linear
	  {
	    valj = cxx[i][j]*ub(xg[j],yg[ii]);
	    valk = cxx[i][k];
	    A[row][S(k,ii)] += valj*valk;
	    if(check && row==0)
	      printf("j=%2d, k=%2d, txx jcon \n", j, k);
	  } // end j constant    
          else if(k==0 || k==nx-1) // k constant, linear
	  {
	    valj = cxx[i][j];
	    valk = cxx[i][k]*ub(xg[k],yg[ii]);
	    A[row][S(j,ii)] += valj*valk;
	    if(check && row==0)
	      printf("j=%2d, k=%2d, txx kcon \n", j, k);
	  } // end k constant
          else // both free nonlinear
	  {
	    valj = cxx[i][j];
	    valk = cxx[i][k];
	    A[row][SS(S(j,ii),S(k,ii))] += valj*valk;
	    if(check && row==0)
              printf("j=%2d, k=%2d, txx nlin \n", j, k);
	  } // end if
        } // end k
      } // end j  end txx*txx

      // tyy*tyy
      for(jj=0; jj<ny; jj++)
      {
        for(kk=0; kk<ny; kk++)
	{
	  if((jj==0 || jj==ny-1) && (kk==0 || kk==ny-1)) // both constant     
	  {
	    valjj = cyy[ii][jj]*ub(xg[i],yg[jj]);
	    valkk = cyy[ii][kk]*ub(xg[i],yg[kk]);
	    F[row] -= valjj*valkk;
	    if(check && row==0)
	      printf("jj=%2d, kk=%2d, tyy 2con \n", jj, kk);
	  } // end constant
	  else if(jj==0 || jj==ny-1) // jj constant, linear                   
	  {
	    valjj = cyy[ii][jj]*ub(xg[i],yg[jj]);
	    valkk = cyy[ii][kk];
	    A[row][S(i,kk)] += valjj*valkk;
	    if(check && row==0)
	      printf("jj=&2d, kk=%2d, tyy jjcon\n", jj, kk);
	  } // end j constant
	  else if(kk==0 || kk==ny-1) // kk constant, linear                   
	  {
	    valjj = cyy[ii][jj];
	    valkk = cyy[ii][kk]*ub(xg[i],yg[kk]);
	    A[row][S(i,jj)] += valjj*valkk;
	    if(check && row==0)
	      printf("jj=%2d, kk=%2d, tyy kkcon \n", jj, kk);
	  } // end k constant
	  else // both free nonlinear                                         
	  {
	    valjj = cyy[ii][jj];
	    valkk = cyy[ii][kk];
	    A[row][SS(S(i,jj),S(i,kk))] += valjj*valkk;
	    if(check && row==0)
	      printf("jj=%2d, kk=%2d, tyy nlin \n", jj, kk);
	  } // end if
	} // end kk
      } // end jj tyy*tyy

      // 2.0^txx*tyy
      for(j=0; j<nx; j++)
      {
        for(jj=0; jj<ny; jj++)
	{
	  if((j==0 || j==nx-1) && (jj==0 || jj==ny-1)) // both constant       
	  {
	    // 4 terms constant
	    valj  = cxx[i][j]*ub(xg[j],yg[ii]);
	    valjj = cyy[ii][jj]*ub(xg[i],yg[jj]);
	    F[row] -=  2.0*valj*valjj;
	    if(check && row==0)
	      printf("j=%d, jj=%d, txy 2con \n", j, jj);
	  }
	  else if(j==0 || j==nx-1) // single unknown
	  {
	    valj  = cxx[i][j]*ub(xg[j],yg[ii]);
	    valjj = cyy[ii][jj];
	    A[row][S(i,jj)] += 2.0*valj*valjj;
	    if(check && row==0)
	      printf("j=%2d, jj=%2d, txy jcon \n", j, jj);
	  }
	  else if(jj==0 || jj==ny-1) // single unknown
          {
            valj  = cxx[i][j];
            valjj = cyy[ii][jj]*ub(xg[i],yg[jj]);
            A[row][S(j,ii)] += 2.0*valj*valjj;
            if(check && row==0)
              printf("j=%2d, jj=%2d, txy jjcon \n", j, jj);
          }
	  else
          {
            // three cases of 2 power
            unk1 = S(i,jj);
            unk2 = S(j,ii);
            nlx = SS(unk1,unk2);
            if(check && (row==0 || row==16 || row==19))
              printf("unk1=%2d,unk2=%2d, nlx=%2d \n", unk1, unk2, nlx);
            if(check && (row==0 || row==16))
              printf("i=%2d, ii=%2d, txy nlin \n", i, ii);
            if(check && (row==0 || row==16))
              printf("j=%2d, jj=%2d, txy nlin \n", j, jj);
            valj  = cxx[i][j];
            valjj = cyy[ii][jj];
            A[row][nlx] += 2.0*valj*valjj; // 2 power
	  }
	} // end jj of 2.0*txx*tyy
      } // end j of 2.0*txx*tyy

      printf("finished row %d \n", row);
      fflush(stdout);
      check_row(row);
    } // end ii
  } // end i
  return;
} // end buildA

static void check_soln(double u[]) // check when solution known, full ub 
{
  // uxx*uxx + 2*uxx*uyy + uyy*uyy = f(x,y) check for each 
  // xg, xy inside boundary 
  int i, j, ii, jj;
  double val, txx, tyy, err, smaxerr;

  printf("check_soln, using just actual points inside boundary \n");
  fflush(stdout);
  smaxerr = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      // ynknown = S(i,ii)
      val = 0.0;
      txx = 0.0;
      tyy = 0.0;
      for(j=0; j<nx; j++)
      {
	if(j==0 || j==nx-1) txx += cxx[i][j]*ub(xg[j],yg[ii]);
	else                txx += cxx[i][j]*u[S(j,ii)];
      }
      for(jj=0; jj<ny; jj++)
      {
	if(jj==0 || jj==ny-1) tyy += cyy[ii][jj]*ub(xg[i],yg[jj]);
	else                  tyy += cyy[ii][jj]*u[S(i,jj)];
      }
      val += txx*txx + 2.0*txx*tyy + tyy*tyy; 
      err = val - f(xg[i],yg[ii]);
      smaxerr = max(smaxerr, abs(err));
      printf("i=%2d, ii=%2d, err=%g \n", i,ii, err);
    }
  }
  printf("check_soln max error=%g \n\n", smaxerr);
  fflush(stdout);
} // end check_soln 

int main(int argc, char *argv[])
{
  int i, j, ii, jj, row, unk, nlk;
  double err, avgerr, maxerr;

  printf("pde_nl22.c running second order, second degree nonlinear \n");
  printf("solve uxx*uxx + 2*uxx*uyy + uyy*uyy = f(x,y)\n");
  printf("  f(x,y) = 4*sin(x-y)^2 \n");
  printf(" \n");
  printf("  0<x<1  0<y<1 \n");
  printf("boundary and Analytic solution ub(x,y) = sin(x-y) \n");
  printf("solve for U A[%d][%d]*U=F /n", neqn, ntot);

  hx = (xmax-xmin)/(double)(nx-1);
  printf("nx=%d, xmin=%f, xmax=%f, hx=%f \n", nx, xmin, xmax, hx);
  for(i=0; i<nx; i++)
  {
    xg[i] = xmin+i*hx; // does not need to be uniform 
    printf("xg[%1d]=%f \n", i, xg[i]);
  } // end i

  hy = (ymax-ymin)/(double)(ny-1);
  printf("ny=%d, ymin=%f, ymax=%f, hy=%f \n", ny, ymin, ymax, hy);
  for(ii=0; ii<ny; ii++)
  {
    yg[ii] = ymin+ii*hy; // does not need to be uniform 
    printf("yg[%1d]=%f \n", ii, yg[ii]);
  } // end ii
  printf("RHS, F forcing function and exact solution \n");
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      row = S(i,ii);
      Ua[row] = ub(xg[i],yg[ii]);
      F[row]  =  f(xg[i],yg[ii]);
      printf("i=%d, ii=%d, F(%7.3f,%7.3f)=%7.3f, exact=%7.3f \n",
	     i, ii, xg[i], yg[ii], F[row], Ua[row]);
    } // end ii 
  } // end i
  printf("\n");

  printf("build cxx[][] and cyy[][] for multiple use \n");
  for(i=0; i<nx; i++)
  {
    nuderiv(2, nx, i, xg, cxx[i]);
    for(j=0; j<nx; j++)
    {
      printf("cxx[%1d][%1d]=%f \n", i, j, cxx[i][j]);
    } // end j
  } // end i
  printf(" \n");
  for(ii=0; ii<ny; ii++)
  {
    nuderiv(2, ny, ii, yg, cyy[ii]);
    for(jj=0; jj<ny; jj++)
    {
      printf("cyy[%1d][%1d]=%f \n", ii, jj, cyy[ii][jj]);
    } // end j
  } // end i
  printf(" \n");

  nlterms(); // build var1, var2, var3
  initU();   // initial guess, starting point, for non linear solution
             // used by check_row in buildA
  
  printf("build A, non linear system of equations A * U = F \n");
  buildA();

  printf("solve non linear equations A * U = F\n");
  fflush(stdout);
  simeq_newton3(neqn, nlin, A, F, var1, var2, var3, U, 14, 0.001, 1);

  // optional solution test
  printf("check solution A*Ua=F against known solution \n");
  check_soln(Ua);
  printf("check solution A*U=F against computed solution \n");
  check_soln(U);


  // check accuracy, print results
  maxerr = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(ii=1; ii<ny-1; ii++)
    {
      row = S(i,ii);
      err = U[row]-Ua[row];
      maxerr = max(maxerr, abs(err));
      printf("i=%2d, ii=%2d, x=%7.3f, y=%7.3f, Ua=%7.3f, U=%7.3f, err=%g\n",
             i, ii, xg[i], yg[ii], Ua[row], U[row], err);
    } // end ii
  } // end i
  printf("pde_nl22.c finished, maxerr=%e \n", maxerr);
  return 0;
} // end main of pde_nl22.c