Jacobian for 6-noded triangle in 3D to calculate the area

I would like to calculate the surface area of a 6-noded triangle element, i.e., the face of a 10-noded tetrahedral element in 3D space. A typical solution is to calculate the surface integral of the Jacobian of such an element, however I can only find derivations of the shape functions and thus the Jacobian for the 2D element.

For a 3-noded triangle, this is not an issue, as all three vertices of the triangle are always co-planar, thus a transformation into a 2D element an be found, as shown here: Finding Shape Functions for a Triangle in 3D coordinate space .

However if mid-nodes are added, they might not necessarily be co-planar with the other vertices anymore, thus I would believe no such transformation can be used.

Is there a special transformation to do this? Do I have to calculate the shape functions and Jacobian for the 3D case on my own or do I maybe need to resort to the shape functions of a 10-noded tetrahedron and integrate over the surfaces?

Edit: I tried to calculate the 3D Variant of the Jacobian on my own, but I do not get to the desired results.

I assumed the following shape functions: $$N_1 = \zeta_1(2\zeta_1 - 1)$$ $$N_2 = \zeta_2(2\zeta_2 - 1)$$ $$N_3 = \zeta_3(2\zeta_3 - 1)$$ $$N_4=4\zeta_1\zeta_2$$ $$N_5=4\zeta_2\zeta_3$$ $$N_6=4\zeta_3\zeta_1$$.

Then I created the partial derivatives w.r.t. the natural coordinates, which gets me the following matrix:

$$\nabla\mathbf{N} = \left(\begin{array}{cccccc} 4 \zeta _1-1 & 0 & 0 & 4 \zeta _2 & 0 & 4 \zeta _3 \\ 0 & 4 \zeta _2-1 & 0 & 4 \zeta _1 & 4 \zeta _3 & 0 \\ 0 & 0 & 4 \zeta _3-1 & 0 & 4 \zeta _2 & 4 \zeta _1 \\ \end{array} \right)$$

I assume the Jacobian matrix can then be calculated as $$\mathbf{J}=\nabla\mathbf{N}\cdot\mathbf{X}$$, where $$\mathbf{X}$$ is the vector of node coordinates: $$\mathbf{J}=4\left( \begin{array}{ccc} \left(\zeta _1-\frac{1}{4}\right) x_1+\zeta _2 x_4+\zeta _3 x_6 & \left(\zeta _1-\frac{1}{4}\right) y_1+\zeta _2 y_4+\zeta _3 y_6 & \left(\zeta _1-\frac{1}{4}\right) z_1+\zeta _2 z_4+\zeta _3 z_6 \\ \zeta _1 x_4+\left(\zeta _2-\frac{1}{4}\right) x_2+\zeta _3 x_5 & \zeta _1 y_4+\left(\zeta _2-\frac{1}{4}\right) y_2+\zeta _3 y_5 & \zeta _1 z_4+\left(\zeta _2-\frac{1}{4}\right) z_2+\zeta _3 z_5 \\ \zeta _1 x_6+\zeta _2 x_5+\left(\zeta _3-\frac{1}{4}\right) x_3 & \zeta _1 y_6+\zeta _2 y_5+\left(\zeta _3-\frac{1}{4}\right) y_3 & \zeta _1 z_6+\zeta _2 z_5+\left(\zeta _3-\frac{1}{4}\right) z_3 \\ \end{array} \right)$$

The area of the triangle should then be given as $$\frac{1}{2}\mathrm{det}(\mathbf{J})$$ evaluated at $$\zeta_1=\zeta_2=\zeta_3=\frac{1}{3}$$, using single point Gauss quadrature. However, for my test cases this does not work out...

What do I miss here? Again, the same question: Can I not use the shape functions derived for the 2D Element?

edit2: Here are the derivations from the answer of Bill Greene:

$$\mathbf{r}_{\xi_1}=\left( \begin{array}{c} 2 \zeta _1 x_1+\left(2 \zeta _1-1\right) x_1+4 \zeta _2 x_4+4 \zeta _3 x_6 \\ 2 \zeta _1 y_1+\left(2 \zeta _1-1\right) y_1+4 \zeta _2 y_4+4 \zeta _3 y_6 \\ 2 \zeta _1 z_1+\left(2 \zeta _1-1\right) z_1+4 \zeta _2 z_4+4 \zeta _3 z_6 \end{array} \right)$$

$$\mathbf{r}_{\xi_2}=\left( \begin{array}{c} 4 \zeta _1 x_4+2 \zeta _2 x_2+\left(2 \zeta _2-1\right) x_2+4 \zeta _3 x_5 \\ 4 \zeta _1 y_4+2 \zeta _2 y_2+\left(2 \zeta _2-1\right) y_2+4 \zeta _3 y_5 \\ 4 \zeta _1 z_4+2 \zeta _2 z_2+\left(2 \zeta _2-1\right) z_2+4 \zeta _3 z_5 \end{array} \right)$$

And some python code:

import numpy as np

# Triangle:
# Nodes 1,2,3 are the corner nodes
X = np.array([
[0,   0, 1],
[2,   0, 1],
[1,   2, 1],
[1,   0, 1],  # 4 between (1, 2)
[1.5, 1, 1],  # 5 between (2, 3)
[0.5, 1, 1],  # 6 between (3, 1)
])

# Calculate area based on corner nodes:
n0, n1, n2 = X[:3, :]
print("Area from corner-nodes:", 0.5 * np.linalg.norm(np.cross(n1 - n0, n2 - n0)))

# Single Gauss Point:
xi1 = xi2 = xi3 = 1/3
w = 1

# Using Jacobian:
NDeriv = np.array([[-1 + 4*xi1, 0,        0,        4*xi2, 0,     4*xi3],
[0,          -1+4*xi2, 0,        4*xi1, 4*xi3, 0    ],
[0,           0,       -1+4*xi3, 0,     4*xi2, 4*xi1]])
Jmat = np.dot(NDeriv, X)
print("Area based on Jacobian:", 0.5 * np.linalg.det(w * Jmat))

# Using tangent vectors
x1, x2, x3, x4, x5, x6 = X[:, 0]
y1, y2, y3, y4, y5, y6 = X[:, 1]
z1, z2, z3, z4, z5, z6 = X[:, 2]
r1 = np.array([
2*x1*xi1 + x1*(2*xi1 - 1) + 4*x4*xi2 + 4*x6*xi3,
2*y1*xi1 + y1*(2*xi1 - 1) + 4*y4*xi2 + 4*y6*xi3,
2*z1*xi1 + z1*(2*xi1 - 1) + 4*z4*xi2 + 4*z6*xi3,
])
r2 = np.array([
4*x4*xi1 + 2*x2*xi2 + x2*(2*xi2 - 1) + 4*x5*xi3,
4*y4*xi1 + 2*y2*xi2 + y2*(2*xi2 - 1) + 4*y5*xi3,
4*z4*xi1 + 2*z2*xi2 + z2*(2*xi2 - 1) + 4*z5*xi3,
])
print("Area based on tangents:", w * np.linalg.norm(np.cross(r1, r2)))


With the result:

Area from corner-nodes: 2.0
Area based on Jacobian: 5.999999999999997
Area based on tangents: 6.5659052011974035


But I noticed that the Jacobian method is dependent on the z-coordiante. If $$z$$ is zero, the area is zero...

Assuming an isoparametric formulation and using the shape functions you show in your post, you can write $$x=N_ix_i, y=N_iy_i$$, and $$x=N_iz_i$$ where $$x_i, y_i$$, and $$z_i$$ are the coordinates of the six nodes and we sum over $$i$$.

From basic analytical geometry, the area of a surface defined in parametric form is given by $$A = \int_S \parallel {\bf r}_{\zeta_1} \times {\bf r}_{\zeta_2} \parallel d\zeta_1 d\zeta_2$$ where $${\bf r}_{\zeta_i}$$ are tangent vectors to the surface along the two parametric directions. As you probably know, only two of the parametric coordinates in your shape functions are independent so write $$\zeta_3=1-\zeta_1-\zeta_2$$ and define

$${\bf r}_{\zeta_1} = \{\partial x/\partial \zeta_1, \partial y/\partial \zeta_1, \partial z/\partial \zeta_1\}$$ and $${\bf r}_{\zeta_2} = \{\partial x/\partial \zeta_2, \partial y/\partial \zeta_2, \partial z/\partial \zeta_2\}$$.

The standard quadrature rules for triangles can be used to evaluate the area integral above numerically.

This simple MATLAB example shows how this algorithm might be implemented.

function cse_triangle_integration
h=2; b=3;
Aanal=b*h/2;
% define nodal coordinates for six node triangle in x-y plane
xyz=[
0,0,0
b,0,0
0,h,0
b/2,0,0
b/2,h/2,0
0,h/2,0
];
A1=calcArea(xyz);

% rotate the triangle by some arbitrary amount
T=rigidRot(pi/3,-pi/6,pi/4);
xyzRot=(T*xyz')';
A2=calcArea(xyzRot);
fprintf('Triangle areas: anal=%f, plane=%f, rotated=%f\n', ...
Aanal, A1, A2);
end

function A=calcArea(xyz)
A=0;
numPts=1;
rs=[1/3 1/3]; wt=.5;
dNdRs=dShapeTri6(rs);
dXyzDrs=dNdRs'*xyz;
dA=norm(cross(dXyzDrs(1,:), dXyzDrs(2,:)));
for i=1:numPts
A = A + wt(i)*dA;
end
end

function ds = dShapeTri6( rs )
L1 = rs(1); L2 = rs(2); L3 = 1 - L1 - L2;
ds = [
1.-4*L3  1.-4*L3
4*L1-1. 0
0   4*L2-1.
4*(1. - L2 - 2.*L1)  -4*L1
4*L2  4*L1
-4*L2  4*(1. - L1 - 2.*L2)];
end

function T=rigidRot(tx,ty,tz)
rotx = @(t) [1 0 0; 0 cos(t) -sin(t) ; 0 sin(t) cos(t)] ;
roty = @(t) [cos(t) 0 sin(t) ; 0 1 0 ; -sin(t) 0  cos(t)] ;
rotz = @(t) [cos(t) -sin(t) 0 ; sin(t) cos(t) 0 ; 0 0 1] ;
T=rotx(tx)*roty(ty)*rotz(tz);
end

• That's an interesting idea but I cannot get this to work and the results are similar to what I calculated with the Jacobian. Maybe I made a typo somewhere... Should the area not also be 1/2 of the integral, as the cross product gives the area of the rectanlge the two vectors form?
– reox
Commented Mar 24 at 13:10
• The 1/2 depends on how you define the weights in your integration rules. Did you try your code on a simple flat, right triangle as a first test? Without seeing the details of your code, I don't have any other suggestions. Do you have some code, preferably matlab, you can post? Commented Mar 24 at 13:19
• I tried it on a flat triangle, with midnodes at 1/2 the distance between the corner nodes. If I calculate the area for the corner nodes only using common triangle equations, I get A=2. However, my Jacobian gives me 11.9 and using the tangent vectors I get 11.03... I'll double check my calculations and update the post above. I calculated the derivates in Mathematica and implemented everything in python.
– reox
Commented Mar 24 at 13:46
• I added a small matlab example to my post. Commented Mar 24 at 14:22
• To be honest, I missed the fact that our definitions 0f $\zeta_1$ and $\zeta_2$ are different ;-) I prefer not to use area coordinates for triangles because, fundamentally, surfaces are parameterized by only two parameters and also I like for triangles and quads to be treated similarly. The main point of my post is that the use of the term "the jacobian" from many FE texts is highly misleading. Better to refer to the fundamental definitions from basic calculus for transforming derivatives from parametric to physical space. Commented Mar 25 at 13:58