I am trying to implement a fem code on tet10 elements. I do not prefer to use open source at the moment as I would like to have the basic feel over the algorithm. I closely follow the lecture notes for tet10 implementation at
http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf
However when I test my code on the example tetrahedron given above I get a very different stiffness matrix than given in the document. I however get a match on the trend of the Eigen values for my stiffness matrix. I have included a test code here to just compute the basic stiffness matrix.
#include<cstdlib>
#include<iostream>
#include<fstream>
#include<ctime>
#include<cmath>
#include<string>
#include<vector>
#include<unistd.h>
#include <iterator>
#include <fstream>
using namespace std;
#include<Eigen/Core>
#include<Eigen/Dense>
#include<Eigen/Eigenvalues>
#include<Eigen/SVD>
#include<Eigen/Sparse>
#include<Eigen/SparseCholesky>
#include<Eigen/SparseLU>
using namespace Eigen;
typedef Eigen::SparseMatrix<double> SpMat;
struct NodeInfo
{
double x,y,z;
};
struct ElementInfo
{
int node[10];
double volume, area[4],nx[4],ny[4],nz[4];
};
struct constants
{
double YoungsModulus, PoissonRatio, Density, gravity;
};
std::istream& operator>>(std::istream& is, ElementInfo& element)
{
is >> element.node[0] >> element.node[1] >> element.node[2] >> element.node[3]
>> element.node[4] >> element.node[5] >> element.node[6] >> element.node[7]
>> element.node[8] >> element.node[9] >> element.volume
>> element.nx[0] >> element.ny[0] >>element.nz[0] >> element.area[0]
>> element.nx[1] >> element.ny[1] >>element.nz[1] >> element.area[1]
>> element.nx[2] >> element.ny[2] >>element.nz[2] >> element.area[2]
>> element.nx[3] >> element.ny[3] >>element.nz[3] >> element.area[3];
return is;
}
std::istream& operator>>(std::istream& is, NodeInfo& node)
{
is >> node.x >>node.y >> node.z ;
return is;
}
void ComputeElasticityMat(Eigen::MatrixXd& El, constants& M);
void ComputeElementStiffness(ElementInfo& ithEle, vector<NodeInfo>& node,
Eigen::MatrixXd& K, Eigen::MatrixXd& E, double (*eta)[4], double weight);
int main()
{
char filename[] = "elements.dat";
char filename1[] ="nodes.dat";
int nelem, nnodes, TotalDOF;
int ndf=3;//number of degrees of freedom at each node
int npe=10;//number of points per element
// Gauss Integration points
double alfa=(5.0+3.0*sqrt(5.0))/20;
double beta =(5.0-sqrt(5.0))/20;
double weight=1.0/4.0;
double GaussPoints[4][4]={{alfa,beta,beta,beta},{beta,alfa,beta,beta}, {beta,beta,alfa,beta},{beta,beta,beta,alfa}};
// To define the element properties
constants MaterialProperties={.YoungsModulus=480.0,.PoissonRatio=1.0/3.0, .Density=2830, .gravity=9.81};
// Read the data file for the elements and the nodes
std::vector<ElementInfo> Elem;
std::vector<NodeInfo> Node;
std::ifstream ifs(filename), ifs1(filename1);
ifs>>nelem;
cout<<"The total number of elements is="<<nelem<<endl;
if (ifs) {
std::copy(std::istream_iterator<ElementInfo>(ifs),
std::istream_iterator<ElementInfo>(),
std::back_inserter(Elem));
}
else {
std::cerr << "Couldn't open " << filename << " for reading\n";
}
cout<<"Element file read sucessfuly "<<endl;
if (ifs1) {
std::copy(std::istream_iterator<NodeInfo>(ifs1),
std::istream_iterator<NodeInfo>(),
std::back_inserter(Node));
}
else {
std::cerr << "Couldn't open " << filename1 << " for reading\n";
}
cout<<"Node file read sucessfuly "<<endl;
cout<<"The total number of Nodes is="<<Node.size() <<endl;
nnodes=Node.size();
//To form the Elasticity Matrix
TotalDOF=nnodes*ndf;
cout<<"The Total number of Degrees of freedom=:\n"<<TotalDOF<<endl;
MatrixXd EelemMat(6,6), KelemMat(30,30), KGl(TotalDOF, TotalDOF);
VectorXd BodyForce(30), BFGl(TotalDOF), SurfaceForce(30);
EelemMat.setZero(6,6);
ComputeElasticityMat(EelemMat,MaterialProperties);
cout<<"Elasticity Matrix=:\n"<<EelemMat<<endl;
//To find the Stiffness Matrix
if(nelem!=Elem.size()){
cout<<"Something wrong in element file...please check!"<<endl;
exit(0);
}
for(int i=0;i<nelem;i++)
{
KelemMat.setZero(30,30);
ComputeElementStiffness(Elem[i],Node,KelemMat,
EelemMat,GaussPoints,weight);
cout<<"Stiffness Matrix="<<KelemMat<<endl;
}
JacobiSVD<MatrixXd> svd(KelemMat, ComputeThinU | ComputeThinV);
return 0;
}
void ComputeElasticityMat(Eigen::MatrixXd& El, constants& M)
{
double K=M.YoungsModulus/((1.0+M.PoissonRatio)* (1.0-2.0*M.PoissonRatio));
El(0,0)=1.0-M.PoissonRatio;
El(0,1)=M.PoissonRatio;
El(0,2)=M.PoissonRatio;
El(1,0)=M.PoissonRatio;
El(1,1)=1.0-M.PoissonRatio;
El(1,2)=M.PoissonRatio;
El(2,0)=M.PoissonRatio;
El(2,1)=M.PoissonRatio;
El(2,2)=1.0-M.PoissonRatio;
El(3,3)=0.5-M.PoissonRatio;
El(4,4)=El(3,3);
El(5,5)=El(3,3);
El=El*K;
return;
}
void ComputeElementStiffness(ElementInfo& ithEle, vector<NodeInfo>& node, Eigen::MatrixXd& K,
Eigen::MatrixXd& E, double (*eta)[4], double weight)
{
double jx1,jx2,jx3,jx4,jy1,jy2,jy3,jy4,jz1,jz2,jz3,jz4,
a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,Jdet;
VectorXd Nfx(10),Nfy(10),Nfz(10),xi(4);
MatrixXd B(6,30), Bt;
MatrixXd J(4,4), P, Iaug(4,3),Jinv(4,4);
double v01,v02,v03,v04,V;
for(int intPoints=0;intPoints<4;intPoints++)
{
for(int j=0;j<4;j++)
xi(j) = eta[intPoints][j];
jx1=4.0*(node[ithEle.node[0]].x*(xi(0)-0.25)+node[ithEle.node[4]].x*xi(1)+node[ithEle.node[6]].x*xi(2)+node[ithEle.node[7]].x*xi(3));
jy1=4.0*(node[ithEle.node[0]].y*(xi(0)-0.25)+node[ithEle.node[4]].y*xi(1)+node[ithEle.node[6]].y*xi(2)+node[ithEle.node[7]].y*xi(3));
jz1=4.0*(node[ithEle.node[0]].z*(xi(0)-0.25)+node[ithEle.node[4]].z*xi(1)+node[ithEle.node[6]].z*xi(2)+node[ithEle.node[7]].z*xi(3));
jx2=4.0*(node[ithEle.node[4]].x*xi(0)+node[ithEle.node[1]].x*(xi(1)-0.25)+node[ithEle.node[5]].x*xi(2)+node[ithEle.node[8]].x*xi(3));
jy2=4.0*(node[ithEle.node[4]].y*xi(0)+node[ithEle.node[1]].y*(xi(1)-0.25)+node[ithEle.node[5]].y*xi(2)+node[ithEle.node[8]].y*xi(3));
jz2=4.0*(node[ithEle.node[4]].z*xi(0)+node[ithEle.node[1]].z*(xi(1)-0.25)+node[ithEle.node[5]].z*xi(2)+node[ithEle.node[8]].z*xi(3));
jx3=4.0*(node[ithEle.node[6]].x*xi(0)+node[ithEle.node[5]].x*xi(1)+node[ithEle.node[2]].x*(xi(2)-0.25)+node[ithEle.node[9]].x*xi(3));
jy3=4.0*(node[ithEle.node[6]].y*xi(0)+node[ithEle.node[5]].y*xi(1)+node[ithEle.node[2]].y*(xi(2)-0.25)+node[ithEle.node[9]].y*xi(3));
jz3=4.0*(node[ithEle.node[6]].z*xi(0)+node[ithEle.node[5]].z*xi(1)+node[ithEle.node[2]].z*(xi(2)-0.25)+node[ithEle.node[9]].z*xi(3));
jx4=4.0*(node[ithEle.node[7]].x*xi(0)+node[ithEle.node[8]].x*xi(1)+node[ithEle.node[9]].x*xi(2)+node[ithEle.node[3]].x*(xi(3)-0.25));
jy4=4.0*(node[ithEle.node[7]].y*xi(0)+node[ithEle.node[8]].y*xi(1)+node[ithEle.node[9]].y*xi(2)+node[ithEle.node[3]].y*(xi(3)-0.25));
jz4=4.0*(node[ithEle.node[7]].z*xi(0)+node[ithEle.node[8]].z*xi(1)+node[ithEle.node[9]].z*xi(2)+node[ithEle.node[3]].z*(xi(3)-0.25));
J.row(0)<<1,1,1,1;
J.row(1)<<jx1,jx2,jx3,jx4;
J.row(2)<<jy1,jy2,jy3,jy4;
J.row(3)<<jz1,jz2,jz3,jz4;
Jdet=J.determinant();
Jinv=J.inverse();
Iaug<<0,0,0,
1,0,0,
0,1,0,
0,0,1;
P=Jinv*Iaug;
a1=P(0,0);
a2=P(1,0);
a3=P(2,0);
a4=P(3,0);
b1=P(0,1);
b2=P(1,1);
b3=P(2,1);
b4=P(3,1);
c1=P(0,2);
c2=P(1,2);
c3=P(2,2);
c4=P(3,2);
Nfx << (4.0*xi(0)-1)*a1, (4.0*xi(1)-1)*a2, (4.0*xi(2)-1)*a3, (4.0*xi(3)-1)*a4,
4.0*(a1*xi(1)+a2*xi(0)), 4.0*(a2*xi(2)+a3*xi(1)), 4.0*(a1*xi(2)+a3*xi(0)),
4.0*(a1*xi(3)+a4*xi(0)), 4.0*(a2*xi(3)+a4*xi(1)), 4.0*(a3*xi(3)+a4*xi(2));
Nfy << (4.0*xi(0)-1)*b1, (4.0*xi(1)-1)*b2, (4.0*xi(2)-1)*b3, (4.0*xi(3)-1)*b4,
4.0*(b1*xi(1)+b2*xi(0)), 4.0*(b2*xi(2)+b3*xi(1)), 4.0*(b1*xi(2)+b3*xi(0)),
4.0*(b1*xi(3)+b4*xi(0)), 4.0*(b2*xi(3)+b4*xi(1)), 4.0*(b3*xi(3)+b4*xi(2));
Nfz << (4.0*xi(0)-1)*c1, (4.0*xi(1)-1)*c2, (4.0*xi(2)-1)*c3, (4.0*xi(3)-1)*c4,
4.0*(c1*xi(1)+c2*xi(0)), 4.0*(c2*xi(2)+c3*xi(1)), 4.0*(c1*xi(2)+c3*xi(0)),
4.0*(c1*xi(3)+c4*xi(0)), 4.0*(c2*xi(3)+c4*xi(1)), 4.0*(c3*xi(3)+c4*xi(2));
B.row(0) << Nfx(0),0,0, Nfx(1),0,0, Nfx(2),0,0, Nfx(3),0,0, Nfx(4),0,0, Nfx(5),0,0, Nfx(6),0,0, Nfx(7),0,0, Nfx(8),0,0, Nfx(9),0,0;
B.row(1) << 0,Nfy(0),0, 0,Nfy(1),0, 0,Nfy(2),0, 0,Nfy(3),0, 0,Nfy(4),0, 0,Nfy(5),0, 0,Nfy(6),0, 0,Nfy(7),0, 0,Nfy(8),0, 0,Nfy(9),0;
B.row(2) << 0,0,Nfz(0), 0,0,Nfz(1), 0,0,Nfz(2), 0,0,Nfz(3), 0,0,Nfz(4), 0,0,Nfz(5), 0,0,Nfz(6), 0,0,Nfz(7), 0,0,Nfz(8), 0,0,Nfz(9);
B.row(3) << Nfy(0),Nfx(0),0, Nfy(1),Nfx(1),0, Nfy(2),Nfx(2),0, Nfy(3),Nfx(3),0, Nfy(4),Nfx(4),0, Nfy(5),Nfx(5),0, Nfy(6),Nfx(6),0,
Nfy(7),Nfx(7),0,Nfy(8),Nfx(8),0, Nfy(9),Nfx(9),0;
B.row(4) << 0,Nfz(0),Nfy(0), 0,Nfz(1),Nfy(1), 0,Nfz(2),Nfy(2), 0,Nfz(3),Nfy(3), 0,Nfz(4),Nfy(4), 0,Nfz(5),Nfy(5), 0,Nfz(6),Nfy(6),
0,Nfz(7),Nfy(7), 0,Nfz(8),Nfy(8), 0,Nfz(9),Nfy(9);
B.row(5) << Nfz(0),0,Nfx(0), Nfz(1),0,Nfx(1), Nfz(2),0,Nfx(2), Nfz(3),0,Nfx(3), Nfz(4),0,Nfx(4), Nfz(5),0,Nfx(5), Nfz(6),0,Nfx(6),
Nfz(7),0,Nfx(7), Nfz(8),0,Nfx(8), Nfz(9),0,Nfx(9);
Bt = B.transpose();
K += weight*(Bt*E*B*Jdet/6.0);
cout<<"Jdet="<<Jdet<<endl;
}
return;
}
The element file (elements.dat) for the problem is
1
0 1 2 3 4 5 6 7 8 9
24.0 -0.6295797707 -0.7709500359 -0.0962567111 0.0018562616
-0.5930703716 0.2172077372 -0.7752988670 0.0051374961
0.9488273288 -0.1721170560 0.2647686142 0.0037463536
0.1956924749 0.2842456298 0.9385674601 0.0033773339
The node file (nodes.dat)
2.0000000000 3.0000000000 4.0000000000
6.0000000000 3.0000000000 2.0000000000
2.0000000000 5.0000000000 1.0000000000
4.0000000000 3.0000000000 6.0000000000
4.0000000000 3.0000000000 3.0000000000
4.0000000000 4.0000000000 1.5000000000
2.0000000000 4.0000000000 2.5000000000
3.0000000000 3.0000000000 5.0000000000
5.0000000000 3.0000000000 4.0000000000
3.0000000000 4.0000000000 3.5000000000
The correct stiffness matrix as given in figure 10.8 is just for the first three values
KelemMat=[447, 324,72.....]
whereas I get
KelemMat=[ 70.617 134.209 -108.193 .... ]
One more thing is it OK to get a negative determinant at a Gaussian point..I do get one