0
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ Are you asking where the bug in the code is? $\endgroup$
    – Wolfgang Bangerth
    Commented Sep 13, 2018 at 15:46
  • $\begingroup$ 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. $\endgroup$
    – Schneider
    Commented Sep 13, 2018 at 19:56

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Schneider
    Commented Sep 14, 2018 at 1:39
  • $\begingroup$ 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. $\endgroup$
    – Schneider
    Commented Sep 14, 2018 at 1:42
  • $\begingroup$ 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. $\endgroup$
    – Biswajit Banerjee
    Commented Sep 14, 2018 at 2:13
  • $\begingroup$ 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. $\endgroup$
    – Biswajit Banerjee
    Commented Sep 14, 2018 at 2:24
  • $\begingroup$ 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. $\endgroup$
    – Schneider
    Commented Sep 14, 2018 at 3:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.