# FEM implementation on tet10 element stiffness matrix:(code snippet provided) [duplicate]

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

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

• This site works a little bit different. You want to post a relevant piece of code, which would allow reproducing the key numerical result. Also, you want to formulate the relevant math in your question using MathJax in the first place. And formulate the question, not "I have a bug in my code, please debug it for me". NB: First, I would check if your stiffness matrix is just scaled differently with some scalar $\alpha$. – Anton Menshov May 8 '17 at 7:31
• Thanks I have provided the code snippet for the same. Regarding scaling it is surprising that it works if I provide a constant function to give me the volume of the tet. – kaush May 8 '17 at 10:49
• If you have a negative determinant it's possible that your ordering is inverted – nicoguaro May 8 '17 at 12:59
• After fixing your problem with the jacobian calculation, as a debug step, add a line for calculation of the element volume to your integration point loop. If you can't calculate the volume of the element, you clearly aren't going to get the correct stiffness matrix. – Bill Greene May 8 '17 at 13:15
• Thanks a lot! It works and I get the correct answer. Have corrected the code above. Also it is worthwhile to do the above calculations with matlab before going to c++. So I am also posting a matlab snippet of the same. – kaush May 9 '17 at 7:38

Here is the matlab snippet for the same.

clc;
node = [2 3 4
6 3 2
2 5 1
4 3 6
4 3 3
4 4 1.5
2 4 2.5
3 3 5
5 3 4
3 4 3.5];
Y=480; %Youngs Modulus
nu=1/3; % Poisson ratio
alfa=(5.0+3.0*sqrt(5.0))/20; %Gauss Points
beta =(5.0-sqrt(5.0))/20;    %Gauss Points
weight=0.25; %weights for integration points
GaussPoints=[alfa beta beta beta;beta alfa beta beta;beta beta alfa  beta;beta beta beta alfa];
%Elasticity Matrix
E=[1-nu nu nu 0 0 0;nu 1-nu nu 0 0 0;nu nu 1-nu 0 0 0;0 0 0 0.5-nu 0 0;0 0 0 0 0.5-nu 0;0 0 0 0 0 0.5-nu];
E=Y/((1+nu)*(1-2*nu))*E;
%Stiffness Matrix
K=zeros(30,30);

for j=1:4
xi=GaussPoints(:,j);
jx1=4.0*(node(1,1)*(xi(1)-0.25)+node(5,1)*xi(2)+node(7,1)*xi(3)+node(8,1)*xi(4));
jy1=4.0*(node(1,2)*(xi(1)-0.25)+node(5,2)*xi(2)+node(7,2)*xi(3)+node(8,2)*xi(4));
jz1=4.0*(node(1,3)*(xi(1)-0.25)+node(5,3)*xi(2)+node(7,3)*xi(3)+node(8,3)*xi(4));

jx2=4.0*(node(5,1)*xi(1)+node(2,1)*(xi(2)-0.25)+node(6,1)*xi(3)+node(9,1)*xi(4));
jy2=4.0*(node(5,2)*xi(1)+node(2,2)*(xi(2)-0.25)+node(6,2)*xi(3)+node(9,2)*xi(4));
jz2=4.0*(node(5,3)*xi(1)+node(2,3)*(xi(2)-0.25)+node(6,3)*xi(3)+node(9,3)*xi(4));

jx3=4.0*(node(7,1)*xi(1)+node(6,1)*xi(2)+node(3,1)*(xi(3)-0.25)+node(10,1)*xi(4));
jy3=4.0*(node(7,2)*xi(1)+node(6,2)*xi(2)+node(3,2)*(xi(3)-0.25)+node(10,2)*xi(4));
jz3=4.0*(node(7,3)*xi(1)+node(6,3)*xi(2)+node(3,3)*(xi(3)-0.25)+node(10,3)*xi(4));

jx4=4.0*(node(8,1)*xi(1)+node(9,1)*xi(2)+node(10,1)*xi(3)+node(4,1)*(xi(4)-0.25));
jy4=4.0*(node(8,2)*xi(1)+node(9,2)*xi(2)+node(10,2)*xi(3)+node(4,2)*(xi(4)-0.25));
jz4=4.0*(node(8,3)*xi(1)+node(9,3)*xi(2)+node(10,3)*xi(3)+node(4,3)*(xi(4)-0.25));

J=[1 1 1 1;jx1 jx2 jx3 jx4;jy1 jy2 jy3 jy4;jz1 jz2 jz3 jz4];

Jdet=det(J);
Jinv=inv(J);
Iaug=[0 0 0;1 0 0;0 1 0;0 0 1];
P=Jinv*Iaug;

a1=P(1,1);
a2=P(2,1);
a3=P(3,1);
a4=P(4,1);

b1=P(1,2);
b2=P(2,2);
b3=P(3,2);
b4=P(4,2);

c1=P(1,3);
c2=P(2,3);
c3=P(3,3);
c4=P(4,3);

Nfx=[(4.0*xi(1)-1)*a1 (4.0*xi(2)-1)*a2  (4.0*xi(3)-1)*a3 (4.0*xi(4)-1)*a4 ...
4.0*(a1*xi(2)+a2*xi(1)) 4.0*(a2*xi(3)+a3*xi(2)) 4.0*(a1*xi(3)+a3*xi(1)) ...
4.0*(a1*xi(4)+a4*xi(1)) 4.0*(a2*xi(4)+a4*xi(2)) 4.0*(a3*xi(4)+a4*xi(3))];

Nfy=[(4.0*xi(1)-1)*b1 (4.0*xi(2)-1)*b2 (4.0*xi(3)-1)*b3 (4.0*xi(4)-1)*b4 ...
4.0*(b1*xi(2)+b2*xi(1)) 4.0*(b2*xi(3)+b3*xi(2)) 4.0*(b1*xi(3)+b3*xi(1)) ...
4.0*(b1*xi(4)+b4*xi(1)) 4.0*(b2*xi(4)+b4*xi(2)) 4.0*(b3*xi(4)+b4*xi(3))];

Nfz=[(4.0*xi(1)-1)*c1 (4.0*xi(2)-1)*c2 (4.0*xi(3)-1)*c3 (4.0*xi(4)-1)*c4 ...
4.0*(c1*xi(2)+c2*xi(1)) 4.0*(c2*xi(3)+c3*xi(2)) 4.0*(c1*xi(3)+c3*xi(1)) ...
4.0*(c1*xi(4)+c4*xi(1)) 4.0*(c2*xi(4)+c4*xi(2)) 4.0*(c3*xi(4)+c4*xi(3))];

B=[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 Nfx(10) 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 0 Nfy(10) 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) 0 0 Nfz(10);
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 Nfy(10) Nfx(10) 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) 0 Nfz(10) Nfy(10);
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) Nfz(10) 0 Nfx(10)];

K=K+weight*(B'*E*B*Jdet/6.0);
Jdet
end

• You can accept your own answer. And it would be really nice if you formulate a sentence of what has been fixed. – Anton Menshov May 9 '17 at 10:18
• It was just the ordering of the most edge nodes used to compute the jacobian which was messed up and that led to the negative volumes. – kaush May 10 '17 at 6:54
• The ordering of the nodes depends on the mesh generator used. – kaush May 10 '17 at 13:01
• Hi kaush and welcome to scicomp! I highly encourage you to include this information as part of your question instead of the "answer" section. As it stands, this is not an "answer" to your own question. – Paul May 23 '17 at 2:02
• I have decided to start a new question as I thought it as a different topic. The above was related to computing the stiffness matrix and the code gives the correct value of the stiffness matrix as the reference. scicomp.stackexchange.com/questions/26957/… – kaush May 23 '17 at 9:23