<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Apr 10, 2014 at 7:57 PM, Adrian Croucher <span dir="ltr"><<a href="mailto:a.croucher@auckland.ac.nz" target="_blank">a.croucher@auckland.ac.nz</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">hi<br>
<br>
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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
Any clues, or could this be a bug?<br></blockquote><div><br></div><div>Thanks for the test. It was a bug in a special case:</div><div><br></div><div><div>-        coords[3] = PetscRealPart(coords[3*dim+0] - coords[0*dim+0]);</div>
<div>+        coords[3] = PetscRealPart(coords[3*dim+1] - coords[0*dim+1]);</div></div><div><br></div><div>I will push it to next</div><div><br></div><div>  Thanks,</div><div><br></div><div>      Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">

Cheers, Adrian<br>
--<br>
<br>
program main<br>
<br>
  ! test using DMPlexComputeCellGeometryFVM() to compute face areas,<br>
  ! on simple mesh of two cubes<br>
<br>
  implicit none<br>
<br>
#include <finclude/petsc.h90><br>
<br>
  DM :: dm, dmi<br>
  PetscInt, parameter :: dim = 3<br>
  PetscInt :: depth = 1<br>
  PetscErrorCode :: ierr<br>
  PetscInt, allocatable :: numPoints(:)<br>
  PetscInt, allocatable :: coneSize(:)<br>
  PetscInt, allocatable :: cones(:)<br>
  PetscInt, allocatable :: coneOrientations(:)<br>
  PetscScalar, allocatable :: vertexCoords(:)<br>
  PetscReal ::  vol = 0.<br>
  PetscReal, target, dimension(dim)  :: centroid = 0., normal = 0.<br>
  PetscReal, pointer :: pcentroid(:), pnormal(:)<br>
  PetscInt :: i, fStart, fEnd, f<br>
  PetscScalar :: facearea<br>
<br>
  numPoints = [12, 2]<br>
  coneSize  = [8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]<br>
  cones = [2,5,4,3,6,7,8,9,  3,4,11,10,7,12,13,8]<br>
  coneOrientations = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]<br>
  vertexCoords = [ &<br>
       -1.0,0.0,0.0, 0.0,0.0,0.0, 0.0,1.0,0.0, -1.0,1.0,0.0, &<br>
       -1.0,0.0,1.0, 0.0,0.0,1.0, 0.0,1.0,1.0, -1.0,1.0,1.0, &<br>
        1.0,0.0,0.0, 1.0,1.0,0.0, 1.0,0.0,1.0,  1.0,1.0,1.0 ]<br>
<br>
  call PetscInitialize(PETSC_NULL_<u></u>CHARACTER,ierr); CHKERRQ(ierr)<br>
<br>
  call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr); CHKERRQ(ierr)<br>
  call DMPlexSetDimension(dm, dim, ierr); CHKERRQ(ierr)<br>
<br>
  call DMPlexCreateFromDAG(dm, depth, numPoints, coneSize, cones, &<br>
       coneOrientations, vertexCoords, ierr); CHKERRQ(ierr)<br>
<br>
  dmi = PETSC_NULL_OBJECT<br>
  call DMPlexInterpolate(dm, dmi, ierr); CHKERRQ(ierr)<br>
  call DMPlexCopyCoordinates(dm, dmi, ierr); CHKERRQ(ierr)<br>
  call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
  dm = dmi<br>
<br>
  pcentroid => centroid<br>
  pnormal => normal<br>
  call DMPlexGetHeightStratum(dm, 1, fStart, fEnd, ierr); CHKERRQ(ierr)<br>
  write(*,'(a2, 1x, 2a10)'), ' f', 'area', 'centroid'<br>
  do f = fStart, fEnd-1<br>
     call DMPlexComputeCellGeometryFVM(<u></u>dm, f, facearea, pcentroid, pnormal, ierr)<br>
     write(*,'(i2, 1x, 4f10.4)'), f, facearea, centroid<br>
  end do<br>
<br>
  call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
  call PetscFinalize(ierr); CHKERRQ(ierr)<br>
<br>
end program main<span class=""><font color="#888888"><br>
<br>
-- <br>
Dr Adrian Croucher<br>
Department of Engineering Science<br>
University of Auckland<br>
New Zealand<br>
tel 64-9-373-7599 ext 84611<br>
<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div>