[petsc-dev] DMPlexComputeCellGeometryFVM from fortran

[Adrian Croucher]

hi,
> They need to be F90 pointers and they need to be associated with an
> allocated array. Take
> a look at how I do it for the arguments to DMPlexCreateSection() in DMPlex
> ex1f90.F
Thanks, I tried that and it nearly works... it gives the correct 
volumes, but it gives centroids all zero, which isn't right. My amended 
code is below.

If I debug into DMPlexComputeGeometryFVM_3D_Internal() using gdb, and 
try 'p *centroid', it says 'Cannot access memory at address 0x0', which 
looks unhealthy.

There is something a bit weird going on, because if I don't initialize 
the centroid and normal arrays, or even if I initialize them in the main 
code rather than in the variable declaration, I get a segfault in 
DMPlexComputeGeometryFVM_3D_Internal().

Guess I'm still doing something subtly wrong here?

- Adrian

--
program main

! setting up 3-D DMPlex using DMPlexCreateFromDAG(), DMPlexInterpolate() and
! DMPlexComputeCellGeometryFVM()

   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

   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 = [ &
        -0.5,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -0.5,1.0,0.0, &
        -0.5,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -0.5,1.0,1.0, &
         0.5,0.0,0.0, 0.5,1.0,0.0, 0.5,0.0,1.0,  0.5,1.0,1.0 ]
   pcentroid => centroid
   pnormal => normal

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

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

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

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

   call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)

   do i = 0, 1
      call DMPlexComputeCellGeometryFVM(dm, i, vol, pcentroid, pnormal, 
ierr)
      CHKERRQ(ierr)
      write(*, '(a, i2, a, f8.4, a, 3(f8.4, 1x))'), &
           'cell: ', i, ' volume: ', vol, ' centroid: ', &
           pcentroid(1), pcentroid(2), pcentroid(3)
   end do

   call DMDestroy(dm, ierr); CHKERRQ(ierr)
   call PetscFinalize(ierr); CHKERRQ(ierr)
   deallocate(numPoints, coneSize, cones, coneOrientations, vertexCoords)

end program main

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