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
