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
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.