// fem_nl12_la.java  Galerkin FEM
// solve D*ux(x) + E(x)*(-ux(x)/u(x)^2 - F*uxx(x) = f(x)
//       f(x) = -2
//       D = 1
//       F = 1
//       E(x) = (x^2+1)^2
//       Analytic solution u(x)=x^2+1
//       inside grid, Gauss Legendre integration

// uses simeq.java
// uses nderiv.java
// uses laphi.java
// uses gaulegf.java
import java.text.*;

class fem_nl12_la
{
  // i, j, nx, xg needed by galk or galf, available at call 
  final int nx = 10;
  double xg[] = new double[nx];
  double hx;
  int i, j;
  laphi L = new laphi();
  nderiv RD = new nderiv();
  double D = 1.0;
  double F = 1.0;
  DecimalFormat f8 = new DecimalFormat("   0.000");

  fem_nl12_la()
  {
    double x;
    final int npx = 96; // Gauss Legendre integration order
    double xx[] = new double[npx+1];
    double wx[] = new double[npx+1]; // for Gauss-Legendre 
    double val, err, avgerr, maxerr;
    int px;

    double xmin = 1.0; // problem parameters 
    double xmax = 2.0;

    double kg[][] = new double[nx][nx];  
    double fg[] = new double[nx];
    double ug[] = new double[nx];
    double Ua[] = new double[nx];


    System.out.println("fem_nl12_la.java running  ");
    System.out.println("Given:");
    System.out.println("D*ux(x)+E(x)*(d/dx 1/u(x)-F*uxx(x)=f(x)");
    System.out.println("D*ux(x)-E(x)*(ux(x)/u(x)^2)-F*uxx(x)=f(x)");
    System.out.println("f(x)=-2");
    System.out.println("D=1, F=1, E(x)=(x^2+1)^2");
    System.out.println("xmin<=x<=xmax  Boundaries  ");
    System.out.println("Analytic solution u(x)=x^2 + 1");

    System.out.println("xmin="+xmin+", xmax="+xmax);
    System.out.println("nx="+nx);
    System.out.println("x grid and analytic solution ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx; // does not need to be uniform
      Ua[i] = uana(xg[i]); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    }  
    System.out.println(" ");

    j=1; x=xg[j];
    System.out.println("observing Lagrange Phi functions,j="+j);
    System.out.println("L.phi(x,j,nx-1,xg)="+L.phi(x,j,nx-1,xg));
    System.out.println("L.phip(x,j,nx-1,xg)="+L.phip(x,j,nx-1,xg));
    System.out.println("L.phipp(x,j,nx-1,xg)="+L.phipp(x,j,nx-1,xg));
    j=2; x=xg[j];
    System.out.println("observing Lagrange Phi functions, j="+j);
    System.out.println("L.phi(x,j,nx-1,xg)="+L.phi(x,j,nx-1,xg));
    System.out.println("L.phip(x,j,nx-1,xg)="+L.phip(x,j,nx-1,xg));
    System.out.println("L.phipp(x,j,nx-1,xg)="+L.phipp(x,j,nx-1,xg));
    j=nx-2; x=xg[j];
    System.out.println("observing Lagrange Phi functions, j=nx-2");
    System.out.println("L.phi(x,j,nx-1,xg)="+L.phi(x,j,nx-1,xg));
    System.out.println("L.phip(x,j,nx-1,xg)="+L.phip(x,j,nx-1,xg));
    System.out.println("L.phipp(x,j,nx-1,xg)="+L.phipp(x,j,nx-1,xg));
    System.out.println(" ");

    System.out.println("calling gauleg xmin="+xmin+", xmax="+xmax+
                       ", npx="+npx);
    new gaulegf(xmin, xmax, xx, wx, npx); // xx and wx set up for integrations 
    System.out.println("xx[1]="+xx[1]+", xx[2]="+xx[2]);
    System.out.println("xx[3]="+xx[3]+", xx[4]="+xx[4]+", xx[npx]="+xx[npx]);
    System.out.println("wx[1]="+wx[1]+", wx[2]="+wx[2]);
    System.out.println("wx[3]="+wx[3]+", wx[4]="+wx[4]+", wx[npx]="+wx[npx]);
    i=1; j=1;
    val = galk(xx[2]);
    System.out.println("galk(xx[2])="+val);
    val = galf(xx[2]);
    System.out.println("galf(xx[2])="+val);
    System.out.println(" ");


    // put solution in Ua vector 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      Ua[i] = uana(x);
      System.out.println("solution at i="+i+",x="+x+
                           " is "+Ua[i]);
    }
    System.out.println(" ");

    // initialize k and f 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        kg[i][j] = 0.0;
      }
    }
    for(i=0; i<nx; i++)
    {
        fg[i] = 0.0;
    }

    // set boundary conditions 
    i = 0;
    x = xg[i];
    kg[i][i] = 1.0;
    fg[i] = uana(x);
    System.out.println("boundary i="+i+" is "+fg[i]);
    i=nx-1;
    x = xg[i];
    kg[i][i] = 1.0;
    fg[i] = uana(x);
    System.out.println("boundary i="+i+" is "+fg[i]);

    // compute entries in stiffness matrix 
    System.out.println("compute stiffness matrix  ");
    for(i=1; i<nx-1; i++)
    {
      for(j=0; j<nx; j++)
      {
        // compute integral for k(i,j)  
        val = 0.0;
        for(px=1; px<=npx; px++)
        {
          val = val + wx[px]*galk(xx[px]);
	}
        kg[i][j] = val;
      } // end j loop 
    } // end i loop 

    // compute integral for f(i)  
    for(i=1; i<nx-1; i++)
    {
      val = 0.0;
      for(px=1; px<=npx; px++)
      {
        val = val + wx[px]*galf(xx[px]);
      }
      System.out.println("Legendre integration="+val+", f at i="+i);
      fg[i] = val;
    }
    System.out.println("k computed stiffness matrix, see above  ");
    System.out.println("f computed forcing function, see above  ");
    new simeq(kg, fg, ug);

    System.out.println(" ");
    System.out.println("check_soln on known solution");
    check_soln(Ua);
    System.out.println("check_soln based on numeric approximation of PDE");
    check_soln(ug);
    System.out.println(" ");

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed Galerkin, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      err = Math.abs(ug[i]-Ua[i]);
      if(err>maxerr) maxerr = err;
      avgerr = avgerr + err;
      System.out.println("ug["+i+"]="+ug[i]+", Ua="+Ua[i]+
                           ", err="+(ug[i]-Ua[i]));
    }
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx));
    System.out.println(" ");
  } // end fem_nl12_la constructor

  // PDE problem definition functions
  double f(double x) // RHS
  {
    return -2.0;
  }

  double uana(double x) // u 
  {
     return x*x+1.0;
  }

  double E(double x)  
  {
     return (x*x+1.0)*(x*x+1.0);
  }

  double galk(double x)
  {                        // Galerkin k stiffness function for this problem 
    return (D*L.phip(x,j,nx-1,xg)
            -E(x)*L.phip(x,j,nx-1,xg)/(L.phi(x,j,nx-1,xg)*L.phi(x,j,nx-1,xg))
	    -F*L.phipp(x,j,nx-1,xg))
            *L.phi(x,i,nx-1,xg);
  } // end galk used by integration 

  double galf(double x) // Galerkin f function for this problem 
  {
    return f(x)*L.phi(x,i,nx-1,xg);
  } // end galf used by integration 

  void check_soln(double U[])
  {
    // check D*ux(x) + E(x)*(-ux(x)/u(x)^2) - F*uxx(x) = f(x)
    double maxerr, err, sumerr, sumsqerr, rmserr, avgerr;
    double cx[] = new double[nx];
    double cxx[] = new double[nx];
    double eval_deriv, x;
    int cnt;

    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    maxerr = 0.0;

    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      // compute derivative coefficients in x direction
      RD.rderiv(1,nx,i,hx,cx);
      RD.rderiv(2,nx,i,hx,cxx);

      eval_deriv = 0.0;

      for(int j=0; j<nx; j++)
      {
	eval_deriv += D*cx[j]*U[j]; // U must include boundary values
        eval_deriv -= E(x)*cx[j]*U[j]/(U[i]*U[i]);
        eval_deriv -= F*cxx[j]*U[j];
      }
     
      err = eval_deriv-f(x);  
      sumerr += Math.abs(err);
      sumsqerr += err*err;
      maxerr = Math.max(Math.abs(err),maxerr);
      cnt++;
    }
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln

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