staticFEM.java
package scu.edu.fem;
import java.math.BigInteger;
class ETNode{ //单元结点结构体
double x=0,y=0; //单元结点坐标
int number=0; //单元结点在总体区域划分中的编号
}
class ElementTriangle{ //三角形单元结构体
ETNode nd[]=new ETNode[]{new ETNode(),new ETNode(),new ETNode()}; //存储相对应的总体结点号
double a[]=new double[3];
double b[]=new double[3];
double c[]=new double[3]; //基函数的系数
double A; //单元面积
// double Aij[3][3]; //单元矩阵
//jimmy 改成6*6的二维
double Aij[][]= new double[6][6]; //单元矩阵(K)
double M[][]=new double[6][6];; //单元矩阵(M)
double C[][]=new double[6][6];; //单元矩阵(C)
double fi[]=new double[3];; //单元上的积分值
}
public class StaticFEM {
// int i; //单元的循环指标
// int j;
// int k;
//int id;
static double dBig=1000000000000000000000000000000f;
//jimmy change for test
//const int nx=1,ny=1; //x,y方向划分网格数,三角形单元个数=nx*ny*2
//jimmy change for test 必须为奇数
final static int nx=1,ny=1; //x,y方向划分网格数,三角形单元个数=nx*ny*2
//const int nx=21,ny=21; //x,y方向划分网格数,三角形单元个数=nx*ny*2
final static double Lx=10; //矩形区域的两边长
final static double Ly=10;
static double D; //单元矩阵行列式值
static double E=1000000f; //模量
static double t=1.0; //厚度
static double u=0.3; //泊松比
static double rou=1.0; //密度
static double Y=0.1; //阻尼系数
static double T=1;
static double dataT=0.001; //时间步长
static double gama=0.5; //参数gama>=0.5
static double bata=1;//0.25*(0.5+gama)*(0.5+gama); //参数bata
static double c0=1/(bata*dataT*dataT); //参数c0
static double c1=gama/(bata*dataT); //参数c1
static double c2=1/(bata*dataT); //参数c2
static double c3=1/(bata*2)-1; //参数c3
static double c4=gama/bata-1; //参数c4
static double c5=dataT/2*(c4-1); //参数c5
static double c6=dataT*(1-gama); //参数c6
static double c7=gama*dataT; //参数c7
static double pMatrix[]; //总体K矩阵指针
static double totalM[]; //总体质量M矩阵指针
static double totalC[]; //总体阻尼C矩阵指针
static double dalta[]; //位移
static double dalta1[]; //速度
static double dalta2[]; //加速度
static double pMf[]; //f向量指针
static double tmp[];
static ElementTriangle[] pE; //单元三角形结构体数组指针
double ai,bi,ci; //基函数的系数
//-----------------------------------------------------
double fn(int n , double x[])//被积函数,高斯积分函数的参数
{
return(ai+bi*x[0]+ci*x[1]);
}
//---------- 全选主元高斯消去法 ------------------------------
// a 体积为n*n的双精度实型二维数组,方程组系数矩阵,返回时将被破坏
// b 长度为n的双精度实型一维数组,方程组右端的常数向量,返回方程组的解向量
// n 整型变量,方程组的阶数
//--------------------------------------------------------------
static int Gauss(double a[],double b[],int n)
{
int []js;
int l,k,i,j,is=0,p,q;
double d,t;
js=new int[n];
l=1;
for(k=0;k<=n-2;k++)
{
d=0.0;
for(i=k;i<=n-1;i++)
for(j=k;j<=n-1;j++)
{
t=Math.abs(a[i*n+j]);
if(t>d) { d=t; js[k]=j; is=i;}
}
if(d+1.0==1.0) l=0;
else
{
if(js[k]!=k)
for(i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if(is!=k)
{
for(j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if(l==0)
{
System.out.println("Gauss funtion failed 1...\n");
return(0);
}
d=a[k*n+k];
for(j=k+1;j<=n-1;j++)
{
p=k*n+j; a[p]=a[p]/d;
}
b[k]=b[k]/d;
for(i=k+1;i<=n-1;i++)
{
for(j=k+1;j<=n-1;j++)
{
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if(Math.abs(d)+1.0==1.0)
{
//free(js);
System.out.println("Gauss funtion failed 2...\n");
return(0);
}
b[n-1]=b[n-1]/d;
for(i=n-2;i>=0;i--)
{
t=0.0;
for(j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for(k=n-1;k>=0;k--)
if(js[k]!=k)
{
t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
}
return 1;
}
public static void main(String[] args) {
//-------------- 全局变量 ---------------------------
StaticFEM staticFEM=new StaticFEM();
final int iNode=(nx+1)*(nx+1); //结点个数
//int main(int argc, char* argv[])
//{
//为总体矩阵,三角形单元数组,f函数向量分配存储内存
pMatrix=new double[2*iNode*2*iNode];
totalM=new double[2*iNode*2*iNode];
totalC=new double[2*iNode*2*iNode];
pE=new ElementTriangle[nx*ny*2];
for(int i=0;i<nx*ny*2;i++)
pE[i]=new ElementTriangle();
//jimmy 修改为2*iNode,因为有x,y两个方向载荷或约束
pMf=new double[2*iNode];
tmp=new double[2*iNode];
dalta=new double[2*iNode];
dalta1=new double[2*iNode];
dalta2=new double[2*iNode];
//初始化值为0,因为下面要累加总体矩阵
for(int i=0;i<2*iNode*2*iNode;i++)
totalC[i]=totalM[i]=pMatrix[i]=0.0;
for(int i=0;i<2*iNode;i++)
dalta2[i]=dalta[i]=dalta1[i]=pMf[i]=0.0;
try{
//------ 计算得到网格的信息 -----------
for(int j=0;j<nx;j++)
for(int i=0;i<ny;i++)
{
//for the first triangle in the rectangle
pE[i*2+j*ny*2].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2].nd[0].number=i+j*(nx+1); //NO.0
pE[i*2+j*ny*2].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2].nd[1].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[1].number=i+1+(nx+1)*(j+1); //NO.1
//加边缘裂纹(当划分为4*4时)在单元i==0&&j==1 ,及i==0, j==2两单元加裂纹
/// if(i==0&&j==1){
// pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
/// pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1);//-0.2;
// pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
// }
// else{
pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
// }
//for the second triangle in the rectangle
pE[i*2+j*ny*2+1].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2+1].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[0].number=i+j*(nx+1); //NO.0
pE[i*2+j*ny*2+1].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[1].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[1].number=i+j*(nx+1)+1; //NO.1
// if(i==0&&j==2){
// pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
// pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1)+0.2;
// pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
// }
// else{
pE[i*2+j*ny*2+1].nd[2].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2+1].nd[2].number=i+1+(nx+1)*(j+1); //NO.2
// }
}
System.out.println("计算基函数系数值...\n");
int j=1,k=2;
for(int id=0;id<nx*ny*2;id++)
{
for( int i=0;i<3;i++)
{
if(i==0) { j=1;k=2;}
else if(i==1){ j=2;k=0;}
else if(i==2){ j=0;k=1;}
pE[id].A=( (pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-(pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x) )/2.0;
D=2.0*pE[id].A;
pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;
pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;
pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;
}
}System.out.println("OK!\n");
System.out.println("计算单元有限元特征式系数矩阵...\n");
//--------------------------------------计算单元有限元特征式系数矩阵可以不再分两个三角形循环
//Aij即是Et/[4*(1-u*u)*A]*一个2*2Matrix
int l,m;
for(int i=0;i<nx*ny*2;i++)
{ for(l=0;l<3;l++)
for(m=0;m<3;m++) // Respaired
{
pE[i].Aij[2*l][2*m]=( pE[i].b[l]*pE[i].b[m]+ pE[i].c[l]*pE[i].c[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
pE[i].Aij[2*l][2*m+1]=( pE[i].b[l]*pE[i].c[m]*u + pE[i].c[l]*pE[i].b[m]*(1-u)/2)*E*t/(2*(1-u*u));
pE[i].Aij[2*l+1][2*m]=( pE[i].c[l]*pE[i].b[m]*u + pE[i].b[l]*pE[i].c[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
pE[i].Aij[2*l+1][2*m+1]=( pE[i].c[l]*pE[i].c[m]*u + pE[i].b[l]*pE[i].b[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
}
System.out.println("\n");
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
for(int i=0;i<nx*ny*2;i++)
{ for(l=0;l<3;l++)
for(m=0;m<3;m++) // Respaired
{
if(l==m)
{
pE[i].M[2*l][2*m]=(rou)*t*pE[i].A/6;
pE[i].M[2*l][2*m+1]=0.0;
pE[i].M[2*l+1][2*m]=0.0;
pE[i].M[2*l+1][2*m+1]=(rou)*t*pE[i].A/6;
pE[i].C[2*l][2*m]=(Y)*t*pE[i].A/6;
pE[i].C[2*l][2*m+1]=0.0;
pE[i].C[2*l+1][2*m]=0.0;
pE[i].C[2*l+1][2*m+1]=(Y)*t*pE[i].A/6;
}else
{
pE[i].M[2*l][2*m]=(rou)*t*pE[i].A/12;
pE[i].M[2*l][2*m+1]=0.0;
pE[i].M[2*l+1][2*m]=0.0;
pE[i].M[2*l+1][2*m+1]=(rou)*t*pE[i].A/12;
pE[i].C[2*l][2*m]=(Y)*t*pE[i].A/12;
pE[i].C[2*l][2*m+1]=0;
pE[i].C[2*l+1][2*m]=0;
pE[i].C[2*l+1][2*m+1]=Y*t*pE[i].A/12;
}
}
System.out.println("\n");
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
System.out.println("打印单元特征式系数矩阵...\n");
//fSystem.out.println(wfp,"计算得各结点上的温度值为:\n");
for(int i=0;i<nx*ny*2;i++)
{
System.out.println("计算得单元"+i+"上的K矩阵为:\n");
for(l=0;l<6;l++)
{
for(m=0;m<6;m++)
{
// if(abs(pE[i].Aij[l][m]-0.0)<0.1)
// System.out.println("0.00\t");
// else
System.out.println("\t"+pE[i].Aij[l][m]);
}
System.out.println("\n");
// System.out.println("%6.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%6.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
}
//fSystem.out.println(wfp,"%d %f\n",i+1,pMf[i]);
// System.out.println("写特征式系数矩阵到文件OK!\n");
// fclose(wfp);
System.out.println("计算积分值,填充到f函数向量数组...\n");
//---------------------------计算三角形单元f函数向量不用多重积分公式,直接代用单元三角形面积的三分之一
int idx=0;
for(int i=0;i<2*nx*ny;i++)
{ for(idx=0;idx<3;idx++)
{
pE[i].fi[idx]=(0.5)*(Lx/nx)*(Ly/ny); // Respaired
}
}
// System.out.println("pE[0].fi[0]=%f\n",pE[0].fi[0]);
// System.out.println("pE[0].fi[1]=%f\n",pE[0].fi[1]);
// System.out.println("pE[0].fi[2]=%f\n",pE[0].fi[2]);
// System.out.println("OK!\n");
//单元矩阵元素累加到总体矩阵相应的位置上,共四个节点,两个单元
// 2 |-------|3
// | I . |
// | . II |
// 0 |-------|1
System.out.println("单元矩阵元素累加到总体矩阵相应的位置上...\n");
for(idx=0;idx<nx*ny*2;idx++)
for(int i=0;i<3;i++)
{
for(j=0;j<3;j++)
{
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].Aij[2*i][2*j];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=pE[idx].Aij[2*i][2*j+1];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].Aij[2*i+1][2*j];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].Aij[2*i+1][2*j+1];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].M[2*i][2*j];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=pE[idx].M[2*i][2*j+1];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].M[2*i+1][2*j];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].M[2*i+1][2*j+1];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].C[2*i][2*j];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=pE[idx].C[2*i][2*j+1];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].C[2*i+1][2*j];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].C[2*i+1][2*j+1];
// pMatrix[pE[idx].nd[i].number*iNode+pE[idx].nd[j].number] += (totalM[pE[idx].nd[i].number*iNode+pE[idx].nd[j].number]+totalC[pE[idx].nd[i].number*iNode+pE[idx].nd[j].number]);
// pMatrix[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+1]+=( totalM[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+1]+totalC[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+1]);
// pMatrix[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+iNode] += (totalM[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+2*iNode]+totalC[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+iNode]);
// pMatrix[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+iNode+1] +=(totalM[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+2*iNode+1]+totalC[(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+iNode+1]);
pMf[pE[idx].nd[i].number ]+=pE[idx].fi[i];
}
}
/*
for(i=0;i<iNode*2;i++)
{
for(j=0;j<iNode*2;j++)
{
System.out.println("%6d\t",pMatrix[i*iNode*2+j]);
}
System.out.println("%6d\t",pMf[i*iNode*2+j]);
System.out.println("\n");
}
*/
// System.out.println("总体矩阵到文件...\n");
//文件操作
// FILE *wfp;
// if((wfp=fopen("gobalMatrix.txt","w"))==NULL)
// System.out.println("Cann't open the file... ");
for(idx=0;idx<nx*ny*2;idx++)
{
System.out.println("\r\n");
for(int i=0;i<3;i++)
{
for(int h=0;h<3;h++)
{
System.out.print("%6.2f\t"+pE[idx].Aij[2*i][2*h]);
System.out.print("%6.2f\t"+pE[idx].Aij[2*i][2*h+1]);
}
System.out.println();
// System.out.println("%6.2f\t"+pE[idx].Aij[2*i+1][2*2]);
// System.out.println("%6.2f\t"+pE[idx].Aij[2*i+1][2*2+1]);
}
}
//fSystem.out.println(wfp,"%d %f\n",i+1,pMf[i]);
System.out.println("\r\n总体矩阵到文件OK!\n");
//fclose(wfp);
//double dBig=pow(10,20); //边界条件对角线扩大法处理所用的大数
// double Ur=0.0; //边界条件1(边界条件2通过Calerkin弱解表达式自动满足)
double Ur=1.23 ; //外力,牛
//jimmy 上下拉伸
//上端2,3节点,Y方向拉力 索引为5,7
for(int i=0;i<nx+1;i++)
{
pMf[(nx+1)*nx*2+2*i+1]=Ur;
}
//下端0,1节点,Y方向拉力为负 ,索引为1,3
// 位移 ,索引
//u0, 0
// v0, 1
/// u1, 2
// v1, 3
// u2, 4
// v2, 5
// u3, 6
// v3, 7
//设置为y方向为0位移,参考王勋成的书
for(j=0;j<nx+1;j++)
{
pMatrix[(2*j)*2*iNode+2*j]*=dBig;
pMatrix[(2*j+1)*2*iNode+2*j+1]*=dBig;
pMf[2*j]=0;
pMf[2*j+1]=0;
}
//右端 for(j=0;j<nx+1;j++)
// { i=(nx+1)*(j+1)-1;
// pMatrix[i*iNode+i]*=dBig;
// pMf[i]=pMatrix[i*iNode+i]*Ur;
//}
System.out.println("调用全选主元高斯消去法函数解方程组...\n");
//if((wfp=fopen("result.txt","w"))==NULL)
// System.out.println("Cann't open the file... ");
//第一步求t=0加速度的,结果在pMf中
//假设t=0,速度和位移为0
for(int i=0;i<2*iNode;i++)
tmp[i]=pMf[i];
Gauss(pMatrix,pMf,2*iNode);
System.out.println("\n节点位移: x方向, y方向 ");
for(int i=0;i<2*iNode;i++)
{
System.out.println("\n"+i/2+" "+pMf[i]+" "+pMf[i+1]);
i++;
//tmp[i]=pMf[i];
}
// fSystem.out.println(wfp,"计算得各结点上的位移值为:\r\n节点,X方向位移,y方向位移\n");
// fSystem.out.println(wfp,"%d\t%8.2f\t%8.2f",iNode,pMf[2*iNode-2],pMf[2*iNode-1]);
// fclose(wfp);
//假设初始加速度为确定值
for(int i=(2*iNode-1);i<2*iNode;i++)
dalta2[i]=0.1;
//增加时间步长的内容
// for(int step=1;step<T/dataT;step++)
// {
// //if((wfp=fopen("result2.txt","w+"))==NULL)
// //System.out.println("Cann't open the file... ");
// for(int i=0;i<2*iNode;i++)
// for(j=0;j<2*iNode;j++)
// {
// pMf[i]=totalM[i*2*iNode+j]*(c0*dalta[j]+c2*dalta1[j]+c3*dalta2[j])+totalC[i*2*iNode+j]*(c1*dalta[j]+c4*dalta1[j]+c5*dalta2[j]);
// }
//
// Gauss(pMatrix,pMf,iNode); //调用全选主元高斯消去法函数解方程组
// for(int i=0;i<2*iNode;i++)
// {
//
// dalta2[i]=c0*(pMf[i]-dalta[i])-c2*dalta1[i]-c3*dalta2[i];
// dalta1[i]=dalta1[i]+c6*dalta2[i]+c7*dalta2[i];
// dalta[i]=pMf[i];
//
// }
//
// System.out.println(iNode+"\t"+pMf[2*iNode-2]+"\t%8.7f"+pMf[2*iNode-1]);
//
// System.out.println("\r\n");
// // fclose(wfp);
//
// System.out.println("step="+step);
//
// }
}
catch(Exception e)
{
e.printStackTrace();
System.out.println("Error occured...\n");
}
//释放总体矩阵和三角形单元数组占用内存
//free(pMf); free(pE); free(pMatrix);
System.out.println("Hello World!\n");
return ;
}
}