Compute outward normal and surface area for 8 noded brick element in FEA

I have a cube which is divided into 8 small cubes by bisecting each edge, I am trying to find out the surface area of each of the faces and the corresponding outward normals for them. This operation is done on a finite element mesh, so I have transformed the cube into the isoparametric form using shape(basis) functions and then tried to extract the area and normals.

Here is the part of the code:

 program polyhedron

DO INPT=1,8  !LOOP OVER 8 SMALL CUBES

! GAUSS POINTS
XII(1,1) = MONE/THREE**HALF
XII(1,2) = MONE/THREE**HALF
XII(1,3) = MONE/THREE**HALF
XII(2,1) = ONE/THREE**HALF
XII(2,2) = MONE/THREE**HALF
XII(2,3) = MONE/THREE**HALF
XII(3,1) = ONE/THREE**HALF
XII(3,2) = ONE/THREE**HALF
XII(3,3) = MONE/THREE**HALF
XII(4,1) = MONE/THREE**HALF
XII(4,2) = ONE/THREE**HALF
XII(4,3) = MONE/THREE**HALF
XII(5,1) = MONE/THREE**HALF
XII(5,2) = MONE/THREE**HALF
XII(5,3) = ONE/THREE**HALF
XII(6,1) = ONE/THREE**HALF
XII(6,2) = MONE/THREE**HALF
XII(6,3) = ONE/THREE**HALF
XII(7,1) = ONE/THREE**HALF
XII(7,2) = ONE/THREE**HALF
XII(7,3) = ONE/THREE**HALF
XII(8,1) = MONE/THREE**HALF
XII(8,2) = ONE/THREE**HALF
XII(8,3) = ONE/THREE**HALF

XI(1) = XII(INPT,1)
XI(2) = XII(INPT,2)
XI(3) = XII(INPT,3)

DO I=1,8
DO J=1,3
dNdxi(I,J) =  ZERO
END DO
END DO

!HEXAHEDRAL SHAPE FUNCTION DERIVATIVES
dNdxi(1,1) =  MONE/EIGHT*(ONE-XI(2))*(ONE-XI(3))
dNdxi(1,2) =  MONE/EIGHT*(ONE-XI(1))*(ONE-XI(3))
dNdxi(1,3) =  MONE/EIGHT*(ONE-XI(1))*(ONE-XI(2))
dNdxi(2,1) =  ONE/EIGHT*(ONE-XI(2))*(ONE-XI(3))
dNdxi(2,2) =  MONE/EIGHT*(ONE+XI(1))*(ONE-XI(3))
dNdxi(2,3) =  MONE/EIGHT*(ONE+XI(1))*(ONE-XI(2))
dNdxi(3,1) =  ONE/EIGHT*(ONE+XI(2))*(ONE-XI(3))
dNdxi(3,2) =  ONE/EIGHT*(ONE+XI(1))*(ONE-XI(3))
dNdxi(3,3) =  MONE/EIGHT*(ONE+XI(1))*(ONE+XI(2))
dNdxi(4,1) =  MONE/EIGHT*(ONE+XI(2))*(ONE-XI(3))
dNdxi(4,2) =  ONE/EIGHT*(ONE-XI(1))*(ONE-XI(3))
dNdxi(4,3) =  MONE/EIGHT*(ONE-XI(1))*(ONE+XI(2))
dNdxi(5,1) =  MONE/EIGHT*(ONE-XI(2))*(ONE+XI(3))
dNdxi(5,2) =  MONE/EIGHT*(ONE-XI(1))*(ONE+XI(3))
dNdxi(5,3) =  ONE/EIGHT*(ONE-XI(1))*(ONE-XI(2))
dNdxi(6,1) =  ONE/EIGHT*(ONE-XI(2))*(ONE+XI(3))
dNdxi(6,2) =  MONE/EIGHT*(ONE+XI(1))*(ONE+XI(3))
dNdxi(6,3) =  ONE/EIGHT*(ONE+XI(1))*(ONE-XI(2))
dNdxi(7,1) =  ONE/EIGHT*(ONE+XI(2))*(ONE+XI(3))
dNdxi(7,2) =  ONE/EIGHT*(ONE+XI(1))*(ONE+XI(3))
dNdxi(7,3) =  ONE/EIGHT*(ONE+XI(1))*(ONE+XI(2))
dNdxi(8,1) =  MONE/EIGHT*(ONE+XI(2))*(ONE+XI(3))
dNdxi(8,2) =  ONE/EIGHT*(ONE-XI(1))*(ONE+XI(3))
dNdxi(8,3) =  ONE/EIGHT*(ONE-XI(1))*(ONE+XI(2))

dXdXi = zero
dXdEta = zero
dXdZeta = zero
dYdXi = zero
dYdEta = zero
dYdZeta = zero
dZdXi = zero
dZdEta = zero
dZdZeta = zero

do k=1,8
dXdXi = dXdXi + dNdxi(k,1)*coords(1,k)
dXdEta = dXdEta + dNdxi(k,2)*coords(1,k)
dXdZeta = dXdZeta + dNdxi(k,3)*coords(1,k)
dYdXi = dYdXi + dNdxi(k,1)*coords(2,k)
dYdEta = dYdEta + dNdxi(k,2)*coords(2,k)
dYdZeta = dYdZeta + dNdxi(k,3)*coords(2,k)
dZdXi = dZdXi + dNdxi(k,1)*coords(3,k)
dZdEta = dZdEta + dNdxi(k,2)*coords(3,k)
dZdZeta = dZdZeta + dNdxi(k,3)*coords(3,k)
enddo

DO face=1,6  !LOOP OVER ALL 6 FACES OF THE SMALL CUBE

! Jacobian of the mapping
!
if((face.eq.1).or.(face.eq.2)) then
! zeta = constant on this face
dA = dsqrt((dYdXi*dZdEta - dYdEta*dZdXi)**two+ (dXdXi*dZdEta -&
&dXdEta*dZdXi)**two+ (dXdXi*dYdEta - dXdEta*dYdXi)**two)
elseif((face.eq.3).or.(face.eq.5)) then
! eta = constant on this face
dA = dsqrt((dYdXi*dZdZeta - dYdZeta*dZdXi)**two+ (dXdXi*dZdZeta -&
&dXdZeta*dZdXi)**two+ (dXdXi*dYdZeta - dXdZeta*dYdXi)**two)
elseif((face.eq.4).or.(face.eq.6)) then
! xi = constant on this face
dA = dsqrt((dYdEta*dZdZeta - dYdZeta*dZdEta)**two+ (dXdEta*dZdZeta -&
dXdZeta*dZdEta)**two+ (dXdEta*dYdZeta - dXdZeta*dYdEta)**two)

endif

write(*,*) 'face',face,'area is dA', dA

if((face.eq.1).or.(face.eq.2)) then
! zeta = constant on this face
normal(1,1) = dYdXi*dZdEta - dYdEta*dZdXi
normal(2,1) = dXdXi*dZdEta - dXdEta*dZdXi
normal(3,1) = dXdXi*dYdEta - dXdEta*dYdXi
if(face.eq.1) normal = -normal
elseif((face.eq.3).or.(face.eq.4)) then
! eta = constant on this face
normal(1,1) = dYdXi*dZdZeta - dYdZeta*dZdXi
normal(2,1) = dXdXi*dZdZeta - dXdZeta*dZdXi
normal(3,1) = dXdXi*dYdZeta - dXdZeta*dYdXi
if(face.eq.3) normal = -normal
elseif((face.eq.5).or.(face.eq.6)) then
! xi = constant on this face
normal(1,1) = dYdEta*dZdZeta - dYdZeta*dZdEta
normal(2,1) = dXdEta*dZdZeta - dXdZeta*dZdEta
normal(3,1) = dXdEta*dYdZeta - dXdZeta*dYdEta
if(face.eq.6) normal = -normal

endif
mag = dsqrt(normal(1,1)**two+normal(2,1)**two+normal(3,1)**two)

normal(1,1) = normal(1,1)/mag
normal(2,1) = normal(2,1)/mag
normal(3,1) = normal(3,1)/mag

write(*,*) 'face', face, 'normal is ',normal

END DO


The area of each surface of each small cube should be .25, for this case but this code it isn't returning so. The normals are also not correctly oriented, here the cube is aligned with the global axes so it should be either 1 or -1.

Note: a. This might always not be a cube, it can be any irregular hexahedron, so the areas will not be equal always thus we need to compute each of them. b. The faces might be oriented along different directions, so the isoparametric transformation is sort of necessary.

Also, I have approched this using vectors of each face, without using the isoparametric transformation, where I use the diagonals of face of the cube and compute their cross product to find the area and outward unit normal, that approch also fails when the element is not regular. That part of the code looks like:

IF (J==1) THEN

xsd1(1:4)=xsd(1:4)
ysd1(1:4)=ysd(1:4)
zsd1(1:4)=zsd(1:4)
b(1)=xsd1(1)-xsd1(3)
b(2)=ysd1(1)-ysd1(3)
b(3)=zsd1(1)-zsd1(3)
c(1)=xsd1(2)-xsd1(4)
c(2)=ysd1(2)-ysd1(4)
c(3)=zsd1(2)-zsd1(4)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
ELSE IF (J==2) THEN

xsd2(1:4)=xsd(5:8)
ysd2(1:4)=ysd(5:8)
zsd2(1:4)=zsd(5:8)
b(1)=xsd2(4)-xsd2(2)
b(2)=ysd2(4)-ysd2(2)
b(3)=zsd2(4)-zsd2(2)
c(1)=xsd2(3)-xsd2(1)
c(2)=ysd2(3)-ysd2(1)
c(3)=zsd2(3)-zsd2(1)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
ELSE IF (J==3) THEN

xsd3(1:2)=xsd(1:2)
xsd3(3)=xsd(6)
xsd3(4)=xsd(5)
ysd3(1:2)=ysd(1:2)
ysd3(3)=ysd(6)
ysd3(4)=ysd(5)
zsd3(1:2)=zsd(1:2)
zsd3(3)=zsd(6)
zsd3(4)=zsd(5)
b(1)=xsd3(4)-xsd3(2)
b(2)=ysd3(4)-ysd3(2)
b(3)=zsd3(4)-zsd3(2)
c(1)=xsd3(3)-xsd3(1)
c(2)=ysd3(3)-ysd3(1)
c(3)=zsd3(3)-zsd3(1)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
ELSE IF (J==4) THEN

xsd4(1:2)=xsd(7:8)
xsd4(3)=xsd(4)
xsd4(4)=xsd(3)
ysd4(1:2)=ysd(7:8)
ysd4(3)=ysd(4)
ysd4(4)=ysd(3)
zsd4(1:2)=zsd(7:8)
zsd4(3)=zsd(4)
zsd4(4)=zsd(3)
b(1)=xsd4(1)-xsd4(3)
b(2)=ysd4(1)-ysd4(3)
b(3)=zsd4(1)-zsd4(3)
c(1)=xsd4(2)-xsd4(4)
c(2)=ysd4(2)-ysd4(4)
c(3)=zsd4(2)-zsd4(4)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
ELSE IF (J==5) THEN

xsd5(1:2)=xsd(2:3)
xsd5(3)=xsd(7)
xsd5(4)=xsd(6)
ysd5(1:2)=ysd(2:3)
ysd5(3)=ysd(7)
ysd5(4)=ysd(6)
zsd5(1:2)=zsd(2:3)
zsd5(3)=zsd(7)
zsd5(4)=zsd(6)
b(1)=xsd5(4)-xsd5(2)
b(2)=ysd5(4)-ysd5(2)
b(3)=zsd5(4)-zsd5(2)
c(1)=xsd5(3)-xsd5(1)
c(2)=ysd5(3)-ysd5(1)
c(3)=zsd5(3)-zsd5(1)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
ELSE IF (J==6) THEN

xsd6(1)=xsd(1)
xsd6(2)=xsd(4)
xsd6(3)=xsd(8)
xsd6(4)=xsd(5)
ysd6(1)=ysd(1)
ysd6(2)=ysd(4)
ysd6(3)=ysd(8)
ysd6(4)=ysd(5)
zsd6(1)=zsd(1)
zsd6(2)=zsd(4)
zsd6(3)=zsd(8)
zsd6(4)=zsd(5)
b(1)=xsd6(3)-xsd6(1)
b(2)=ysd6(3)-ysd6(1)
b(3)=zsd6(3)-zsd6(1)
c(1)=xsd6(4)-xsd6(2)
c(2)=ysd6(4)-ysd6(2)
c(3)=zsd6(4)-zsd6(2)
call cross(c, b,crossproduct)
norm(:)=crossproduct/(crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5
area_isd=((crossproduct(1)**2+crossproduct(2)**2+crossproduct(3)**2)**0.5)/2
END IF


J being the faces, I loop over the points in the respective faces and find cross product of the diagonals.

I would like to know about what is wrong in both these approches and the way to correct it. I would be glad if some one could give me a way forward.

• Are you asking where the bug in the code is? Commented Sep 13, 2018 at 15:46
• Not really. I just posted the code to make my algorithm clear. I want to if my theory/algorithm is correct. If not, what will be the best possible way to do this. Commented Sep 13, 2018 at 19:56

To debug situations like this, it's easiest if you break up the problem into several parts. Since you don't really need an isoparametric map to compute areas and normals, I'll ignore that aspect of your approach.

Step 1: Find the orientation of the nodes on the face. They should follow a right-hand rule (counterclockwise sequence). You can use the approach in https://en.wikipedia.org/wiki/Curve_orientation (with the appropriate projection on to the xy-, yz-, or zx-planes). For example

function isCounterClockwise(x_coords, y_coords, z_coords)
real   :: isCounterClockwise
real   :: x_coords(4)
real   :: y_coords(4)
real   :: z_coords(4)

isCounterClockwise = ....
return


Step 2: If the points are not counter-clockwise, reverse the order.

subroutine reverseOrder(x_coords, y_coords, z_coords)
....


Step 3: Compute cross product

subroutine computeAreaAndNormal(x_coords, y_coords, z_coords, area, normal)

u1 = x_coords(2) - x_coords(1)
u2 = y_coords(2) - y_coords(1)
u3 = z_coords(2) - z_coords(1)

v1 = x_coords(4) - x_coords(1)
v2 = y_coords(4) - y_coords(1)
v3 = z_coords(4) - z_coords(1)

normal = crossProduct(u1, u2, u3, v1, v2, v3)
area = sqrt(normal(1)**2 + normal(2)**2 + normal(3)**2)
normal = normal/area


If you write your code this way it will be much easier to debug. Caveat: Note that the above is not compileable Fortran.

Update:

For arbitrary planar quadrilaterals with nodes (1,2,3,4) arranged clockwise, find the areas of the triangles (1,2,3) and (1,3,4) (half of the norm of the cross product) and add them.

• Thanks for the detailed explanation.I was confident that the vector approach should work, but it didn't. My code has a similar structure as your suggestion, where the only difference was that I transform the given hexahedron into another hex, whose face and vertices are oriented/numbered as I want them to be(anti-clockwise in this case). Thus the anti-clockwise check and re-order subroutine are redundant.This works for a simple cube oriented along the global axes, but not for randomly oriented hexahedrons.I cannot find a flaw in this logic, as to why it shouldn't work for a general hex. Commented Sep 14, 2018 at 1:39
• I am assuming the the points of a surface remain on a plane, by that assumption the vector product of the sides (numbered in anti-clockwise fashion) should give the area of any arbitrary hex and also the outward unit normal. Commented Sep 14, 2018 at 1:42
• That indicates there's a bug in your program. Break up your code into small functions and test each function using unit tests. The bug should be easy to find. Commented Sep 14, 2018 at 2:13
• Also, your results suggest that there may be something wrong with your isoparametric map. Bugs are caused when one assumes that what they have coded is what they imagined. It's usually the trivial bugs that are the most difficult to find (hence the orientation check). Try the direct method first, before the isoparametric method. Commented Sep 14, 2018 at 2:24
• Thanks for the tips. Just to mention, the results I mentioned in the comments are computed directly using the vector product, no isoparametric map is involved. I will definitely be very glad if I can solve this without using the isoparametric transformation, which would unnecessarily make the code clumsier. Commented Sep 14, 2018 at 3:00