////////////  常應變三角形單元解彈性力學平面問題的有限元程序  //////////////////
//                                                                            //
//      File name:	ttfem.cpp                                             // 
//      Compiler:	Turbo C 2.0                                           //
//      Author:		LiuH                                                  // 
//      Date:		2004-5-4                                              //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////  INCLUDE  /////////////////////////////////////

#include "elem.h"
#include "matrix.h"

#include<iostream>
#include<fstream>
#include<cstdlib>
using namespace std;

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

extern inline void error(char* v);       // 輔助函數

void readData(Node*&, Elem*&, Mtx<double>*&);     //讀取數據并生成整體剛度矩陣 
void solveEq(Node*&, Mtx<double>*&);              //解方程子程序
void stress(Elem*&);                              //求應力子程序
void output(Node*&, Elem*&);                      //輸出位移及應力子程序  

////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////

int main()
{
	Node* nptr = 0; //節點指針 
	Elem* eptr = 0; //單元指針 
	Mtx<double>* kptr = 0;               //整體剛度矩陣指針，初始值為0 	
	
	readData(nptr, eptr, kptr);//讀取節點數、單元數及常量Mu, E, t 
	solveEq(nptr, kptr);
	stress(eptr);
 	output(nptr, eptr); 
 	
  	delete[] nptr;
 	delete[] eptr;
 	delete kptr;
 	
 	cout << "\n  ...program ended.\n\n\n\n";
  	system("pause");
 	return 0;
}

///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////

void readData(Node*& ndp, Elem*& elp, Mtx<double>*& kp)
{
	char fname[20];
  	cout << "Please input the name of data file: ";
  	cin >> fname;
  	  	
  	ifstream infile;
   	infile.open(fname, ios_base::in);
   	if(!infile)
   	{	
   		cout << fname << ": Can not open this file" << endl;
   		error("  ");
 	}
 	
   	double dtmp1, dtmp2, dtmp3;     //定義輔助變量 
	int itmp1, itmp2, itmp3, itmp4;  	
	infile >> dtmp1 >> dtmp2 >> dtmp3;
 	infile >> itmp1 >> itmp2 >> itmp3;
	Elem::setStatics(dtmp1, dtmp2, dtmp3, itmp1, itmp2, itmp3); 

	const int nNum = Elem::getNdTN();
 	const int eNum = Elem::getElTN();
	
	ndp = new Node[nNum];                  //定義節點數組 
	elp = new Elem[eNum];                  //定義單元數組 
	kp = new Mtx<double>(2*nNum, 2*nNum);  //定義整體剛度矩陣 
  	
   	int count;   	
   	infile >> count;        //檢查結束標志-1
	if(count != -1) error("wrong constants number");
        
	// 讀取節點數據 
 	count = 0;            
	while(1)
	{
   		infile >> itmp1;   		
   		if(itmp1 == -1) break;
   		
		if(itmp1 < 0 || itmp1 >= nNum)        ////檢查讀入數據的合法性 
		{
			cout << itmp1;
			error(": wrong node number");
		}   		
		infile >> dtmp2 >> dtmp3;
     	
		ndp[itmp1].setPoint(itmp1, dtmp2, dtmp3);
		count++;  		
	}
	if(count != nNum) error("wrong nodes number");
   
	//讀取單元數據 
	count = 0;    
	while(1)
	{
		infile >> itmp1;
		if(itmp1 == -1) break;
    	
		if(itmp1 < 0 || itmp1 >= eNum)      //檢查讀入數據的合法性 
  		{
  			cout << itmp1;
  			error(": wrong Elem number");
   		}   
		infile >> itmp2 >> itmp3 >> itmp4;
		if(itmp2 < 0 || itmp2 >= nNum || itmp2 < 0 || itmp2 >= nNum
		   || itmp2 < 0 || itmp2 >= nNum)
  		{
			cout << "wrong Node number for Elem " << itmp1;
  			error(" ");
   		}   
    	
		elp[itmp1].setElem(itmp1, &ndp[itmp2], &ndp[itmp3], &ndp[itmp4]);
		count++;
   	}
   	if(count != eNum) error("wrong elems number");
  	
   	//生成整體剛度矩陣 
	for(int i = 0; i < eNum; i ++)
   		elp[i].setMtxk(*kp);
   	
   	//讀入節點等效荷載
   	while(1)
   	{
   		infile >> itmp1;
		if(itmp1 == -1) break;
     
		if(itmp1 < 0 || itmp1 >= nNum)      //檢查讀入數據的合法性 
		{
			cout << itmp1;
			error(": wrong Node number");
		}
        
		infile >> dtmp2 >> dtmp3;
		ndp[itmp1].setLoad(dtmp2, dtmp3);        	
	}
         
	//讀取支撐條件，并對整體剛度矩陣做相應的改動 
	const double BigNum = 1000000000;  //設置一個大系數 
	
	while(1)                           //讀取x方向支撐條件 
   	{
   		infile >> itmp1;
		if(itmp1 == -1) break;
     
		if(itmp1 < 0 || itmp1 >= nNum)      //檢查讀入數據的合法性 
		{
			cout << itmp1;
			error(": wrong Node number");
		}
        
		infile >> dtmp2; 
		if(dtmp2 == 0)                  //水平位移為0 
  		{
  			ndp[itmp1].setLoadX(0);
			for(int i = 0; i < 2*nNum; i++)
			{
				if(i == 2*itmp1)
					(*kp)(i, 2*itmp1) = 1;
   		 		else
				{
 					(*kp)(i, 2*itmp1) = 0;
  					(*kp)(2*itmp1, i) = 0;
				}
			}
		}
		else                            //水平位移不為0  
		{
			(*kp)(2*itmp1, 2*itmp1) *= BigNum;
   			ndp[itmp1].setLoadX((*kp)(2*itmp1, 2*itmp1)*dtmp2); 
		}	
   	}
   	
   	while(1)                                //讀取y方向支撐條件 
   	{
   		infile >> itmp1;
		if(itmp1 == -1) break;
     
		if(itmp1 < 0 || itmp1 >= nNum)      //檢查讀入數據的合法性 
		{
			cout << itmp1;
			error(": wrong Node number");
		}
        
		infile >> dtmp3;
		        	
		if(dtmp3 == 0)                  //垂直位移為0 
		{
			ndp[itmp1].setLoadY(0);
			for(int i = 0; i < 2*nNum; i++)
			{
				if(i == 2*itmp1 + 1)
					(*kp)(i, 2*itmp1 +1) = 1;
				else
   				{
   					(*kp)(i, 2*itmp1 +1) = 0;
   					(*kp)(2*itmp1 +1, i) = 0;
   				}
			}
		}
		else                            //垂直位移不為0 
		{
			(*kp)(2*itmp1 +1, 2*itmp1 +1) *= BigNum;
			ndp[itmp1].setLoadY((*kp)(2*itmp1 +1, 2*itmp1 +1)*dtmp3);
		}
   	}
   	   	
	infile.close();
}

void solveEq(Node*& ndp, Mtx<double>*& aa)
{
	const int nNum = Elem::getNdTN();
	const int eNum = Elem::getElTN();
	const int tnNum = 2*nNum;
 		
	double* bb = new double[tnNum];
	double* xx = new double[tnNum];
 	for(int i = 0; i < nNum; i++)
  	{
  		bb[2*i] = ndp[i].getLoadX();
  		bb[2*i +1] = ndp[i].getLoadY();
   	} 
   	
   	//for(int i = 0; i < tnNum; i++)   //輔助調試 
   	//	cout << bb[i] << "\t";
 	//cout << endl;
     	
	for(int k = 0; k < tnNum -1; k++)                  //高斯消元法消元 
	{
     		double aet = abs((*aa)(k, k));             //挑選主元素 
     		int pc = k;                 
     		for(int i = k + 1; i < tnNum; i++)
     		{
     			if(abs((*aa)(i, k)) > aet)
        		{
        			aet = abs((*aa)(i, k));
        			pc = i;
          		}	
     		}
     		if(!aet)
     			error("pivot is zero in solveEq()");
     			
		if(pc != k)                                //交換行 
		{
			double tmp;
			for(int i = 0; i < tnNum; i++)
			{
				tmp = (*aa)(k, i);
				(*aa)(k, i) = (*aa)(pc, i);
				(*aa)(pc, i) = tmp;
			}
			tmp = bb[k];
			bb[k] = bb[pc];
			bb[pc] = tmp;
		}
		
		for(int i = k + 1; i < tnNum; i++)
		{
			(*aa)(i, k) /= (*aa)(k, k);
			bb[i] -= (*aa)(i, k)*bb[k];
  			for(int j = k +1; j < tnNum; j++)
     				(*aa)(i, j) -= (*aa)(i, k)*(*aa)(k, j);
          	}	
	}
	
	for(int i = tnNum - 1; i >= 0; i--)                //回代求解 
	{
		double tmp = 0;
		for(int j = i + 1; j < tnNum; j++)
			bb[i] -= (*aa)(i, j)*xx[j];
		xx[i] = bb[i]/(*aa)(i, i);
	}
	
	for(int i = 0; i < nNum; i++)
  	{
  		 ndp[i].setUx(xx[2*i]);
  		 ndp[i].setVy(xx[2*i +1]);
   	} 
   	
   	//for(int i = 0; i < tnNum; i++)   //輔助調試 
   	//	cout << xx[i] << "\t";
 	//cout << endl;
}

void stress(Elem*& clp)                              //求應力子程序
{
	const int eNum = Elem::getElTN();
	for(int i = 0; i < eNum; i++)
		clp[i].setStress();
}

void output(Node*& ndp, Elem*& elp)                      //輸出位移及應力子程序  
{
	char fname[20];
  	cout << "\nPlease input the name of output data file: ";
  	cin >> fname;
  	  	
  	ofstream outfile;
   	outfile.open(fname, ios_base::out);
   	if(!outfile)
   	{	
   		cout << fname << ": Can not open this file" << endl;
   		error("  ");
 	}
 	
 	const int nNum = Elem::getNdTN();
 	const int eNum = Elem::getElTN();
 	
 	cout << endl;
 	outfile << endl;
 	for(int i = 0; i < nNum; i++)
 	{
 		outfile << ndp[i];
 		cout << ndp[i];
 	}
 	
 	cout << endl; 
  	outfile << endl;	
 	for(int i = 0; i < eNum; i++)
 	{
 		outfile << elp[i];
 		cout << elp[i];
 	}
 	
 	cout << "\n  ...output successfully.\n";
 	outfile.close();
}

///////////////////////////  END  /////////////////////////////////////////////