elem.h

一个简单的三角形单元有限元计算程序,是标准的C++源程序

//文件名：elem.h
//类型：头文件
//内容：节点类及单元类的定义 

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

#ifndef _ELEM_H
#define _ELEM_H

#include "matrix.h"

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


class Node;
class Elem;

extern inline void error(char* v);       // 辅助函数
/*
inline void error(char* v)       // 辅助函数
{
	cout << v << "...program exited\n";
	system("pause");
	exit(1);            // 包括<stdlib.h>
}
*/
////////////////////////  INCLUDE ENDED  ///////////////////////////////////////

//////////////////////////  节点类定义  ////////////////////////////////////////

class Node
{
	double x;         //节点坐标 
   	double y;         // 
   	int ndNum;        //节点编号 
   	
 	double u;         //节点位移 
 	double v;         //
 	
 	double loadX;     //节点等效荷载 
 	double loadY;     //
        
public:
	Node();
	~Node(){ };
	
	void setPoint(int, double, double);
	void setUx(double uu) { u = uu; };
	void setVy(double vv) { v = vv; };
	void setLoad(double lx, double ly);
	void setLoadX(double lx) { loadX = lx; }
	void setLoadY(double ly) { loadY = ly; }
	double getLoadX() const { return loadX; }
	double getLoadY() const { return loadY; }
	int getNum() const { return ndNum; }
	
	friend ostream& operator<<(ostream&, const Node&);
	friend class Elem;//定义单元类Elem为其友元类	
};

Node::Node()
{
	x = 0;
	y = 0;
	loadX = 0;
	loadY = 0;
	u = 0;
	v = 0;
	ndNum = -1;
}

void Node::setPoint(int nn, double xx, double yy)
{
	if(ndNum != -1)
	{
		cout<< "multi-definded Node: " << nn << endl;
		error(" ");
	}
 	ndNum = nn;
 	x = xx;
	y = yy;
}

void Node::setLoad(double lx, double ly)
{
	loadX += lx;    // 注意是“ += ” 
	loadY += ly;
}

ostream& operator<<(ostream& s, const Node& nd)
{
	s << "Node number: " << nd.ndNum 
 	  << "  (" << nd.x << ", " << nd.y << ")  ";
	s << "u = " << nd.u << "  v = " << nd.v <<endl;
	//s << "loadX = " << nd.loadX << "  loadY = " << nd.loadY <<endl;
	
 	return s;
}
////////////////////////////////////////////////////////////////////////////////

///////////////////////////////  单元类定义  ///////////////////////////////////

class Elem
{
	Node* ndp[3];         //节点指针 
	
	int elNum;
	
	double sgmX;   //单元应力 
	double sgmY;
	double tauXY;
	double sgmMax;   //主应力 
	double sgmMin;   //
	double theta;    //应力主角 
	
	double area;      //单元面积
 	double a[3];        //计算中间变量 
  	double b[3];
  	double c[3];
  	
  	static double Mu;   //泊松比 
  	static double E;    //弹性模量 
  	static double t;    //厚度 
  	static int ndTotNum;//总节点数 
  	static int elTotNum;//总单元数   	
  	
  	void getAbc();      //计算变量数组a[],b[],c[]及面积elArea  
  	
public:
	Elem();
	~Elem() { };
    
	void setElem(int, Node*, Node*, Node*);                 //设置成员变量值 
	void setMtxk(Mtx<double>&);     //生成整体刚度矩阵
	void setStress();                  //计算单元应力 
	
	static void setStatics(double, double, double, int, int, int);
 	static int getNdTN() { return ndTotNum; }; 
  	static int getElTN() { return elTotNum; };	
	
	friend ostream& operator<<(ostream&, const Elem&);
};


double Elem::Mu = 0;     //初始化静态数据成员 
double Elem::E = 1;      // 
double Elem::t = 1;      //
int Elem::ndTotNum = 0;
int Elem::elTotNum = 0;

                
inline void Elem::setStatics(double vMu, double vE, double vt, 
                             int flag, int ntn, int ctn)
                             //flag为零则为平面应力问题，否则为平面应变问题 
{                                                             
	if(vMu < 0 || vMu >=1) error("wrong value of Mu");
	if(vE <= 0) error("wrong value of E");
	if(vt <= 0) error("wrong value of t");
	if(flag != 0 && flag != 1) error("wrong value of flag");
	if(ntn <= 0) error("wrong Nodes number");
	if(ctn <= 0) error("wrong Elems number");
	
	if(flag) 
 		cout << "\n平面应变问题" << endl; 
	else
		cout << "\n平面应力问题" << endl; 
		
	cout << "Mu = " << vMu << "  E = " << vE << "  t = " << vt << endl;
	cout << "Nodes number: " << ntn << "  Elems number: " << ctn << endl;
	
 	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;
	double SgmX = 0;     //单元应力 
	double SgmY = 0;
	double 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)
	{
		cout<< "multi-definded Elem: " << n << endl;
		error(" ");
	}	
	if(nd1 == nd2 || nd1 == nd3 || nd2 == nd3)
	{
		cout << "same nodes in Elem " << elNum;
		error(" ");
	}
     
 	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;
 	
  	if(area < 0)
   	{
   		cout << "Bad order of these nodes in Elem " << elNum << " : " 
             << ndp[0]-> ndNum << "  " << ndp[1]-> ndNum << "  " 
             << ndp[2]-> ndNum << endl;             
		error("");
	}	
}

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);	
}

ostream& operator<<(ostream& s, const Elem& el)
{	
  	s << "Elem number: " << el.elNum;
 	s << "    included nodes: " << el.ndp[0]->getNum() << "  "
	               << el.ndp[1]->getNum() << "  "
	               << el.ndp[2]->getNum();
	s << "\tarea = " << el.area << endl;	
	s << "\tsgmX = " << el.sgmX 
 	  << "\tsgmY = " << el.sgmY 
	  << "\ttauXY = " << el.tauXY <<endl;
	s << "\tsgmMax = " << el.sgmMax 
	  << "\tsgmMin = " << el.sgmMin 
	  << "\ttheta = " << el.theta <<endl;
	
 	return s;
}

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

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

#endif 

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