// elem1.cpp: implementation of the elem class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "四接點兩單元的平面靜力問題.h"
#include "elem1.h"
#include "math.h"
#include<cstdlib>
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

double Elem::Mu = 0;     //初始化靜態數據成員 
double Elem::E = 1;      // 
double Elem::t = 1;      //
int Elem::ndTotNum = 0;
int Elem::elTotNum = 0;

                
 void Elem::setStatics(double vMu, double vE, double vt, int flag, int ntn, int ctn)
{                                                             
	t = vt;	
	if(flag)
	{
		Mu = vMu/(1 - vMu);
		E = vE/(1 - vMu * vMu);		
		return;
	}
	Mu = vMu;
	E = vE;
	ndTotNum = ntn;
	elTotNum = ctn;	
}

Elem::Elem()
{
	ndp[0] = 0;
	ndp[1] = 0;
	ndp[2] = 0;
	elNum = -1;
    sgmX = 0;     //單元應力 
    sgmY = 0;
    tauXY = 0;
    double	 SgmMax = 0;   //主應力 
    double	 SgmMin = 0;   //
    double	 Theta = 0;    //應力主角 
}

void Elem::setElem(int n, Node* nd1, Node* nd2, Node* nd3)
{ 
	if(elNum != -1)
	{
	}	
	if(nd1 == nd2 || nd1 == nd3 || nd2 == nd3)
	{
	
	}
     
 	ndp[0] = nd1;
	ndp[1] = nd2;
	ndp[2] = nd3;
	elNum = n;
	getAbc();
}

void Elem::getAbc()      //計算變量數組a[],b[],c[]及面積elArea
{	
	for(int i = 0; i < 3; i++)
	{
		int j = (i + 1)%3;
		int m = (i + 2)%3;
		a[i] = (ndp[j]->x) * (ndp[m]->y) - (ndp[m]->x) * (ndp[j]->y);
		b[i] = (ndp[j]->y) - (ndp[m]->y);
		c[i] = (ndp[m]->x) - (ndp[j]->x);
	}
 	area = (a[0] + a[1] + a[2])/2;
 	
  
}

void Elem::setMtxk(Mtx<double>& mk)
{
	double md = E*t/(4*(1 - Mu*Mu)*area);
 	for(int i = 0; i < 3; i++)
  		for(int j = 0; j < 3; j++)
  		{
  			int r = 2 *(ndp[i]->ndNum);
  			int s = 2 *(ndp[j]->ndNum);
  						
  			mk(r, s) += md*(b[i]* b[j] + (1 - Mu)*c[i]*c[j]/2);  			
  			mk(r, s+1) += md*(Mu*b[i]*c[j] + (1 - Mu)*c[i]*b[j]/2);  			
  			mk(r+1, s) += md*(Mu*c[i]*b[j] + (1 - Mu)*b[i]*c[j]/2);  			
  			mk(r+1, s+1) += md*(c[i]*c[j] + (1 - Mu)*b[i]*b[j]/2);
		}
}

void Elem::setStress()        //計算單元應力 
{	
   
	for(int i = 0; i < 3; i++)
	{
		sgmX += (b[i]*(ndp[i]->u) + Mu*c[i]*(ndp[i]->v));
  		sgmY += (Mu*b[i]*(ndp[i]->u) + c[i]*(ndp[i]->v));
    	tauXY +=(1 - Mu)*(c[i]*(ndp[i]->u) + b[i]*(ndp[i]->v))/2;	
	}
	sgmX *= E/(2*area*(1 - Mu*Mu));
	sgmY *= E/(2*area*(1 - Mu*Mu));
	tauXY *= E/(2*area*(1 - Mu*Mu));
	
	double tmp = sqrt((sgmX - sgmY)*(sgmX - sgmY)/4 + tauXY*tauXY);
 	sgmMax = (sgmX + sgmY)/2 + tmp;
	sgmMin = (sgmX + sgmY)/2 - tmp;
	theta = 90*atan(tauXY/(sgmY - sgmMin))/asin(1.0);	
}

void Elem::output()
{
  CString ss;
  ss.Format("%f",sgmX);
  AfxMessageBox(ss);
}