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

[Matthew Knepley]

On Thu, Apr 10, 2014 at 7:57 PM, Adrian Croucher
<a.croucher at auckland.ac.nz>wrote:

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

Thanks for the test. It was a bug in a special case:

-        coords[3] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);
+        coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]);

I will push it to next

  Thanks,

      Matt


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


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20140412/95405449/attachment.html>