// fem2D.java      problem definition and Galerkin phi functions:
//                elliptic equation 
//                Dirichlet boundary conditions
//                uxx(x,y) + uyy(x,y) + r(x,y)*u(x,y) = f(x,y)
//                u(x,y)=x*x+2*x*y+3*y*y+4*x+5*y+6
//                r(x,y)=x+2*y   requires integration be cubic
//                f(x,y)=8+(x+2*y)*(x*x+2*x*y+3*y*y+4*x+5*y+6)
//
//  uses phi2D and
//  provides intgk(T1,T2), intgf(T1) and bound(x,y) for solver

public class fem2D
{
  fem2D()
  {
    System.out.println("fem2D called");
    System.out.println("uxx(x,y)+uyy(x,y)+r(x,y)*u(x,y)=f(x,y)");
    System.out.println("u(x,y)=x*x+2*x*y+3*y*y+4*x+5*y+6");
    System.out.println("r(x,y)=x+2*y");
    System.out.println("f(x,y)=8+(x+2*y)*(x*x+2*x*y+3*y*y+4*x+5*y+6)");
    new triquad();
    new phi2D();
  } // end fem2D constructor

  // functions specific to this problem

  public static double f(double x, double y)
  {
      return 8.0+(x+2.0*y)*(x*x+2.0*x*y+3.0*y*y+4.0*x+5.0*y+6.0);
  }

  public static double r(double x, double y)
  {
    return x+2.0*y;
  }

  public static double bound(double x, double y) 
  {
      return x*x+2.0*x*y+3.0*y*y+4.0*x+5.0*y+6.0;
  }

  public static double galk(double TT[], double x, double y)
  {
    return phi2D.phixx(TT,x,y)+phi2D.phiyy(TT,x,y)+r(x, y)*phi2D.phi(TT,x,y);
  }

  public static double intgk(double TT1[],double TT2[])
  {
    // integrate galk(TT1) * phi(TT2) 
    int np=6;
    int nn;
    double X[] = new double[np*np];
    double Y[] = new double[np*np];
    double W[] = new double[np*np];

    nn = triquad.tint(np,TT1,X,Y,W);
    double u=0;
    for(int i=0; i<nn; i++)
       u=u+W[i]*galk(TT1,X[i],Y[i])*phi2D.phi(TT2,X[i],Y[i]);
    return u;
  }

  public static double intgf(double TT2[])
  {
    // integrate f * phi(TT2)
    int np=6;
    int nn;
    double X[] = new double[np*np];
    double Y[] = new double[np*np];
    double W[] = new double[np*np];

    nn = triquad.tint(np,TT2,X,Y,W);
    double u=0;
    for(int i=0; i<nn; i++)
	u = u + W[i]*f(X[i],Y[i])*phi2D.phi(TT2,X[i],Y[i]);
    return u;
  }
} // end class fem2D