import javax.swing.*;
import org.math.plot.*;
import static java.lang.Math.*;
import static org.math.array.DoubleArray.*;
import visad.*;
import visad.java3d.DisplayImplJ3D;
import visad.java2d.DisplayImplJ2D;
import java.rmi.RemoteException;
import java.awt.*;
import javax.swing.*;
import java.awt.event.*;

public class triangular_grid_xyt
{
public triangular_element te[];
public setElement se[];
public int P,N,M,K,K1,K2,Nm1,Mm1;
public double x[][];
public double y[][];
public double T[][];
public double Ti[];
public double xi[];
public double yi[];
public boolean bi[];
public double dx,dy;
public int Bx0[],BxM[],By0[],ByN[];
public double L,W;

public triangular_grid_xyt(double Li, double Wi,int Ni, int Mi,double Ti[][],double dt,double t)
{generateTriangularGrid(Li,Wi,Ni,Mi,Ti,dt,t);}

public triangular_grid_xyt(double xi[][],double yi[][],double Ti[][],double dt,double t)
{ x=xi;
  y=yi;
  T=Ti;
  N=x.length-1;
  M=x[0].length-1;
  L=max(xi);
  W=max(yi);
  setElements1();
}
public void generateTriangularGrid(double Li, double Wi,int Ni, int Mi,double zi[][],double dt,double t)
{int i,j;
 N=Ni;M=Mi;
 L=Li;
 W=Wi;
 x=new double[N+1][M+1];
 y=new double[N+1][M+1];
 T=zi;
 setElements();
}

public void generateTriangularGrid(double xi[][],double yi[][],double zi[][],double dt,double t)
{ x=xi;
  y=yi;
  T=zi;
  N=x.length-1;
  M=x[0].length-1;
  L=max(xi);
  W=max(yi);
  setElements1();
}

public void setElements()
{ int i,j;
 dx=L/M;
 dy=W/N;
 K=M*N;
 K1=K*2;
 K2=K1+1;
 //grid points
 xi=new double[K2];
 yi=new double[K2];
 Ti=new double[K2];
 bi=new boolean[K2];
 se=new setElement[K2];
 for(i=0;i<K2;i++) {se[i]=new setElement();}
 te=new triangular_element[K1]; 
 int k1=-1,k=0,n=0,m=0;
 for(n=0;n<=N;n++)
 { for(m=0;m<=M;m++)
  {k=(N+1)*n+m;
   x[n][m]=m*dx;y[n][m]=n*dy;
   xi[k]=x[n][m];
   yi[k]=y[n][m];
   Ti[k]=T[n][m];
   bi[k]=false;
  }
}
for(n=0;n<N;n++)
{for(m=0;m<M;m++)
  { k=(N+1)*n+m;
	int kp1=(N+1)*(n+1)+m;
	int kp1p1= kp1+1; 
	double xi1[]={x[n][m],x[n][m+1],x[n+1][m]};
	double yi1[]={y[n][m],y[n][m+1],y[n+1][m]};
    double zi1[]={T[n][m],T[n][m+1],T[n+1][m]};
	k1++; 
	int n_point1[]={k,(k+1),kp1};
	se[k].add(k1,0);
	se[k+1].add(k1,1);
	se[kp1].add(k1,2);
	te[k1]=new triangular_element(k1,n_point1,xi1,yi1);  	
	double xi2[]={x[n][m+1],x[n+1][m+1],x[n+1][m]};
	double yi2[]={y[n][m+1],y[n+1][m+1],y[n+1][m]};
	double zi2[]={T[n][m+1],T[n+1][m+1],T[n+1][m]};
	te[k1].setT(zi1);
	k1++;
	int n_point2[]={k+1,kp1p1,kp1};
    se[k+1].add(k1,0);
	se[kp1p1].add(k1,1);
	se[kp1].add(k1,2);
	te[k1]=new triangular_element(k1,n_point2,xi2,yi2);
	te[k1].setT(zi2);
  } 
}	
}
public void setElements1()
{ int i,j;
 dx=L/M;
 dy=W/N;
 K=M*N;
 K1=K*2;
 K2=K1+1;
 //grid points
 xi=new double[K2];
 yi=new double[K2];
 bi=new boolean[K2];
 se=new setElement[K2];
 for(i=0;i<K2;i++) {se[i]=new setElement();}
 te=new triangular_element[K1]; 
 int k1=-1,k=0,n=0,m=0;
 for(n=0;n<=N;n++)
 { for(m=0;m<=M;m++)
  {k=(N+1)*n+m;
   xi[k]=x[n][m];
   yi[k]=y[n][m];
   Ti[k]=T[n][m];
   bi[k]=false;
  }
}
for(n=0;n<N;n++)
{for(m=0;m<M;m++)
  { k=(N+1)*n+m;
	int kp1=(N+1)*(n+1)+m;
	int kp1p1= kp1+1; 
	double xi1[]={x[n][m],x[n][m+1],x[n+1][m]};
	double yi1[]={y[n][m],y[n][m+1],y[n+1][m]};
    double zi1[]={T[n][m],T[n][m+1],T[n+1][m]};
	k1++; 
	int n_point1[]={k,(k+1),kp1};
	se[k].add(k1,0);
	se[k+1].add(k1,1);
	se[kp1].add(k1,2);
	te[k1]=new triangular_element(k1,n_point1,xi1,yi1);  	
	double xi2[]={x[n][m+1],x[n+1][m+1],x[n+1][m]};
	double yi2[]={y[n][m+1],y[n+1][m+1],y[n+1][m]};
	double zi2[]={T[n][m+1],T[n+1][m+1],T[n+1][m]};
	te[k1].setT(zi1);
	k1++;
	int n_point2[]={k+1,kp1p1,kp1};
    se[k+1].add(k1,0);
	se[kp1p1].add(k1,1);
	se[kp1].add(k1,2);
	te[k1]=new triangular_element(k1,n_point2,xi2,yi2);
	te[k1].setT(zi2);
  } 
}	
}

public double min(double xi[][])
{ double min=1e200;
  for(int i=0;i<xi.length;i++)
   {for(int j=0;j<xi[0].length;j++)
     {if(xi[i][j]<min) min=xi[i][j];}
 }
 return min;
}
public double max(double xi[][])
{double max=-1e200;
 for(int i=0;i<xi.length;i++)
 {for(int j=0;j<xi[0].length;j++)
   {if(xi[i][j]<max) max=xi[i][j];}
 }
 return max;
}
public int find_element(double z[])
{ return find_element(z[0],z[1]);}
public int find_element(double x,double y)
{double R1=Math.sqrt(x*x+y*y);
 int kmin=0;
 double min=1.0e50;
 for(int k=0;k<K1;k++)
 {   te[k].setN(x,y);
     double r0=Math.sqrt((te[k].x[0]-x)*(te[k].x[0]-x)+(te[k].y[0]-y)*(te[k].y[0]-y));
     double r1=Math.sqrt((te[k].x[1]-x)*(te[k].x[1]-x)+(te[k].y[1]-y)*(te[k].y[1]-y));
     double r2=Math.sqrt((te[k].x[2]-x)*(te[k].x[2]-x)+(te[k].y[2]-y)*(te[k].y[2]-y));
     double dr=Math.sqrt(r0*r0+r1*r1+r2*r2);
	 boolean b1=(x>=te[k].x[0]) && (x<=te[k].x[2]);
	 boolean b2=(y>=te[k].y[0]) && (y<=te[k].y[1]);
	 boolean b= b1 && b2;
     if(dr<min) {min=dr;kmin=k;}
 }
 return kmin;
}
public void setConvection(int cond,double L0,double L,double W0,double W, double hi, double tinfi)
{//cond 0: x==L0
 //cond 1: x==L
 //cond 2: y==W0
 //cond 3: y==W
 int j=0;
 for(int i=0;i<K1;i++)
 {  j=te[i].setI(cond,L0,L,W0,W);
	te[i].setConvection(j,hi,tinfi);
 }
}
public void setHeatFlux(int cond,double L0,double L,double W0,double W, double qi)
{//cond 0: x==L0
 //cond 1: x==L
 //cond 2: y==W0
 //cond 3: y==W
 int j=-1;
 for(int i=0;i<K1;i++)
 {  j=te[i].setI(cond,L0,L,W0,W);
	te[i].setHeatFlux(j,qi);
 }
}
public void setTemperature(int cond,double L0,double L,double W0,double W, double TTi)
{//cond 0: x==L0
 //cond 1: x==L
 //cond 2: y==W0
 //cond 3: y==W
 for(int i=0;i<K2;i++)
 {  if(cond==0)
    {if(xi[i]==L0) {Ti[i]=TTi;bi[i]=true;}}
    else if(cond==1)
    {if(xi[i]==L) {Ti[i]=TTi;bi[i]=true;}}
    else if(cond==2)
    {if(yi[i]==W0) {Ti[i]=TTi;bi[i]=true;}}
    else if(cond==3)
    {if(yi[i]==W) {Ti[i]=TTi;bi[i]=true;}}
 }
}
public void setHeatGeneration(int node[],double qgeni)
{  int n=node.length;
   for(int i=0;i<n;i++)
   {te[node[i]].setHeatGeneration(qgeni);}
}

public static double[][] TT1(int N,int M,double tt[])
{
int k1=0;
double TTT[][]=new double[N][M];
for(int n=0;n<N;n++)
{ for(int m=0;m<M;m++)
  {k1=(N+1)*n+m;
  TTT[n][m]=tt[k1];
  }
}
return TTT;  
}

public static double[] TT2(double tt[][])
{
	int N=tt.length;
	int M=tt[0].length;
	int K1=N*M;
	System.out.println("N="+N+"M="+M+"K1="+K1);
	double TT[]=new double[K1];
    int k1=0;
    for(int n=0;n<N;n++)
    { for(int m=0;m<M;m++)
       {k1=N*n+m;
        TT[k1]=tt[n][m];
       }
    }
return TT;  
}
public double[][] solve(double fi,double ti,double k[],double roi,double Cpi)
{
double K[][]=new double[K2][K2];
double KK[][]=new double[K2][K2];
double ff[]=new double[K2];
double KN[][]=new double[3][3];
double f[]=new double[K2];
double fN[]=new double[3];
for(int i=0;i<K1;i++)
{KN=te[i].AA(fi,ti,k,roi,Cpi);
 fN=te[i].BB(fi,ti,k,roi,Cpi);
 int n1=te[i].n_element;
 int m[]=te[i].n_point;
 for(int j=0;j<3;j++)
   { for(int k1=0;k1<3;k1++)
      {K[m[j]][m[k1]]+=KN[j][k1];KK[m[j]][m[k1]]=K[m[j]][m[k1]];}
       f[m[j]]+=fN[j];       
       ff[m[j]]=f[m[j]];
   }
}
//redesign K matrix and f vector if any known 
//temperature boundary is given
//int j=0;
for(int j=0;j<K2;j++)
{   
	for(int i=0;i<K2;i++)
     { if(bi[i])
       {
	   if(i==j)
	   {f[i]=Ti[i];K[i][j]=1.0;}
	   else
	   {
	   K[i][j]=0.0;}
	   //System.out.println("i="+i+"j="+j+"bi="+bi[i]+"Ti="+Ti[i]+"fi="+f[i]);
       }
       }//end of for
} //end of for
//System.out.println("K =\n"+Matrix.toString(K));
//System.out.println("f =\n"+Matrix.toStringT(f));
double TT[]=AXB(K,f);
//System.out.println("TT =\n"+Matrix.toStringT(TT));

int k1=0;
for(int n=0;n<=N;n++)
{ for(int m=0;m<=M;m++)
  {k1=(N+1)*n+m;
  T[n][m]=TT[k1];
  }
}
return T; 
}
	
//LU decomposition method
public static double[][] LU(double c[][],int indx[],int d[])
{
//returns LU decomposition of matrix c and index indx
double a[][];
int n=c.length;
a=new double[n][n];
double vv[];
vv=new double[n];
double sum,dum,big,temp;
int i,j,k;
int imax;
int nmax=100;
double tiny=1.0e-40;
imax=0;
for(i=1;i<=n;i++)
{
  for(j=1;j<=n;j++)
     a[i-1][j-1]=c[i-1][j-1];
}
d[0]=1;
for(i=1;i<=n;i++)
  {
  big=0.0;
  for(j=1;j<=n;j++)
    {
    if(Math.abs(a[i-1][j-1])>big) big=Math.abs(a[i-1][j-1]);
    }
  if(big==0) {System.out.println("singular matrix");return a;}
  vv[i-1]=1.0/big;
  }
for(j=1;j<=n;j++)
{
  for(i=1;i<j;i++)
    {
      sum=a[i-1][j-1];
      for(k=1;k<i;k++)
	{
	  sum-=a[i-1][k-1]*a[k-1][j-1];
	}
    a[i-1][j-1]=sum;
    }
  big=0;
  for(i=j;i<=n;i++)
    {
    sum=a[i-1][j-1];
    for(k=1;k<j;k++)
      {
      sum-=a[i-1][k-1]*a[k-1][j-1];
      }
    a[i-1][j-1]=sum;
    dum=vv[i-1]*Math.abs(sum);
    if(dum>=big)
      {
      imax=i;
      big=dum;
      }
    } //end of i=0
if(j != imax)
    {
    for(k=1;k<=n;k++)
      {
      dum=a[imax-1][k-1];
      a[imax-1][k-1]=a[j-1][k-1];
      a[j-1][k-1]=dum;
      }
    d[0]=-d[0];
    vv[imax-1]=vv[j-1];
    } //end of if
indx[j-1]=imax;
if(a[j-1][j-1]==0) a[j-1][j-1]=tiny;
if(j!=n)
  {
  dum=1.0/a[j-1][j-1];
  for(i=j+1;i<=n;i++)
    a[i-1][j-1]*=dum;
  }//endif
} //end for j=
return a;
}

public static double[] LUaxb(double a[][],double x[],int indx[])
{
//solves AX=B system of linear equation of LU decomposed matrix a
//(calculated by method LU)
int ii=0;
int i,j,ll=0;
double sum=0;
int n=a.length;
double b[];
b=new double[n];
for(i=1;i<=n;i++)
{
    b[i-1]=x[i-1];
}
for(i=1;i<=n;i++)
  {
  ll=indx[i-1];
  sum=b[ll-1];
  b[ll-1]=b[i-1];
  if(ii!=0)
     {
     for(j=ii;j<=(i-1);j++)
       {
       sum-=a[i-1][j-1]*b[j-1];
       }
     }
  else if(sum!=0) ii=i;
  b[i-1]=sum;
  }
for(i=n;i>=1;i--)
  {
  sum=b[i-1];
  if(i<n)
    {
    for(j=(i+1);j<=n;j++)
      {
      sum-=a[i-1][j-1]*b[j-1];
      }
    }
  b[i-1]=sum/a[i-1][i-1];
  }
return b;
}
public void ContourPlot()
{
ContourPlot pp=new ContourPlot( x,y,T,16,true,"Finite element","x","y");
pp.plot();  
}
public void Plot3D(double L0,double L,double W0,double W) throws RemoteException, VisADException
{
Plot3D.plot3D_CD1(T,L0,L,W0,W,200,"x","y","T degree C");  
}

public void GridPlot(double L0,double L,double W0,double W)
{Plot3DPanel plot = new Plot3DPanel("SOUTH");
 double dx=(L-L0)/M;
 double dy=(W-W0)/N;
 double[] xx = increment(L0, dx,L); // x = 0.0:0.1:1.0
 double[] yy = increment(W0,dy,W);
 plot.addGridPlot("Finite element Grid", xx, yy, T);
 JFrame frame = new JFrame("a plot panel");
 frame.setSize(600, 600);
 frame.setContentPane(plot);
 frame.setVisible(true); 
}

public static double[] AXB(double a[][],double b[])
{
//Solution of system of linear equations by LU method
// note that the same calculation can be done by divide method.
int n=a.length;
double c[]=new double[n];
int d[]={1};
int indx[]=new int[n];
double e[][]=new double[n][n];
e=LU(a,indx,d);
c=LUaxb(e,b,indx);
return c;
} //end of AXB
public static void main(String arg[]) throws RemoteException, VisADException
{double L=1.0;
 double W=1.0;
 int N=5,M=5;
 double hi=100.0;
 double dt=1.0;
 double t=2;
 double ti[][]=new double[N+1][M+1];
 for(int i=0;i<=N;i++)
  {for(int j=0;j<=M;j++)
    {ti[i][j]=100.0;}
  }
 //double ti[][]={{20.0,20.0,20.0},{20.0,20.0,20.0},{20.0,20.0,20.0}};
 double tti[]=triangular_grid_xyt.TT2(ti);
 System.out.println("n="+tti.length);
 triangular_grid_xyt tg=new triangular_grid_xyt(L,W,M,N,ti,dt,t);
 System.out.println("x=\n"+Matrix.toString(tg.x));
  System.out.println("y=\n"+Matrix.toString(tg.y));
   System.out.println("T=\n"+Matrix.toString(tg.T));
 int cond0=0;//x==0
 int cond1=1;//x==L
 int cond2=2;//y==0
 int cond3=3;//y==W
 double Ti0=80.0;
 double Ti1=80.0;
 double Ti2=80.0;
 double Ti3=200.0;
 double qi1=0.0;
 
  //tg.setHeatFlux(cond2,0.0,L,0.0,W,qi);
  //tg.setHeatFlux(cond3,0.0,L,0.0,W,qi);
  //tg.setConvection(cond0,0.0,L,0.0,W,hi,Ti0);
  tg.setTemperature(cond0,0.0,L,0.0,W,Ti0);
  tg.setTemperature(cond1,0.0,L,0.0,W,Ti1);
  tg.setTemperature(cond2,0.0,L,0.0,W,Ti2);
  tg.setTemperature(cond3,0.0,L,0.0,W,Ti3);
  //tg.setConvection(cond1,0.0,L,0.0,W,hi,Ti1);
 //tg.setTemperature(cond2,0.0,L,0.0,W,Ti2);
 //tg.setConvection(cond3,0.0,L,0.0,W,hi,Ti3);
 //tg.setTemperature(cond3,0.0,L,0.0,W,Ti3);
 //tg.setConvection(cond3,0.0,L,0.0,W,hi,Ti3);
 //tg.setConvection(cond1,0.0,L,0.0,W,hi,Ti2);
 //tg.setConvection(cond3,0.0,L,0.0,W,hi,Ti3);
 //tg.setTemperature(cond1,0.0,L,0.0,W,Ti1);
 //tg.setTemperature(cond2,0.0,L,0.0,W,Ti2);
 //tg.setTemperature(cond3,0.0,L,0.0,W,Ti3);
 //tg.setConvection(cond3,0.0,L,0.0,W,hi,Ti3);
 //tg.setHeatFlux(cond3,0.0,L,0.0,W,qi);
 double k[]={50,50};
 double fi=0.5;
 double ro=7800.0;
 double Cp=1200.0;
 for(double tt=0;tt<=t;tt+=dt)
  {double TT[][]=tg.solve(fi,t,k,ro,Cp);
   System.out.println("t="+tt+"\n"+Matrix.toString(TT));
    tg.setTemperature(cond0,0.0,L,0.0,W,Ti0);
    tg.setTemperature(cond1,0.0,L,0.0,W,Ti1);
    tg.setTemperature(cond2,0.0,L,0.0,W,Ti2);
    tg.setTemperature(cond3,0.0,L,0.0,W,Ti3);
  }
 
 /*Text.print(tg.T);
 tg.ContourPlot();
 tg.Plot3D(0.0,L,0.0,W);
 tg.GridPlot(0.0,L,0.0,W);
 */
}
}