Calculate surface area and normal for each face of an arbitrary hexahedron

Question

I am trying to find out the surface area of each of the faces of a cube and the corresponding outward unit normals. This operation is done on a finite element mesh, so I have transformed each surface of the cube into the isoparametric form using shape(basis) functions and then tried to extract the area and normals.

Here is the code:

program polyhedron

IMPLICIT NONE

real(8)  coord(3,8)
INTEGER  face, INPT, I, ino,k, ii  
REAL(8)  XI(3), dNdxi(8,3), ZERO, ONE, MONE, EIGHT
REAL(8) VJACOB(3,3),XII(8,3),norm(3), TWO, THREE, FOUR, HALF, SIX
PARAMETER(ZERO=0.D0,ONE=1.D0,MONE=-1.D0,EIGHT=8.D0)
PARAMETER(TWO=2.D0,THREE=3.D0,FOUR=4.D0,HALF=0.5D0,SIX=6.D0)
REAL(8) dXdXi,dXdEta,dXdZeta,dYdXi,dYdEta,dYdZeta,dZdXi,dZdEta,area_isd
REAL(8) dZdZeta, dA, mag, normal(3,1),xLocal(4),yLocal(4),zLocal(4)


!COORDINATES OF THE CUBE

 coord(1,1)=1.00
 coord(2,1)=1.00
 coord(3,1)=1.00
 coord(1,2)=1.00
 coord(2,2)=0.00
 coord(3,2)=1.00
 coord(1,3)=1.00
 coord(2,3)=1.00
 coord(3,3)=0.00
 coord(1,4)=1.00
 coord(2,4)=0.00
 coord(3,4)=0.00
 coord(1,5)=0.00
 coord(2,5)=1.00
 coord(3,5)=1.00
 coord(1,6)=0.00
 coord(2,6)=0.00
 coord(3,6)=1.00
 coord(1,7)=0.00
 coord(2,7)=1.00
 coord(3,7)=0.00
 coord(1,8)=0.00
 coord(2,8)=0.00
 coord(3,8)=0.00

do face=1,6   !Loop over the faces
area_isd=0.0
call xintSurf3D4pt(face,xLocal,yLocal,zLocal) !get local points
do ii=1,4
call computeSurf3D(xLocal(ii),yLocal(ii),zLocal(ii),face,coord,dA,norm) !compute area and normal
area_isd=area_isd+dA
end do
write(*,*) 'face', face, 'area', area_isd
write(*,*) 'norm', norm
end do
end program polyhedron

The subroutine to calculate the local jacobian and the normals are:

subroutine computeSurf3D(xLocal,yLocal,zLocal,face,coords,dA,normal)



IMPLICIT NONE

integer face,stat,i,j,k

real(8) xLocal,yLocal,zLocal,dA,dshxi(8,3),zero,dsh(8,3),one
real(8) coords(3,8),two,eighth,mapJ(3,3),mag,normal(3)

real(8) dXdXi,dXdEta,dXdZeta,dYdXi,dYdEta,dYdZeta,dZdXi,dZdEta
real(8) dZdZeta

parameter(one=1.d0,two=2.d0,eighth=1.d0/8.d0,zero=0.d0)

!Hex shape function derivatives  
dshxi(1,1) = -eighth*(one - yLocal)*(one - zLocal)
dshxi(1,2) = -eighth*(one - xLocal)*(one - zLocal)
dshxi(1,3) = -eighth*(one - xLocal)*(one - yLocal)
dshxi(2,1) = eighth*(one - yLocal)*(one - zLocal)
dshxi(2,2) = -eighth*(one + xLocal)*(one - zLocal)
dshxi(2,3) = -eighth*(one + xLocal)*(one - yLocal)
dshxi(3,1) = eighth*(one + yLocal)*(one - zLocal)
dshxi(3,2) = eighth*(one + xLocal)*(one - zLocal)
dshxi(3,3) = -eighth*(one + xLocal)*(one + yLocal)
dshxi(4,1) = -eighth*(one + yLocal)*(one - zLocal)
dshxi(4,2) = eighth*(one - xLocal)*(one - zLocal)
dshxi(4,3) = -eighth*(one - xLocal)*(one + yLocal)
dshxi(5,1) = -eighth*(one - yLocal)*(one + zLocal)
dshxi(5,2) = -eighth*(one - xLocal)*(one + zLocal)
dshxi(5,3) = eighth*(one - xLocal)*(one - yLocal)
dshxi(6,1) = eighth*(one - yLocal)*(one + zLocal)
dshxi(6,2) = -eighth*(one + xLocal)*(one + zLocal)
dshxi(6,3) = eighth*(one + xLocal)*(one - yLocal)
dshxi(7,1) = eighth*(one + yLocal)*(one + zLocal)
dshxi(7,2) = eighth*(one + xLocal)*(one + zLocal)
dshxi(7,3) = eighth*(one + xLocal)*(one + yLocal)
dshxi(8,1) = -eighth*(one + yLocal)*(one + zLocal)
dshxi(8,2) = eighth*(one - xLocal)*(one + zLocal)
dshxi(8,3) = eighth*(one - xLocal)*(one + yLocal)


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


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



      !
      if((face.eq.1).or.(face.eq.2)) then
         ! zeta = constant on this face
         normal(1) = dYdXi*dZdEta - dYdEta*dZdXi
         normal(2) = dXdXi*dZdEta - dXdEta*dZdXi
         normal(3) = 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) = dYdXi*dZdZeta - dYdZeta*dZdXi
         normal(2) = dXdXi*dZdZeta - dXdZeta*dZdXi
         normal(3) = 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) = dYdEta*dZdZeta - dYdZeta*dZdEta
         normal(2) = dXdEta*dZdZeta - dXdZeta*dZdEta
         normal(3) = dXdEta*dYdZeta - dXdZeta*dYdEta
         if(face.eq.5) normal = -normal
      endif
      mag = dsqrt(normal(1)**two+normal(2)**two+normal(3)**two)
      normal(1) = normal(1)/mag
      normal(2) = normal(2)/mag
      normal(3) = normal(3)/mag


end subroutine computeSurf3D

The local Gauss points are obtained from this subroutine:

subroutine xintSurf3D4pt(face,xLocal,yLocal,zLocal)



      implicit none

integer face
real(8) xLocal(4),yLocal(4),zLocal(4),w(4),one,three
parameter(one=1.d0,three=3.d0)




      ! Gauss pt locations in master element
      !
      if(face.eq.1) then
         xLocal(1) = -dsqrt(one/three)
         yLocal(1) = -dsqrt(one/three)
         zLocal(1) = -one
         xLocal(2) = dsqrt(one/three)
         yLocal(2) = -dsqrt(one/three)
         zLocal(2) = -one
         xLocal(3) = dsqrt(one/three)
         yLocal(3) = dsqrt(one/three)
         zLocal(3) = -one
         xLocal(4) = -dsqrt(one/three)
         yLocal(4) = dsqrt(one/three)
         zLocal(4) = -one
      elseif(face.eq.2) then
         xLocal(1) = -dsqrt(one/three)
         yLocal(1) = -dsqrt(one/three)
         zLocal(1) = one
         xLocal(2) = dsqrt(one/three)
         yLocal(2) = -dsqrt(one/three)
         zLocal(2) = one
         xLocal(3) = dsqrt(one/three)
         yLocal(3) = dsqrt(one/three)
         zLocal(3) = one
         xLocal(4) = -dsqrt(one/three)
         yLocal(4) = dsqrt(one/three)
         zLocal(4) = one
      elseif(face.eq.3) then
         xLocal(1) = -dsqrt(one/three)
         yLocal(1) = -one
         zLocal(1) = -dsqrt(one/three)
         xLocal(2) = dsqrt(one/three)
         yLocal(2) = -one
         zLocal(2) = -dsqrt(one/three)
         xLocal(3) = dsqrt(one/three)
         yLocal(3) = -one
         zLocal(3) = dsqrt(one/three)
         xLocal(4) = -dsqrt(one/three)
         yLocal(4) = -one
         zLocal(4) = dsqrt(one/three)
      elseif(face.eq.4) then
         xLocal(1) = one
         yLocal(1) = -dsqrt(one/three)
         zLocal(1) = -dsqrt(one/three)
         xLocal(2) = one
         yLocal(2) = dsqrt(one/three)
         zLocal(2) = -dsqrt(one/three)
         xLocal(3) = one
         yLocal(3) = dsqrt(one/three)
         zLocal(3) = dsqrt(one/three)
         xLocal(4) = one
         yLocal(4) = -dsqrt(one/three)
         zLocal(4) = dsqrt(one/three)
      elseif(face.eq.5) then
         xLocal(1) = -dsqrt(one/three)
         yLocal(1) = one
         zLocal(1) = -dsqrt(one/three)
         xLocal(2) = dsqrt(one/three)
         yLocal(2) = one
         zLocal(2) = -dsqrt(one/three)
         xLocal(3) = dsqrt(one/three)
         yLocal(3) = one
         zLocal(3) = dsqrt(one/three)
         xLocal(4) = -dsqrt(one/three)
         yLocal(4) = one
         zLocal(4) = dsqrt(one/three)
      elseif(face.eq.6) then
         xLocal(1) = -one
         yLocal(1) = -dsqrt(one/three)
         zLocal(1) = -dsqrt(one/three)
         xLocal(2) = -one
         yLocal(2) = dsqrt(one/three)
         zLocal(2) = -dsqrt(one/three)
         xLocal(3) = -one
         yLocal(3) = dsqrt(one/three)
         zLocal(3) = dsqrt(one/three)
         xLocal(4) = -one
         yLocal(4) = -dsqrt(one/three)
         zLocal(4) = dsqrt(one/three)
      endif

      end subroutine xintSurf3D4pt

The area of each surface should be 1 unit, for this case but this code it isn't returning so. The normals are are also incorrect. The output is

    face           1 area  0.57735026918962573     
 norm   1.0000000000000000        0.0000000000000000E+000   0.0000000000000000E+000
 face           2 area  0.57735026918962573     
 norm  -1.0000000000000000       -0.0000000000000000E+000  -0.0000000000000000E+000
 face           3 area   1.0000000000000000     
 norm   0.0000000000000000E+000  -0.0000000000000000E+000   1.0000000000000000     
 face           4 area  0.57735026918962573     
 norm  -0.0000000000000000E+000   0.0000000000000000E+000  -1.0000000000000000     
 face           5 area   1.1547005383792515     
 norm  -0.0000000000000000E+000  0.86602540378443871       0.50000000000000000     
 face           6 area   1.4142135623730951     
 norm   0.0000000000000000E+000 -0.70710678118654746      -0.70710678118654746

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.

Is this the correct way to approach this issue? Will be glad if someone can help me figure this out. I have also tried using vector product of diagonals for the each of the faces to calculate area and unit normal, but they don't work when the structure is irregular. Here is an example picture of an irregular hexahedron 1. A regular rube roughly cube looks something like this. 2

Could you please provide a [mcve], limited just to the part related to this question? Also, show us the output you get (any error or result, even if it is wrong) — Rodrigo Rodrigues, Sep 13 '18 at 00:34
hi, I made it as short and complete as possible. This code calculates the surface area and normals for six faces of a regular cube. — Schneider, Sep 13 '18 at 01:13
Thanks for the comments. Not setting a value to the for loop was indeed a mistake and I corrected it, still can't get final result. I have updated the code and the output. The code is a bit longer now but probably explains the problem better. I also added `implicit none`. — Schneider, Sep 13 '18 at 19:50
This looks like a geometry problem and not a programming problem, now. — Ross, Sep 13 '18 at 20:37
Actually it would be great to see some figure/plot of the shape you are thinking about in this problem. — msi_gerva, Sep 13 '18 at 20:40
Could you kindly explain what do you think about the geometry problem? It is just a simple cubic structure @Ross — Schneider, Sep 13 '18 at 21:08
@msi_gerva : For this example I am just using a unit cube. It can be any tetrahedron, I will try to post a picture. I am using the co-ordinate transfomration to handle arbitrary shapes. — Schneider, Sep 13 '18 at 21:16
I would also want to add that if I change the sequence of the coordinates of the 8 points, the result changes. Which shouldn't be the case, given 8 points the area and outward normal of a hexahedron should be constant. This leads to believe that the code is wrong. I just can't figure out where. — Schneider, Sep 13 '18 at 22:40
Given the answer you accepted I think it is now clear to you that it indeed is a geometry problem and not a programming problem. — Vladimir F Героям слава, Sep 14 '18 at 11:33
I still don't see it though. Because this is just another approach to calculate area, whose theory is correct. So this should also give me the correct result unless there is a programming bug. I accepted the vector product answer because that probably will make the code more compact. — Schneider, Sep 14 '18 at 12:58

score 2 · Accepted Answer · answered Sep 14 '18 at 08:10

2

Normal: take three vertices and compute the cross product of the vectors they form.

Area: apply the shoelace formula on XY, YZ, and ZX, then take the Euclidean norm of the three results.

answered Sep 14 '18 at 08:10

Just to confirm, this would work for any surface of an arbitrary hexahedron, right? My knowledge of vector says so, but still, want to confirm. Because I tried this approach and it is not working, so if the logic is right then I can be sure that there is a bug in the code and not the logic. – Schneider Sep 14 '18 at 12:28
@Sauradeep: the faces must be oriented consitently. – Sep 14 '18 at 13:51
Yeah, right. An anticlockwise or right-handed orientation of the vertices in a particular face is necessary. – Schneider Sep 14 '18 at 14:13

Calculate surface area and normal for each face of an arbitrary hexahedron

1 Answers1