// pde_nl22c.java  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
 
public class pde_nl22c
{
  int nx = 7; // problem parameters
  int ny = 6;
  // number of equations and number of linear unknowns 
  int neqn = (nx-2)*(ny-2);
  // number of nonlinear terms, all pairs of unknowns 
  int nlin = neqn*neqn;
  // total number of terms needed for nonlinear solution
  int ntot = neqn+nlin;

  double A[][] = new double[neqn][ntot]; // nonlinear system of equations 
  double F[]  =  new double[neqn];       // RHS of nonlinear equations 
  double xg[] =  new double[nx];         // xgrid 
  double yg[] =  new double[ny];         // ygrid 
  double U[]  =  new double[neqn+nlin];  // solution being computed, big 
  double Ua[] =  new double[neqn];  // solution known, single variable 
  double cxx[][] = new double[nx][nx];   // numerical second derivative i,j 
  double cyy[][] = new double[ny][ny];   // numerical second derivative ii,jj 
  int var1[] = new int[ntot];    // unknown variable index in U
  int var2[] = new int[ntot];    // unknown variable index in U second power
  int var3[] = new int[ntot];    // unknown variable index in U third power
  double var[] = new double[neqn]; // for initial guess, middle of range
  double xmin = 0.0;
  double xmax = 0.6;
  double ymin = 0.0;
  double ymax = 0.4;
  double hx, hy;
  boolean check = true; // set true for debugging

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

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

  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 

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

  void nlterms() // build var1, var2, var3 for simeq_newton3
  {
    int unk, nlk;
    // linear, then nonlinear terms neqn, nlin
    System.out.println(neqn+" linear terms, then "+nlin+" nonlinear terms");
    System.out.println("unk = S(i,ii)");
    for(int i=1; i<nx-1; i++)
    {
      for(int 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) System.out.println("unk=S("+i+","+ii+")= "+unk+
                                     ", u(xg["+(i)+"],yg["+(ii)+"])");
      } // end ii
    } // end i
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        unk = S(i,ii);   // unknown zero based index

        for(int j=1; j<nx-1; j++)
        {
          for(int 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) System.out.println("unk="+unk+", nlk="+nlk+
                             ", SS(unk,nlk)="+SS(unk,nlk)+
                             ", u(xg["+(i)+"],yg["+(ii)+"])*u(xg["+
                             j+"],yg["+jj+"])");
	  } // end jj
        } // end j
      } // end ii
    } // end i
    if(check)
    {
      System.out.println("     var1 var2 var3 ");
      for(int i=0; i<ntot; i++)
      {
        System.out.println("i="+i+"  "+var1[i]+" "+var2[i]+" "+var3[i]);
      }
    } // end check
  } // end nlterms

  void initU() // initial guess for non linear solution
  {
    int unk;
    // linear, then nonlinear terms neqn, nlin
    for(int i=1; i<nx-1; i++)
    {
      for(int 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 average of boundary values
      } // end ii
    } // end i
  } // end initU

  void check_row(int row) // needs initU and nlterms
  {
    double err = 0.0;
    double val1, val2, val = 0.0;
    int unk, nlk, nlx;
    // linear, then nonlinear terms neqn, nlin
    for(int i=1; i<nx-1; i++)
    {
      for(int 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))
          System.out.println("unk="+unk+", A="+A[row][unk]+", val="+val);
        // for each linear term, accumulate for nonlinear terms
        for(int j=1; j<nx-1; j++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
	    nlk = S(j,jj);     // non linear term
            nlx = SS(unk,nlk); // non linear term index in A row
            if(var1[nlx]<0 || var1[nlx]>neqn)
              System.out.println("error var1="+var1[nlx]+" at nlx="+ nlx);
            if(var2[nlx]<0 || var2[nlx]>neqn)
              System.out.println("error var2="+var2[nlx]+" at nlx="+ nlx);
            if(var3[nlx] != -1)
               System.out.println("error var3="+var3[nlx]+" at nlx="+nlx);
	    val += A[row][nlx]*Ua[var1[nlx]]*Ua[var2[nlx]]; 
            if(check && (row==0 || row==19))
              System.out.println("nlx="+nlx+", A="+A[row][nlx]+", val="+val);
	  } // end jj
        } // end j
      } // end ii
    } // end i
    err = val - F[row];
    System.out.println("check row="+row+", F[row]="+F[row]+", err="+err);
  } // end check_row

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

    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        row = S(i,ii); // building a row of A
        System.out.println("buildA working on row "+row+", i="+i+", ii="+ii);
        F[row] = f(xg[i],yg[ii]);
        val = 0.0;
        // loop through terms in row, 25  txx*txx + 2.0*txx*tyy + tyy*tyy 
        // txx*txx
        for(int j=0; j<nx; j++) // x
        {
	  for(int 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) 
                System.out.println("j="+j+", k="+k+", txx 2con");
	    } // 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) 
                System.out.println("j="+j+", k="+k+", txx jcon");
	    } // 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)
                System.out.println("j="+j+", k="+k+", txx kcon");
	    } // 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)
                System.out.println("j="+j+", k="+k+", txx nlin");
	    } // end if
	  } // end k
	} // end j  end txx*txx

        // tyy*tyy
        for(int jj=0; jj<ny; jj++)
	{
	  for(int 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)
                System.out.println("jj="+jj+", kk="+kk+", tyy 2con");
	    } // 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)
                System.out.println("jj="+jj+", kk="+kk+", tyy jjcon");
	    } // 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) 
                System.out.println("jj="+jj+", kk="+kk+", tyy kkcon");
	    } // 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)
                System.out.println("jj="+jj+", kk="+kk+", tyy nlin");
	    } // end if
	  } // end kk
	} // end jj tyy*tyy

        // 2.0^txx*tyy
        for(int j=0; j<nx; j++)
        {
          for(int 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)
                System.out.println("j="+j+", jj="+jj+", txy 2con");
	    }
            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)
                System.out.println("j="+j+", jj="+jj+", txy jcon");
	    }
            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)
                System.out.println("j="+j+", jj="+jj+", txy jjcon");
	    }
            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))
                System.out.println("unk1="+unk1+",unk2="+unk2+", nlx="+nlx);
              if(check && (row==0 || row==16)) 
                System.out.println("i="+i+", ii="+ii+", txy nlin");
              if(check && (row==0 || row==16))
                System.out.println("j="+j+", jj="+jj+", txy nlin");
              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

        System.out.println("finished row "+row);
        check_row(row);
      } // end ii
    } // end i
  } // end buildA

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

    smaxerr = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        // check at = S(i,ii)
        val = 0.0;
        txx = 0.0;
        tyy = 0.0;
        for(int 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(int 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 = Math.max(smaxerr, Math.abs(err));
        System.out.println("i="+i+", ii="+ii+", err="+err);
      } // end ii on Y
    } // end i on X
    System.out.println("check_soln max error="+smaxerr);
  } // end check_soln 

  public pde_nl22c()
  {
    int row, unk, nlk;
    double err, avgerr, maxerr;

    System.out.println("pde_nl22c.java running second order, second degree nonlinear ");
    System.out.println("solve uxx*uxx + 2*uxx*uyy + uyy*uyy = f(x,y)");
    System.out.println("  f(x,y) = 4*sin(x-y)^2 ");
    System.out.println(" ");
    System.out.println("  0<x<1  0<y<1 ");
    System.out.println("boundary and Analytic solution ub(x,y) = sin(x-y) ");
    System.out.println("solve for U A["+neqn+"]["+ntot+"]*U=F");

    hx = (xmax-xmin)/(double)(nx-1);
    var[0] = (xmax-xmin)/2.0;
    System.out.println("nx="+nx+", xmin="+xmin+", xmax="+xmax+", hx="+hx);
    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx; // does not need to be uniform 
      System.out.println("xg["+i+"]="+xg[i]);
    } // end i

    hy = (ymax-ymin)/(double)(ny-1);
    var[1] = (ymax-ymin)/2.0;
    System.out.println("ny="+ny+", ymin="+ymin+", ymax="+ymax+", hy="+hy);
    for(int ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy; // does not need to be uniform 
      System.out.println("yg["+ii+"]="+yg[ii]);
    } // end ii
    System.out.println("RHS, F forcing function and exact solution ");
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        row = S(i,ii);
        Ua[row] = ub(xg[i],yg[ii]);
        F[row]  =  f(xg[i],yg[ii]);
        if(check)
          System.out.println("i="+i+", ii="+ii+", F("+xg[i]+","+yg[ii]+")="+
                             F[row]+", exact="+Ua[row]);
      } // end ii 
    } // end i
    System.out.println("");

    System.out.println("build cxx[][] and cyy[][] for multiple use ");
    for(int i=0; i<nx; i++)
    {
      new nuderiv(2, nx, i, xg, cxx[i]);
      if(check) for(int j=0; j<nx; j++)
      {
        System.out.println("cxx["+i+"]["+j+"]="+cxx[i][j]);
      } // end j
    } // end i
    System.out.println(" ");
    for(int ii=0; ii<ny; ii++)
    {
      new nuderiv(2, ny, ii, yg, cyy[ii]);
      if(check) for(int jj=0; jj<ny; jj++)
      {
        System.out.println("cyy["+ii+"]["+jj+"]="+cyy[ii][jj]);
      } // end j
    } // end i
    System.out.println(" ");

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

    System.out.println("solve non linear equations A * U = F");
    new simeq_newton3(neqn, nlin, A, F, var1, var2, var3, U, 14, 1.0, 0.0001, 2);

    // optional solution test
    System.out.println("check solution A*Ua=F against known solution");
    check_soln(Ua);
    System.out.println("check solution A*U=F against computed solution");
    check_soln(U);

    // check accuracy, print results
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        row = S(i,ii);
        err = U[row]-Ua[row];
        maxerr = Math.max(maxerr, Math.abs(err));
        System.out.println("i="+i+", ii="+ii+", x="+xg[i]+", y="+yg[ii]);
        System.out.println("Ua="+Ua[row]+", U="+U[row]+", err="+err);
      } // end ii
    } // end i
    System.out.println("pde_nl22c.java finished, maxerr="+maxerr);
  } // pde_nl22c 

  public static void main (String[] args)
  {
    new pde_nl22c(); // construct and execute
  }
} // end class pde_nl22c