public class triangular_element
{
public double AA[][],BB[];
public double t;         //Thickness of two dimensional element m
public double a[],b[],c[];
public double dt;
public int P;
public double x[],y[];   //coordinates of data points m
public double xc,yc;     //center point m
public double A;         //Area m^2
public double N[],B[][]; // Shape factor and derivative
public double Ti[];      // Temperature nodes
public int n_point[];    // point numbers
public int n_element;    // element number
public double l[];       // triangular element side length m
public boolean isconv[];
public boolean isheatflux[];
public boolean isheatgen;
public double h[];       // convective heat transfer coefficient W/m^2K
public double Tinf[];    // environmental temperature degree C (degree K)
public double qflux[];   // heat flux W/m^2
public double qgen;      // heat generation W/m^3
public double ro;        // density kg/m^3
public double Cp;        // Specific heat KJ/kgK
public triangular_element(int n_elementi,int n_pointi[],double xi[],double yi[])
{dt=1.0;
 P=1;
 x=new double[3];
 y=new double[3];
 Ti=new double[3];
 l=new double[3];//0:lij,1:ljk,2:lik
 isconv=new boolean[3];
 isheatflux=new boolean[3];
 n_point=new int[3];
 xc=0;yc=0;
 isheatgen=false;
 h=new double[3];
 Tinf=new double[3];
 qflux=new double[3];
 n_element=n_elementi;
 for(int i=0;i<3;i++)
 {n_point[i]=n_pointi[i];
  x[i]=xi[i];y[i]=yi[i];
  xc+=x[i];yc+=y[i];
  isconv[i]=false;
  isheatflux[i]=false;
  h[i]=0.0;
  Tinf[i]=0.0;
  qflux[i]=0.0;
  qgen=0.0;
 }
 xc/=3.0;
 yc/=3.0;
 a=new double[3];
 b=new double[3];
 c=new double[3];
 N=new double[3];
 B=new double[2][3];
 a[0]=x[1]*y[2]-x[2]*y[1];
 a[1]=x[2]*y[0]-x[0]*y[2];
 a[2]=x[0]*y[1]-x[1]*y[0];
 b[0]=y[1]-y[2];
 b[1]=y[2]-y[0];
 b[2]=y[0]-y[1];
 c[0]=x[2]-x[1];
 c[1]=x[0]-x[2];
 c[2]=x[1]-x[0];
 A=0.5*((x[0]*y[1]-x[1]*y[0])+(x[2]*y[0]-x[0]*y[2])+(x[1]*y[2]-x[2]*y[1]));
 for(int i=0;i<3;i++)
 { B[0][i]=1.0/(2.0*A)*b[i];
   B[1][i]=1.0/(2.0*A)*c[i];	 
 }
 l[0]=Math.sqrt(b[2]*b[2]+c[2]*c[2]);
 l[1]=Math.sqrt(b[0]*b[0]+c[0]*c[0]);
 l[2]=Math.sqrt(b[1]*b[1]+c[1]*c[1]);
}
public triangular_element(int n_elementi,int n_pointi[],double xi[],double yi[],double dti,int Pi)
{
 x=new double[3];
 y=new double[3];
 Ti=new double[3];
 l=new double[3];//0:lij,1:ljk,2:lik
 isconv=new boolean[3];
 isheatflux=new boolean[3];
 n_point=new int[3];
 xc=0;yc=0;
 dt=dti;
 P=Pi;
 isheatgen=false;
 h=new double[3];
 Tinf=new double[3];
 qflux=new double[3];
 n_element=n_elementi;
 for(int i=0;i<3;i++)
 {n_point[i]=n_pointi[i];
  x[i]=xi[i];y[i]=yi[i];
  xc+=x[i];yc+=y[i];
  isconv[i]=false;
  isheatflux[i]=false;
  h[i]=0.0;
  Tinf[i]=0.0;
  qflux[i]=0.0;
  qgen=0.0;
 }
 xc/=3.0;
 yc/=3.0;
 a=new double[3];
 b=new double[3];
 c=new double[3];
 N=new double[3];
 B=new double[2][3];
 a[0]=x[1]*y[2]-x[2]*y[1];
 a[1]=x[2]*y[0]-x[0]*y[2];
 a[2]=x[0]*y[1]-x[1]*y[0];
 b[0]=y[1]-y[2];
 b[1]=y[2]-y[0];
 b[2]=y[0]-y[1];
 c[0]=x[2]-x[1];
 c[1]=x[0]-x[2];
 c[2]=x[1]-x[0];
 A=0.5*((x[0]*y[1]-x[1]*y[0])+(x[2]*y[0]-x[0]*y[2])+(x[1]*y[2]-x[2]*y[1]));
 for(int i=0;i<3;i++)
 { B[0][i]=1.0/(2.0*A)*b[i];
   B[1][i]=1.0/(2.0*A)*c[i];	 
 }
 l[0]=Math.sqrt(b[2]*b[2]+c[2]*c[2]);
 l[1]=Math.sqrt(b[0]*b[0]+c[0]*c[0]);
 l[2]=Math.sqrt(b[1]*b[1]+c[1]*c[1]);
}

public double[][] C(double ti,double roi,double Cpi)
{Cp=Cpi;
ro=roi;
return C(ti);
}

public double[][] C(double ti)
{t=ti;
 double a0=ro*Cp*A*t/12.0;
 double C1[][]={{4.0*a0,a0,a0},{a0,4.0*a0,a0},{a0,a0,4.0*a0}};
return C1;
}

public double[][] K(double ti,double k[])
{t=ti;
 double a0=t/(4.0*A)*k[0];
 double a1=t/(4.0*A)*k[1];
 double b0=t*h[0]*l[0]/6.0;
 double b1=t*h[1]*l[1]/6.0;
 double b2=t*h[2]*l[2]/6.0;
 double K1[][]={{(a0*b[0]*b[0]+a1*c[0]*c[0]+2*b0+2.0*b2),
                 (a0*b[0]*b[1]+a1*c[0]*c[1]+b0),
                 (a0*b[0]*b[2]+a1*c[0]*c[2]+b2)},
                {(a0*b[1]*b[0]+a1*c[1]*c[0]+b0),
                 (a0*b[1]*b[1]+a1*c[1]*c[1]+2*b0+2*b1),
                 (a0*b[1]*b[2]+a1*c[1]*c[2]+b1)},
                {(a0*b[2]*b[0]+a1*c[2]*c[0]+b2),
                 (a0*b[2]*b[1]+a1*c[2]*c[1]+b1),
                 (a0*b[2]*b[2]+a1*c[2]*c[2]+2*b1+2*b2)}
                 };
return K1;
}
public double[][] AA(double fi,double ti,double k[],double roi,double Cpi)
{ double B1[][]=Matrix.multiply(K(ti,k),dt*(1.0-fi));
  double B2[][]=Matrix.add(B1,C(ti,roi,Cpi));
  return B2;
}

public double[] BB(double fi,double ti,double k[],double roi,double Cpi)
{ double B1[][]=Matrix.multiply(K(ti,k),(-dt*fi));
  double B2[][]=Matrix.add(B1,C(ti,roi,Cpi));
  double B3[]=Matrix.multiply(B2,Ti);
  double B4[]=Matrix.multiply(dt,f(ti,k));
  double B5[]=Matrix.add(B3,B4);
  return B5;
}

public double[] f(double ti,double k[])
{
t=ti;
double a0=qgen*A*t/3.0;
double b0=-qflux[0]*t*l[0]/2.0;
double b1=-qflux[1]*t*l[1]/2.0;
double b2=-qflux[2]*t*l[2]/2.0;
double c0= h[0]* Tinf[0]*t*l[0]/2.0;
double c1= h[1]* Tinf[1]*t*l[1]/2.0;
double c2= h[2]* Tinf[2]*t*l[2]/2.0;
double f1[]={(a0+b0+b2+c0+c2),(a0+b0+b1+c0+c1),(a0+b1+b2+c1+c2)};	
return f1;
}

public int setI(int cond,double L0,double L,double W0,double W)
{  
   int i=100;
	if(cond==0) //x=0
	{if(x[0]==L0 && x[1]==L0) i=0;
	 else if(x[1]==L0 && x[2]==L0) i=1;
	 else if(x[0]==L0 && x[2]==L0) i=2;
	 else i=-1;
	}
	else if(cond==1) //x=L
	{if(x[0]==L && x[1]==L) i=0;
	 else if(x[1]==L && x[2]==L) i=1;
	 else if(x[0]==L && x[2]==L) i=2;
	 else i=-1;
	}
	else if(cond==2) //y=0
	{if(y[0]==W0 && y[1]==W0) i=0;
	 else if(y[1]==W0 && y[2]==W0) i=1;
	 else if(y[0]==W0 && y[2]==W0) i=2;
	 else i=-1;
	}
	else if(cond==3) //y=W
	{if(y[0]==W && y[1]==W) i=0;
	 else if(y[1]==W && y[2]==W) i=1;
	 else if(y[0]==W && y[2]==W) i=2;
	 else i=-1;
	}
	return i;	
}

public void setConvection(int i, double hi, double tinfi)
{if(i>=0 && i<=2) 
{  
	h[i]=hi;Tinf[i]=tinfi;isconv[i]=true;
}
}
public void setRo(double roi)
{ro=roi;}
public void setdT(double dti)
{dt=dti;}
public void setCp(double Cpi)
{Cp=Cpi;}

public void setHeatFlux(int i, double qi)
{if(i>=0 && i<=2 ){qflux[i]=qi;isheatflux[i]=true;}	}
	
public void setHeatGeneration(double qgeni)
{qgen=qgeni;isheatgen=true;}
public void setN(double x,double y)
{for(int i=0;i<3;i++)
  {N[i]=1.0/(2.0*A)*(a[i]+b[i]*x+c[i]*y);}
}
public double[]  q(double k)
{ double qi[]=new double[2];
	for(int i=0;i<3;i++)
  {qi[0]+=-k*B[0][i]*Ti[i];
   qi[1]+=-k*B[1][i]*Ti[i];
  }
  return qi;
}

public double[]  q(double k,double Ti[])
{  //heat transfer
	setT(Ti);
   double qi[]=new double[2];
   for(int i=0;i<3;i++)
     {qi[0]+=-k*B[0][i]*Ti[i];
      qi[1]+=-k*B[1][i]*Ti[i];
     }
  return qi;
}

public void setN(double z[])
{setN(z[0],z[1]);}
public void setT(double TT[])
{for(int i=0;i<3;i++)
  {Ti[i]=TT[i];}
}
public double T(double x,double y)
{ //interpolation
	double TT=0;
 setN(x,y);
 for(int i=0;i<3;i++)
  {TT+=N[i]*Ti[i];}
   return TT;
}

public double T(double T[],double x,double y)
{ //set grid point temperatures and get interpolated value
	setT(T);
  setN(x,y);
  double TT=0;
 for(int i=0;i<3;i++)
  {TT+=N[i]*Ti[i];}
  return TT;
}
public double T(double T[],double z[])
{ return T(T,z[0],z[1]);}

public static void main(String arg[])
{int n_elementi=0;
int n_pointi[]={0,1,2};
double xi[]={0.0,4,0};
double yi[]={0,0,2.5};
double dt=1.0;
int P=100;
triangular_element fe=new triangular_element(n_elementi,n_pointi,xi,yi,dt,P);
double p[]={2.0,1.0};
fe.setN(p);
double Ti[]={50,70,100};
fe.setT(Ti);
double T=fe.T(Ti,p[0],p[1]);
double q[]=fe.q(200.0);
System.out.println("T="+T);
}

}