// pde_nl22.java  non linear PDE, second order, second degree, two dimension
//
// The PDE to be solved is:
//  a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) + 
//  c1(x,y)*u(x,y)uyy(x,y) = f(x,y)
//
//  f(x,y) = 0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)*
//          (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)*
//          (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
//          5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+
//          (8.0-2.0*exp(x)-2.0*exp(0.5*y))*
//          (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
//          5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x)
//
//  a1(x,y) = exp(x/2)*exp(y)/2
//  b1(x,y) = 0.7/(x^2*y^2+0.5)
//  c1(x,y) = 2*(4-exp(x)-exp(y/2))
//
// Boundary conditions computed from
// ub(x,y) = x^4 + 2 y^4 + 3 x^3 y + 4 x y^3 +
//           5 x^2 y^2 + 6 x^2 + 7 y^2 + 8
//

import java.io.*;
import java.text.*;

class pde_nl22
{
  // number of points for discretization         
  int nx = 5;
  int ny = 5;
  int dof = (nx-2)*(ny-2);  // Degrees of Freedom
  int nlin = dof*dof; // all possible, non boundary combinations
  int ntot = dof+nlin; // total number of terms, linear and nonlinear
  double hx, hy;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double ug[][] = new double[nx][ny]; // computed solution, incl boundary
  double uA[][] = new double[nx][ny]; // actual solution
  double us[] = new double[ntot]; // single subscript solution
  int var1[]  = new int[ntot];
  int var2[]  = new int[ntot];
  int var3[]  = new int[ntot];
  int vari1[] = new int[ntot];
  int vari2[] = new int[ntot];
  simeq_newton5 SN = new simeq_newton5(1.0e-3, 1.0, 24, false);

  // Define the region for the solution:                               
  // xmin <= x <= xmax, ymin <= y <= ymax  using  nx by ny points      

  double xmin = -1.0;
  double xmax =  1.0;
  double ymin = -1.0;
  double ymax =  1.0;

  // coefficients of derivatives
  double cx[][] = new double[nx][nx];
  double cxx[][] = new double[nx][nx];
  double cy[][] = new double[ny][ny];
  double cyy[][] = new double[ny][ny];
  int k = 0;
  DecimalFormat f4 = new DecimalFormat("0.000");
  double start_time, end_time;
  boolean debug = false;
  String Xname[];

  pde_nl22() // construct and run
  {
    double x, y;

    System.out.println("pde_nl22.java running");
    System.out.println("nonlinear second order, second degree, two dimension");
    System.out.println("a1(x,y)*uxx(x,y)*uyy(x,y) + b1(x,y)*u(x,y)*uxx(x,y) +");
    System.out.println("c1(x,y)*u(x,y)uyy(x,y) = f(x,y)      ");
    System.out.println(" ");
    System.out.println("f(x,y)=0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)*");
    System.out.println("       (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)*");
    System.out.println("       (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+");
    System.out.println("       5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+");
    System.out.println("       (8.0-2.0*exp(x)-2.0*exp(0.5*y))*");
    System.out.println("       (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+");
    System.out.println("       5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x);");
    System.out.println("a1(x,y) = exp(x/2.0)*exp(y)/2.0");
    System.out.println("b1(x,y) = 0.7/(x*x*y*y+0.5)");
    System.out.println("c1(x,y) = (4.0 - exp(x) - exp(y/2.0))*2.0");
    System.out.println(" ");
    System.out.println("Boundary values ");
    System.out.println("ub(x,y)= x^4 + 2 y^4 + 3 x^3 y + 4 x y^3 +");
    System.out.println("         5 x^2 y^2 + 6 x^2 + 7 y^2 + 8");
    System.out.println(" ");
    start_time = System.currentTimeMillis()/1000.0;

    // initialize hx, hy, xg, yg
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    System.out.println("xmin="+xmin+", xmax="+xmax+", hx="+hx+", nx="+nx);
    System.out.println("ymin="+ymin+", ymax="+ymax+", hy="+hy+", ny="+ny);
    for(int i=0; i<nx; i++) xg[i] = xmin+hx*(double)i;
    for(int ii=0; ii<ny; ii++) yg[ii] = ymin+hy*(double)ii;

    // initialize ug with boundary values, uA all values
    System.out.println("Dirchilet Boundary values on square");
    System.out.println("x             y           ub(x,y)");
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
        y = yg[ii];
        uA[i][ii] = ub(x,y);
        System.out.println("xg="+f4.format(x)+", yg="+f4.format(y)+
                           ", uA["+i+"]["+ii+"]="+f4.format(uA[i][ii])+
                           ", f(x,y)="+f4.format(f(x,y)));
        if(i==0 || ii==0 || i==nx-1 || ii==ny-1)
	{
          ug[i][ii] = ub(x,y);
	}
	else
	{
          ug[i][ii] = 0.0;
	}
      }
    }

    // initialize nonlinear system of equations
    System.out.println("dof="+dof+", nlin="+nlin+", dof+nlin="+(dof+nlin));

    for(int i=1; i<nx-1; i++) // linear variables
    {
      for(int ii=1; ii<ny-1; ii++)
      {
	k = S(i,ii);
	var1[k] = k; // linearized double subscript
	var2[k] = -1;
	var3[k] = -1;
	vari1[k] = -1;
	vari2[k] = -1;
      }
    }
    for(int i=0; i<dof; i++) // nonlinear terms
    {
      for(int j=0; j<dof; j++)
      {
	k = dof+i*dof+j; // a bit redundant, coefficients will be zero
        var1[k] = i;
	var2[k] = j;
	var3[k] = -1;
	vari1[k] = -1;
	vari2[k] = -1;
      }
    }

    // for development and debug print
    Xname = new String[dof+nlin];

    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        k = S(i,ii);
        Xname[k] = new String("X"+Integer.toString(i)+Integer.toString(ii));
        // System.out.println("Xname["+k+"]= "+Xname[k]);
      }
    }
    for(int j=0; j<dof; j++)
    {
      for(int jj=0; jj<dof; jj++)
      {
	k = dof + j*dof+jj;
	Xname[k] = new String(Xname[j]+"*"+Xname[jj]);
        // System.out.println("Xname["+k+"]= "+Xname[k]);
      }
    } // end for development and debug print used below

    print_eqn(dof, nlin, var1, var2, var3, vari1, vari2);

    // compute coefficients of derivatives
    for(int i=0; i<nx; i++)
    {
      // compute derivative coefficients in x direction
      new rderiv(1,nx,i,hx,cx[i]);    // Ux
      new rderiv(2,nx,i,hx,cxx[i]);   // Uxx
    }
    for(int ii=0; ii<ny; ii++)
    {
      // compute derivative coefficients in y direction
      new rderiv(1,ny,ii,hy,cy[ii]);    // Uy
      new rderiv(2,ny,ii,hy,cyy[ii]);   // Uyy
    }
    if(debug)
    {
      for(int i=0; i<nx; i++)
        System.out.println("cx["+i+"] "+cx[i][0]+","+cx[i][1]+","+
                           cx[i][2]+","+cx[i][3]+","+cx[i][4]);
      for(int i=0; i<nx; i++)
        System.out.println("cxx["+i+"] "+cxx[i][0]+","+cxx[i][1]+","+
                           cxx[i][2]+","+cxx[i][3]+","+cxx[i][4]);
      for(int ii=0; ii<ny; ii++)
        System.out.println("cy["+ii+"] "+cy[ii][0]+","+cy[ii][1]+","+
                            cy[ii][2]+","+cy[ii][3]+","+cy[ii][4]);
      for(int ii=0; ii<ny; ii++)
        System.out.println("cyy["+ii+"] "+cyy[ii][0]+","+cyy[ii][1]+","+
                           cyy[ii][2]+","+cyy[ii][3]+","+cyy[ii][4]);
    }

    solve_pde();
    for(int i=1; i<nx-1; i++)
      for(int ii=1; ii<ny-1; ii++)
	System.out.println("at "+i+","+ii+" computed "+f4.format(us[S(i,ii)])+
                           "  exact="+f4.format(uA[i][ii])+"  error="+
                           (uA[i][ii]-us[S(i,ii)]));
    System.out.println(" ");

    check_soln(ug);
    System.out.println(" ");

    System.out.println("check check_soln using expected solution");
    check_soln(uA);

    System.out.println(" ");
    end_time = System.currentTimeMillis()/1000.0;
    System.out.println("pde_nl22.java finished in "+(end_time-start_time)+
                       " seconds");
  } // end pde_nl22
 
  double a1(double x, double y)
  {
    return Math.exp(x/2.0)*Math.exp(y)/2.0;
  }

  double b1(double x, double y)
  {
    return 0.7/(x*x*y*y+0.5);
  }

  double c1(double x, double y)
  {
    return (4.0 - Math.exp(x) - Math.exp(y/2.0))*2.0;
  }

  // must compute boundary.
  // also solution for checking 
  double ub(double x, double y)
  {
    return x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*y +
           6.0*x + 7.0*y + 8.0;
  }

  // RHS, for testing, this becomes defined by u and a1, b1, c1  

  double f(double x, double y)
  {

    return  0.5*Math.exp(x/2.0)*Math.exp(y)*(6.0*x+6.0*y)*
            (12.0*y+8.0*x)+0.7/(x*x*y*y+0.5)*
            (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
            5.0*x*y+6.0*x+7.0*y+8.0)*(6.0*x+6.0*y)+
            (8.0-2.0*Math.exp(x)-2.0*Math.exp(0.5*y))*
            (x*x*x+2.0*y*y*y+3.0*x*x*y+4.0*x*y*y+
	    5.0*x*y+6.0*x+7.0*y+8.0)*(12.0*y+8.0*x);
  } // end f

  void solve_pde()
  {
    double A[][] = new double[dof][dof+nlin];
    double Y[] = new double[dof];
    double x, y, a1xy, b1xy, c1xy, t;
    int row;

    // clear A
    for(int i=0; i<dof; i++)
    {
      for(int j=0; j<dof+nlin; j++)
      {
	A[i][j] = 0.0;
      }
    }

    for(int i=1; i<nx-1; i++)  // i and j are x direction
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++) // ii and jj are y direction
      {
	y = yg[ii];
        row = S(i, ii); // equation index, row
        a1xy = a1(x,y);
        b1xy = b1(x,y);
        c1xy = c1(x,y);
        Y[row] = f(x,y);

        // discretize each term for each equation
        // term  a1(x,y) * uxx * uyy
        for(int j=0; j<nx; j++) // x direction
	{
	  if(j==0 || j==nx-1) // j boundary
	  {
	    for(int jj=0; jj<ny; jj++) // y direction
	    {
	      if(jj==0 || jj==ny-1) // j and jj boundary 
	      {
		// ug is set for boundary conditions
                t = a1xy*(cxx[i][j]*ug[j][ii])*(cyy[ii][jj]*ug[i][jj]);
		Y[row] -= t;
	      }
              else // j boundary, jj not boundary
	      {
		k = S(i,jj);
                t = a1xy*cyy[ii][jj]*(cxx[i][j]*ug[j][ii]);
                A[row][k] += t;
	      }
	    } // end jj
	  }
          else // j non boundary
	  {
            for(int jj=0; jj<ny; jj++)
	    {
	      if(jj==0 || jj==ny-1) // j not boundary, jj boundary
	      {
		k = S(j,ii);
                t = a1xy*cxx[i][j]*(cyy[ii][jj]*ug[i][jj]);
                A[row][k] += t;
	      }
	      else // neither boundary (non linear both)
	      {
		  k = dof+dof*S(j,ii)+S(i,jj); // nonlinear position
                t = a1xy*cxx[i][j]*cyy[ii][jj];
                A[row][k] += t;
	      }
	    } // end jj
	  } // end if
	} // end j

        // term  b1(x,y) * u * uxx
        for(int j=0; j<nx; j++) // x direction
	{
	  if(j==0 || j==nx-1) // j boundary
	  {
	    // ug is set for boundary conditions
	    k = S(i,ii);
            t = b1xy*(cxx[i][j]*ug[j][ii]);
	    A[row][k] += t;
	  }
          else // j not boundary
	  {
	    k = dof + dof*S(i,ii) +S(j,ii);
            t = b1xy*cxx[i][j];
            A[row][k] += t;
	  } // end if
	} // end j

        // term  c1(x,y) * u * uyy
        for(int jj=0; jj<ny; jj++)
	{
	  if(jj==0 || jj==ny-1) // j boundary
	  {
	    // ug is set for boundary conditions
	    k = S(i,ii);
            t = c1xy*(cyy[ii][jj]*ug[i][jj]);
	    A[row][k] += t;
	  }
          else // jj not boundary
	  {
	    k = dof + dof*S(i,ii) + S(i,jj);
            t = c1xy*cyy[ii][jj];
            A[row][k] += t;
	  } // end if
	} // end jj

      } // end y
    } // end x

    System.out.println(" ");
    System.out.println("A matrix, and Y RHS, matrix values of zero not printed");
    for(int i=0; i<dof; i++)
    {
      for(int ii=0; ii<ntot; ii++)
	if(A[i][ii]!=0.0)
	  System.out.println("A["+i+"]["+ii+"]="+A[i][ii]+"  at "+Xname[ii]);
      System.out.println("Y["+i+"]="+Y[i]);
      System.out.println(" ");
    }
    // testing guess, should be very close
    for(int i=1; i<nx-1; i++)
      for(int ii=1; ii<ny-1; ii++)
	us[S(i,ii)] = uA[i][ii];
    System.out.println("Check by giving solution as initial guess");
    SN.simeq(dof, nlin, A, Y, var1, var2, var3, vari1, vari2, us);
    System.out.println("no output expected from SN.simeq ");
    for(int i=1; i<nx-1; i++)
      for(int ii=1; ii<ny-1; ii++)
	System.out.println("at "+i+","+ii+" expected "+f4.format(uA[i][ii])+
                             "  got "+f4.format(us[S(i,ii)]));

    // initial guess
    for(int i=0; i<dof; i++) us[i] = 5.0;
    System.out.println("initial guess, all="+us[0]);

    SN.simeq(dof, nlin, A, Y, var1, var2, var3, vari1, vari2, us);
    for(int i=1; i<nx-1; i++)
      for(int ii=1; ii<ny-1; ii++) 
	  ug[i][ii] = us[S(i,ii)]; // final solution of linear variables
    System.out.println(" ");

  } // end solve_pde

  int S(int i, int ii)
  {
    return (i-1)*(ny-2)+(ii-1); // linear index, not including boundary
  }

  void check_soln(double Ux[][]) // check when solution unknown
  {
    double err, sumerr, sumsqerr, rmserr, maxerr, avgerr;
    double x, y, u, uxx, uyy, eval_deriv, term; 
    int cnt;

    System.out.println("check solution against PDE");
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;

    for(int i=1; i<nx-1; i++) // inside boundary
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
        y = yg[ii];

        // compute derivatives numerically for this x,y
        u = Ux[i][ii]; 
        uxx = 0.0;
        for(int j=0; j<nx; j++)
	{
          uxx  += cxx[i][j]*Ux[j][ii]; 
	}
        uyy = 0.0;
        for(int jj=0; jj<ny; jj++)
        {
	  uyy  += cyy[ii][jj]*Ux[i][jj];
	}

        eval_deriv = 0.0; // based on numeric derivatives using ug[][]
        // term   a1(x,y)*uxx(x,y)*uyy(x,y) +
        term  = a1(x,y)*uxx*uyy;
        eval_deriv += term;
	// term   b1(x,y)*u(x,y)*uxx(x,y) + 
        term = b1(x,y)*u*uxx; 
        eval_deriv += term;
        // term   c1(x,y)*u(x,y)uyy(x,y) = f(x,y)      
        term = c1(x,y)*u*uyy;
        eval_deriv += term;
        err = eval_deriv - f(x,y);
        err = Math.abs(err);
        sumerr += err;
        sumsqerr += err*err;
        maxerr = Math.max(err,maxerr);
        cnt++;  
      } // end y
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln 

  void print_eqn(int n, int nlin, int var1[], int var2[],
                 int var3[], int vari1[], int vari2[])
  {
    String Xname[] = new String[n];

    System.out.println("solve system of equations A X = Y for X");
    System.out.println("the equations are for i=0,"+(n-1));
    for(int i=1; i<nx-1; i++)
    {
      for(int ii=1; ii<ny-1; ii++)
      {
        k = S(i,ii);
        Xname[k] = new String("X"+Integer.toString(i)+Integer.toString(ii));
        System.out.println("A[i]["+k+"]*"+Xname[k]+"+");
      }
    }
    for(int i=n; i<n+nlin; i++)
    {
      if(var1[i]>=n || var2[i]>=n || var3[i]>=n ||
         vari1[i]>=n || vari2[i]>=n)
	  System.out.println("var error "+var1[i]+" "+var2[i]+" "+var3[i]+
                             " "+vari1[i]+" "+vari2[i]);
      System.out.print("A[i]["+i+"]");
      if(var1[i]>=0) System.out.print("*"+Xname[var1[i]]);
      if(var2[i]>=0) System.out.print("*"+Xname[var2[i]]);
      if(var3[i]>=0) System.out.print("*"+Xname[var3[i]]);
      if(vari1[i]>=0 || vari2[i]>=0) System.out.print("/(");
      if(vari1[i]>=0) System.out.print(Xname[vari1[i]]);
      if(vari2[i]>=0) System.out.print("*"+Xname[vari2[i]]);
      if(vari1[i]>=0 || vari2[i]>=0) System.out.print(")");
      if(i==n+nlin-1) System.out.println(" ");
      else       System.out.println("+");
    }
    System.out.println("  = Y[i]");
    System.out.println(" ");
  } // end print_eqn

  public static void main(String args[])
  {
    new pde_nl22();
  } // end main

} // end class pde_nl22