//3次元四面体要素　非定常　CG法　Ｌ２ノルム　    
public class tet_CG {
	static double Lx = 1.0; //X方向の長さ
    static double Ly = 1.0; //Y方向の長さ
    static double Lz = 1.0; //Z方向の長さ
    static int Nx = 8; //X方向の分割数
    static int Ny = 8; //Y方向の分割数
    static int Nz = 8; //Z方向の分割数
    static double sx = -Lx/2;
    static double sy = -Ly/2;
    static double sz = -Lz/2;
    
    static double dx = Lx/Nx; //X方向の微小長さ
    static double dy = Ly/Ny; //Y方向の微小長さ
    static double dz = Lz/Nz; //Z方向の微小長さ
    static int Nn = ((Nx+1)*(Ny+1)*(Nz+1)); //節点数
    static double []x = new double[Nn+1];
    static double []y = new double[Nn+1];
    static double []z = new double[Nn+1];
    int ii = 0;
   
    static int NN = Nx*Ny*Nz; 
    static int N = NN*6; //要素数
    
    static int [][]ie = new int[N+1][4+1];
    
    static double [] V = new double[N+1];
	static double Vn;
	static double []X1 = new double[5];
	static double []Y1 = new double[5];
	static double []Z1 = new double[5];
	static int [] ia = new int[5];
	
	static int I = 0;
	static int J = 0;
	static int K = 0;
	static int L = 0;
	
	static double [][]M = new double[5][5];
	static double [][]SM = new double[Nn+1][Nn+1];
	static double [][]SMM = new double[Nn+1][Nn+1];
	static double [][]MM = new double[Nn+1][Nn+1];
	static double [][]MMM = new double[Nn+1][Nn+1];
	
	static double [][]a1 = new double[N+1][5];
	static double [][]b1 = new double[N+1][5];
	static double [][]c1 = new double[N+1][5];
	static double [][]d1 = new double[N+1][5];
	
	static double [][]Km = new double[5][5];
	static double [][]SK = new double[Nn+1][Nn+1];
	
	static int Nm = (Nx+1)*(Nz+1)*2 + (Ny-1)*2*(Nx+Nz); //境界条件数
	static int []Ko = new int[Nm+1];
	
	static double []f = new double[Nn+1];
	static double []phi_0 = new double[Nn+1];
	static double []F = new double[Nn+1];
    
	static double step = 10000;
	static double dt = 10;
	
	static double []f_r = new double[Nn+1];
	
	static double []F_old = new double[Nn+1];
	static double []F_new = new double[Nn+1];
	
	static double esp_1 = 1/Math.pow(10,16);
	static double esp_2 = 0;
	
	static double [] r = new double[Nn + 1];
	static double [] p = new double[Nn + 1];
	static double [] ap = new double[Nn + 1];
	static double au,al,bu;
	static double fnorm,rnorm;
	static double alpha,beta,eps;
	static int round;
	
	static double [] ff = new double[Nn + 1];
	static double [] phi = new double[Nn + 1];
	
    public static void main(String[] arg){
    	tet_CG yuugenyousohou = new tet_CG();
    	yuugenyousohou.run();
    }
    
    public void run(){
    	init();
    	solve();
    	volume();
    	taisekizahyoukeisuu();
    	kakusangyouretu();
    	kousizahyou();
    	shokijyoukenn();
    	situryougyouretu();
    	kousijyoukenn();
    	uhen();
    	for (int n = 1;n <= step;n++){
    		keisann();
    		cgm();
    		l2();
    		System.out.println(n);
    		if(esp_2 <= esp_1){
        		n = (int)step;
        	}
    		irekae();
    	}
    	output();
    	output2();
    }
    
    public void init(){
    	System.out.println("1");
    	System.out.println("data");
    	System.out.println("step1");
    	System.out.println(((Nx + 1)*(Ny + 1)*(Nz + 1))+" "+(NN*6));
 //   	System.out.println();    	
    	
    	double a = 1;
    	
    	for (int i = 0;i <= Ny;i++){
    		for (int j = 0;j <= Nz;j++){
    			for (int k = 0;k <= Nx;k++){
    				ii = ii + 1;
/*    				y[ii] = Ly/2*Math.tanh(-Math.PI/2 + i*dy*Math.PI);
    				x[ii] = Lx/2*Math.tanh(-Math.PI/2 + k*dx*Math.PI);
    				z[ii] = Lz/2*Math.tanh(-Math.PI/2 + j*dz*Math.PI);
*/    				
    				y[ii] = Ly*(Math.exp(a*(-Math.PI/2 + i*dy*Math.PI))-Math.exp(-a*(-Math.PI/2 + i*dy*Math.PI)))
    				/(Math.exp(a*(-Math.PI/2 + i*dy*Math.PI))+Math.exp(-a*(-Math.PI/2 + i*dy*Math.PI)));
    				x[ii] = Lx*(Math.exp(a*(-Math.PI/2 + k*dx*Math.PI))-Math.exp(-a*(-Math.PI/2 + k*dx*Math.PI)))
    				/(Math.exp(a*(-Math.PI/2 + k*dx*Math.PI))+Math.exp(-a*(-Math.PI/2 + k*dx*Math.PI)));
    				z[ii] = Lz*(Math.exp(a*(-Math.PI/2 + j*dz*Math.PI))-Math.exp(-a*(-Math.PI/2 + j*dz*Math.PI)))
    				/(Math.exp(a*(-Math.PI/2 + j*dz*Math.PI))+Math.exp(-a*(-Math.PI/2 + j*dz*Math.PI)));
    			}
    		}
		}
    	
    	for(int i = 1;i <= Nn;i++){
    		y[i]=y[i]/y[Nn]/2;
    		x[i]=x[i]/x[Nn]/2;
    		z[i]=z[i]/z[Nn]/2;
    	}
    	for(int i = 1;i <= Nn;i++){ //座標の表示
   			System.out.println(i+" "+x[i]+" "+y[i]+" "+z[i]);
    	}
 //   	System.out.println();
    	
    }
    
    public void solve(){
    	int a = 0;
    	int m = 0;
    	int No = 0;
    	int b = 0;
    	int e = 0;
    	int f = 1;
    	int g = 0;
    	int count = 1;
    	for(int M = 1;M <= NN;M++){
  //  		count = count + 1;
    		if(No == Nx*Nz){
    			a = 0;
    			No = 1;
    			b = b+(Nx+1)*(Nz+1);
    	    	f = 1;
    	    	count = 1;
    		}else if(No == count*Nx){
    			a = a + (Nx+1);
    			e = -1;
    			f = f + 1;
    			No = No + 1;
    			count = count + 1;
        	}else{
    			No = No + 1;
    		}
    		
    		int d;
    		d = M - M/2*2;
    		g = f - f/2*2;
    		
 //   		if(count == 6 && No<=Nx){
    			e = e+1;
 //   			count = 0;
 //   		}
    			if(No == 1){
    				e = 0;
    			}
    		
    		if (g == 1){
    			if(d==0){//奇数
    				m = m + 1;
    				for(int n=1;n<=4;n++){//1
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//2
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}

    				m = m + 1;
    				for(int n=1;n<=4;n++){//3
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//4
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//5
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//6
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    			}else if(d==1){//偶数
    				m = m + 1;
    				for(int n=1;n<=4;n++){//1
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//2
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//3
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//4
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//5
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//6
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}	
    				}
    				
    			}else{
    				System.out.println("Error");
    			}
    		}else if(g == 0){
    			if(d==1){//偶数
    				m = m + 1;
    				for(int n=1;n<=4;n++){//1
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//2
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}

    				m = m + 1;
    				for(int n=1;n<=4;n++){//3
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//4
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//5
    	   	    		if(n==1){
    	   	    			ie[m][n]=1+(Nx+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    				m = m + 1;
    				for(int n=1;n<=4;n++){//6
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=1+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+b+e+a;
    	   	    		}	
    				}
    				
    			}else if(d==0){//奇数
    				m = m + 1;
    				for(int n=1;n<=4;n++){//1
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//2
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//3
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//4
    	   	    		if(n==1){
    	   	    			ie[m][n]=2+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//5
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=3+(Nx+1)+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}	
    				}
    				m = m + 1;
    				for(int n=1;n<=4;n++){//6
    	   	    		if(n==1){
    	   	    			ie[m][n]=3+(Nx+1)-1+b+e+a;
    	   	    		}else if(n==2){
    	   	    			ie[m][n]=2-1+b+e+a;
    	   	    		}else if(n==3){
    	   	    			ie[m][n]=3-1+b+e+a;
    	   	    		}else if(n==4){
    	   	    			ie[m][n]=3+(Nx+1)*(Nz+1)-1+b+e+a;
    	   	    		}	
    				}
    				
    			}else{
    				System.out.println("Error");
    			}}}
/*    	for(int M = 1;M <= N;M++){ //要素番号と要素接点番号の表示
    		System.out.println(M+" "+"1"+" "+"tet"+" "+ie[M][1]+" "+ie[M][2]+" "+ie[M][3]+" "+ie[M][4]);
    	}
*/    	
    	
    	
    }

    public void volume(){
    	Vn = 0.0;
    	for(int m = 1;m <= N;m++){
    		V[m]=0;
    		for(int n = 1;n <= 4;n++){
    			ia[n] = ie[m][n];
    			X1[n] = x[ia[n]];
    			Y1[n] = y[ia[n]];
    			Z1[n] = z[ia[n]];
    			
    		}
    		V[m] = ((X1[2]-X1[1])*((Y1[4]-Y1[1])*(Z1[3]-Z1[1])-(Y1[3]-Y1[1])*(Z1[4]-Z1[1]))
    				+(X1[4]-X1[1])*((Y1[3]-Y1[1])*(Z1[2]-Z1[1])-(Y1[2]-Y1[1])*(Z1[3]-Z1[1]))
    				+(X1[3]-X1[1])*((Y1[2]-Y1[1])*(Z1[4]-Z1[1])-(Y1[4]-Y1[1])*(Z1[2]-Z1[1])))/6;
    		Vn = Vn + V[m];    
    	}
    }
    
    public void taisekizahyoukeisuu(){
    	for(int m = 1;m <= N;m++){
    		for(int n = 1;n <= 4;n++){
    			ia[n] = ie[m][n];
    			X1[n] = x[ia[n]];
    			Y1[n] = y[ia[n]];
    			Z1[n] = z[ia[n]];
    		}
    		for(int i=1;i<=4;i++){
    		
    			if(i==1){
    				I = 1;J = 2;K = 3;L = 4;
    			}else if(i==2){
    				I = 2;J = 3;K = 4;L = 1;
    			}else if(i==3){
    				I = 3;J = 4;K = 1;L = 2;
    			}else if(i==4){
    				I = 4;J = 1;K = 2;L = 3;
    			}
    			
    			a1[m][i] = Math.pow(-1,i)*(Y1[K]*Z1[L] + Y1[L]*Z1[J] + Y1[J]*Z1[K]
    			                         - Y1[K]*Z1[J] - Y1[L]*Z1[K] - Y1[J]*Z1[L]);
    			b1[m][i] = Math.pow(-1,i)*(X1[J]*Z1[L] + X1[K]*Z1[J] + X1[L]*Z1[K]
    			                         - X1[L]*Z1[J] - X1[J]*Z1[K] - X1[K]*Z1[L]);
    			c1[m][i] = Math.pow(-1,i)*(X1[J]*Y1[K] + X1[K]*Y1[L] + X1[L]*Y1[J]
    			                         - X1[L]*Y1[K] - X1[J]*Y1[L] - X1[K]*Y1[J]);
    			d1[m][i] = Math.pow(-1,i+1)*(X1[J]*Y1[K]*Z1[L] + X1[K]*Y1[L]*Z1[J] + X1[L]*Y1[J]*Z1[K]
    			                           - X1[L]*Y1[K]*Z1[J] - X1[J]*Y1[L]*Z1[K] - X1[K]*Y1[J]*Z1[L]);
    		}
    	
    	}
    }
    
    public void kakusangyouretu(){
    	for(int i = 1;i <= Nn;i++){
    		for(int j = 1;j <= Nn;j++){
    			SK[i][j] = 0;
    		}
    	}
    	
    	I = 1;J = 2;K = 3;L = 4;
    	for(int m = 1;m <= N;m++){
    		Km[1][1] = (Math.pow(a1[m][I],2) + Math.pow(b1[m][I],2) + Math.pow(c1[m][I],2))/(6.0*V[m]);
    		Km[1][2] = (a1[m][I]*a1[m][J] + b1[m][I]*b1[m][J] + c1[m][I]*c1[m][J])/(6.0*V[m]);
    		Km[1][3] = (a1[m][I]*a1[m][K] + b1[m][I]*b1[m][K] + c1[m][I]*c1[m][K])/(6.0*V[m]);
    		Km[1][4] = (a1[m][I]*a1[m][L] + b1[m][I]*b1[m][L] + c1[m][I]*c1[m][L])/(6.0*V[m]);
    		
    		Km[2][1] = (a1[m][J]*a1[m][I] + b1[m][J]*b1[m][I] + c1[m][J]*c1[m][I])/(6.0*V[m]);
    		Km[2][2] = (Math.pow(a1[m][J],2) + Math.pow(b1[m][J],2) + Math.pow(c1[m][J],2))/(6.0*V[m]);
    		Km[2][3] = (a1[m][J]*a1[m][K] + b1[m][J]*b1[m][K] + c1[m][J]*c1[m][K])/(6.0*V[m]);
    		Km[2][4] = (a1[m][J]*a1[m][L] + b1[m][J]*b1[m][L] + c1[m][J]*c1[m][L])/(6.0*V[m]);
    		
    		Km[3][1] = (a1[m][K]*a1[m][I] + b1[m][K]*b1[m][I] + c1[m][K]*c1[m][I])/(6.0*V[m]);
    		Km[3][2] = (a1[m][K]*a1[m][J] + b1[m][K]*b1[m][J] + c1[m][K]*c1[m][J])/(6.0*V[m]);
    		Km[3][3] = (Math.pow(a1[m][K],2) + Math.pow(b1[m][K],2) + Math.pow(c1[m][K],2))/(6.0*V[m]);
    		Km[3][4] = (a1[m][K]*a1[m][L] + b1[m][K]*b1[m][L] + c1[m][K]*c1[m][L])/(6.0*V[m]);
    		
    		Km[4][1] = (a1[m][L]*a1[m][I] + b1[m][L]*b1[m][I] + c1[m][L]*c1[m][I])/(6.0*V[m]);
    		Km[4][2] = (a1[m][L]*a1[m][J] + b1[m][L]*b1[m][J] + c1[m][L]*c1[m][J])/(6.0*V[m]);
    		Km[4][3] = (a1[m][L]*a1[m][K] + b1[m][L]*b1[m][K] + c1[m][L]*c1[m][K])/(6.0*V[m]);
    		Km[4][4] = (Math.pow(a1[m][L],2) + Math.pow(b1[m][L],2) + Math.pow(c1[m][L],2))/(6.0*V[m]);
    		    		
    		for(int j = 1;j <= 4;j++){
    			int iii = ie[m][j];
    			for(int k = 1;k <= 4;k++){
    				int jjj = ie [m][k];
    				SK[iii][jjj] = SK[iii][jjj] + Km[j][k];
    			}
    		}
    	}
    }
    
    public void kousizahyou(){
    	int j=0;
    	
    	for(int i = 1;i <= Nn;i++){
    		if(y[i] == y[1]  || y[i] == y[Nn]){
    			j=j+1;
    			Ko[j] = i;
    		}else if(x[i] == x[1] || x[i] == x[Nn]){
    			j=j+1;
    			Ko[j] = i;
    		}else if(z[i] == z[1] || z[i] == z[Nn]){
    			j=j+1;
    			Ko[j] = i;
    		}else{}
    		
    	}
    }
    
    public void shokijyoukenn(){
    	for(int i = 1;i <= Nn;i++){
    		if (i <= (Nx+1)*(Nz+1)){
    			phi_0[i] = 10.0;
    		}else{
    			phi_0[i] = 0.0;
    		}
//    		System.out.println("phi_0 = " + phi_0[i]);
    		F_old[i] = phi_0[i];
    		F[i] = phi_0[i];
    	}
    	 
    }
    
    public void situryougyouretu(){
    	for(int i = 1;i <= Nn;i++){
    		for(int j = 1;j <= Nn;j++){
    			SM[i][j] = 0;
    		}
    	}
    	I = 1;J = 2;K = 3;
    	for(int m = 1;m <= N;m++){
    		M[1][1] = V[m]/10.0;
    		M[1][2] = V[m]/20.0;
    		M[1][3] = V[m]/20.0;
    		M[1][4] = V[m]/20.0;
    		
    		M[2][1] = V[m]/20.0;
    		M[2][2] = V[m]/10.0;
    		M[2][3] = V[m]/20.0;
    		M[2][4] = V[m]/20.0;
    		
    		M[3][1] = V[m]/20.0;
    		M[3][2] = V[m]/20.0;
    		M[3][3] = V[m]/10.0;
    		M[3][4] = V[m]/20.0;
    		
    		M[4][1] = V[m]/20.0;
    		M[4][2] = V[m]/20.0;
    		M[4][3] = V[m]/20.0;
    		M[4][4] = V[m]/10.0;
    		
    		for(int j = 1;j <= 4;j++){
    			int ii = ie[m][j];
    			for(int k = 1;k <= 4;k++){
    				int jj = ie [m][k];
    				MM[ii][jj] = MM[ii][jj] + M[j][k];
    			}
    		}
    	}
    	
    	for(int i = 1;i <=Nn;i++){
    		for(int j = 1;j <= Nn;j++){
    			SM[i][j] = MM[i][j] + dt*SK[i][j];
    	    }
    	}
 
    	for(int i = 1;i <= Nn;i++){
    		for(int j = 1;j <= Nn;j++){
    			MMM[i][j] = MM[i][j];
    			SMM[i][j] = SM[i][j];
    		}
    	}
    }
    
    public void kousijyoukenn(){
    	for(int i = 1;i <= Nm;i++){
    		int j = Ko[i];
    		if (j <= (Nx+1)*(Nz+1)){
    			f[j] = 10.0;
    		}else{
    			f[j] = 0.0;
    		}
    	}
    }
    
    public void uhen(){
    	for(int i = 1;i <= Nn;i++){
    		f_r[i] = 0;
    	}
    	for(int i = 1;i <= Nm;i++){
    		int ii = Ko[i];
    		for (int j = 1;j <= Nn;j++){
    			f_r[j] = f_r[j] - SM[ii][j]*f[i];
    			SM[ii][j] = 0.0;
    			SM[j][ii] = 0.0;
    		}
    		SM[ii][ii] = 1.0;
    	}
    	for(int i = 1;i <= Nm;i++){
    		ii = (int)Ko[i];
    		f_r[ii] = f[i]; 
    	}
    }
    
    public void keisann(){
    	for(int i = 1;i <= Nn;i++){
    		ff[i] = f_r[i];
    	}
    	for(int m = 1;m <= Nn;m++){
    		for(int i = 1;i <= Nn;i++){
    			MMM[m][i] = MM[m][i];
    			ff[m] = ff[m] + MMM[m][i]*phi_0[i];
    		}
    	}
    	
    	for(int i = 1;i <= Nm;i++){
    		ii = (int)Ko[i];
    		ff[ii] = f_r[ii];
    	}
    }
    
    public void l2(){
    	double A = 0;
    	double B = 0;
    	esp_2 = 0;
    	for(int i = 1;i <= Nn;i++){
    		A = A + ((F_new[i] - F_old[i])*(F_new[i] - F_old[i]));
    		B = B + (F_old[i]*F_old[i]);
    	}
    	esp_2 = A/B;
    }
    
    public void irekae(){
    	for(int i = 1;i <= Nn;i++){
    		F_old[i] = phi_0[i];
    		phi_0[i] = F[i];
//    		System.out.println("F[" + i + "] = " + F[i]);
//    		System.out.println("phi[" + i + "] = " + phi_0[i]);
    	}
    	for(int i = 1;i <= Nn;i++){
    		for(int j = 1;j <= Nn;j++){
    			SM[i][j] = SMM[i][j];
    		}
    	}
    	
    	for(int i = 1;i <= Nm;i++){
    		int ii = Ko[i];
    		for (int j = 1;j <= Nn;j++){
    			SM[ii][j] = 0.0;
    			SM[j][ii] = 0.0;
    		}
    		SM[ii][ii] = 1.0;
    	}
    }
    
    public void cgm(){
    	beta = 0.0;
    	fnorm = 0.0;
    	
    	prospx(phi_0,ap);
    	
    	for(int i = 1;i <= Nn;i++){
    		r[i] = ff[i] - ap[i];
			p[i] = r[i];
			fnorm = fnorm + ff[i]*ff[i]; 
    	}
    	
    	do{
			prospx(p,ap);
			
			au = 0.0;
			al = 0.0;
			
			for(int i = 1;i <= Nn;i++){
				au = au + r[i]*r[i];
				al = al + p[i]*ap[i];
			}
			
			alpha = au/al;
			bu = 0.0;
			
			for(int i = 1;i <= Nn;i++){
				F[i] = F[i] + alpha*p[i];
				r[i] = r[i] - alpha*ap[i];
				bu = bu + r[i]*r[i];
				F_new[i] = F[i];
			}
			
			beta = bu/au;
			rnorm = bu;
			eps = Math.sqrt(rnorm/fnorm);
			
			for(int i = 1;i <= Nn;i++){
				p[i] = r[i] + beta*p[i];
			}
			
/*			for(int i = 1;i <= Nn;i++){
				System.out.println("F[" + i + "]=" + F[i]);
			}
*/			round++;
		
    	}while(eps > Math.pow(10,-3));
		

    }
    
    public double[] prospx(double[] x,double[] res){
		
    	for(int i = 1;i <= Nn;i++){
			double tmp = 0.0;
			for(int j = 1;j <= Nn;j++){
				tmp = tmp + SM[i][j]*x[j];
			}
			res[i] = tmp;
		}
    	
		return res;
	}
    
    
    public void output(){
    	for(int m = 1;m <= N;m++){ //要素番号と要素接点番号の表示
    		System.out.println(m+" "+"1"+" "+"tet"+" "+ie[m][1]+" "+ie[m][2]+" "+ie[m][3]+" "+ie[m][4]);
    	}
		System.out.println("1 0");
		System.out.println("1 1");
		System.out.println("temperature,");
    	    	
    	    }
    
    public void output2(){
   	    	
    	for(int i = 1;i <= Nn;i++){ 
    		System.out.println(i+" "+F[i]);
    	}
    }
}