[petsc-dev] mesh face areas with DMPlexComputeCellGeometryFVM()

[Adrian Croucher]

hi

I am getting some slightly weird results from 
DMPlexComputeCellGeometryFVM() when I use it to compute mesh face areas 
and centroids. The Fortran code below demonstrates the problem.

It constructs a very simple mesh with 2 cube cells, length 1.0 along 
each side. So the face areas should all be 1.0.

Most of them are, but the two bottom faces (z = 0) are returning areas 
of 0.5, and the wrong centroids. It appears they're being computed based 
on triangles with one of the face corners left out.

Any clues, or could this be a bug?

Cheers, Adrian
--

program main

   ! test using DMPlexComputeCellGeometryFVM() to compute face areas,
   ! on simple mesh of two cubes

   implicit none

#include <finclude/petsc.h90>

   DM :: dm, dmi
   PetscInt, parameter :: dim = 3
   PetscInt :: depth = 1
   PetscErrorCode :: ierr
   PetscInt, allocatable :: numPoints(:)
   PetscInt, allocatable :: coneSize(:)
   PetscInt, allocatable :: cones(:)
   PetscInt, allocatable :: coneOrientations(:)
   PetscScalar, allocatable :: vertexCoords(:)
   PetscReal ::  vol = 0.
   PetscReal, target, dimension(dim)  :: centroid = 0., normal = 0.
   PetscReal, pointer :: pcentroid(:), pnormal(:)
   PetscInt :: i, fStart, fEnd, f
   PetscScalar :: facearea

   numPoints = [12, 2]
   coneSize  = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
   cones = [2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8]
   coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
   vertexCoords = [ &
        -1.0,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -1.0,1.0,0.0, &
        -1.0,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -1.0,1.0,1.0, &
         1.0,0.0,0.0, 1.0,1.0,0.0, 1.0,0.0,1.0,  1.0,1.0,1.0 ]

   call PetscInitialize(PETSC_NULL_CHARACTER,ierr); CHKERRQ(ierr)

   call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr); CHKERRQ(ierr)
   call DMPlexSetDimension(dm, dim, ierr); CHKERRQ(ierr)

   call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, &
        coneOrientations, vertexCoords, ierr); CHKERRQ(ierr)

   dmi = PETSC_NULL_OBJECT
   call DMPlexInterpolate(dm, dmi, ierr); CHKERRQ(ierr)
   call DMPlexCopyCoordinates(dm, dmi, ierr); CHKERRQ(ierr)
   call DMDestroy(dm, ierr); CHKERRQ(ierr)
   dm = dmi

   pcentroid => centroid
   pnormal => normal
   call DMPlexGetHeightStratum(dm, 1, fStart, fEnd, ierr); CHKERRQ(ierr)
   write(*,'(a2, 1x, 2a10)'), ' f', 'area', 'centroid'
   do f = fStart, fEnd-1
      call DMPlexComputeCellGeometryFVM(dm, f, facearea, pcentroid, 
pnormal, ierr)
      write(*,'(i2, 1x, 4f10.4)'), f, facearea, centroid
   end do

   call DMDestroy(dm, ierr); CHKERRQ(ierr)
   call PetscFinalize(ierr); CHKERRQ(ierr)

end program main

-- 
Dr Adrian Croucher
Department of Engineering Science
University of Auckland
New Zealand
tel 64-9-373-7599 ext 84611