On Thu, Apr 10, 2014 at 7:57 PM, Adrian Croucher
<a.crouc...@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
