<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sun, Apr 6, 2014 at 9:43 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
  
    
  
  <div text="#000000" bgcolor="#FFFFFF">
    <br>
    <div>On 05/04/14 02:31, Matthew Knepley
      wrote:<br>
    </div>
    <blockquote type="cite">
      
      <div dir="ltr">
        <div class="gmail_extra">
          <div class="gmail_quote"><br>
            <div>Instead of this, use</div>
            <div><br>
            </div>
            <div>  call DMPlexGetDefaultSection(cell_dm, section, ierr);
              (Move outside of loop)</div>
            <div>  call PetscSectionGetOffset(section, c, off, ierr);</div>
            <div>  part(off) =rank</div>
          </div>
        </div>
      </div>
    </blockquote>
    <br>
    OK thanks. The code compiles now.<br>
    <br>
    However, at runtime I'm getting another error, when it calls
    DMSetPointSF():<br>
    <br>
    [0]PETSC ERROR: #1 DMSetPointSF() line 3204 in
    /home/acro018/software/PETSc/code/src/dm/interface/dm.c<br>
    [1]PETSC ERROR: --------------------- Error Message
    --------------------------------------------------------------<br>
    [1]PETSC ERROR: Invalid argument<br>
    [1]PETSC ERROR: Wrong type of object: Parameter # 1<br>
    [1]PETSC ERROR: See
    <a href="http://http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for
    trouble shooting.<br>
    [1]PETSC ERROR: Petsc Development GIT revision:
    v3.4.4-5306-g9849faf  GIT Date: 2014-04-06 15:35:41 -0500<br>
    <br>
    This seems to work fine in C, so I'm wondering if I need to do
    something different in fortran, or if the interface isn't quite
    working right. My test code is below. The offending call to
    DMSetPointSF() is in the routine create_partition_vector(). Any tips
    much appreciated!<br></div></blockquote><div><br></div><div>1) You do not need those calls since DMClone() now does this automatically.</div><div><br></div><div>2) However, it looks like DMClone() is not functioning correctly for you since cell_dm should be alright.</div>
<div><br></div><div>I will try and run your code today.</div><div><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div text="#000000" bgcolor="#FFFFFF">
    It also occurred to me there might just be something wrong with my
    mesh, so I also tested it using a mesh created with
    DMPLexCreateBoxMesh()- but this gave the same error (and also
    indicated that the fortran interface to DMPlexSetRefinementLimit()
    is probably also missing in action!).<br></div></blockquote><div><br></div><div>Fixed. Will push it today.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
    - Adrian<br>
    <br>
    --<br>
    <br>
    program testplex<br>
    <br>
      ! Testing PETSc dmplex from fortran, using a line of blocks in
    3-D.<br>
      ! Create the DMPlex by passing in the cell definitions to
    DMPlexCreateFromDAG(),<br>
      ! interpolate the result to form the faces and edges, distribute
    over the<br>
      ! processors, and visualize the partitions in VTK.<br>
    <br>
      implicit none<br>
    <br>
    #include <finclude/petsc.h90><br>
    <br>
      PetscInt, parameter  :: num_cells = 20<br>
      PetscReal, parameter :: cell_size = 10.<br>
      PetscErrorCode :: ierr<br>
      DM :: dm, dist_dm, cell_dm<br>
      character(5) :: partitioner = "chaco"<br>
      PetscInt :: overlap = 0<br>
      Vec :: partition<br>
    <br>
      call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)<br>
    <br>
      call create_1d_grid(num_cells, cell_size, dm, ierr); CHKERRQ(ierr)<br>
    <br>
      call DMPlexDistribute(dm, partitioner, overlap, PETSC_NULL_OBJECT,
    dist_dm, ierr)<br>
      CHKERRQ(ierr)<br>
      call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
      dm = dist_dm<br>
    <br>
      call DMView(dm, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr)<br>
    <br>
      call create_partition_vector(dm, cell_dm, partition, ierr);
    CHKERRQ(ierr)<br>
    <br>
      ! Visualize partition vector here...<br>
    <br>
      call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
      call PetscFinalize(ierr); CHKERRQ(ierr)<br>
    <br>
    end program testplex<br>
    <br>
!------------------------------------------------------------------------<br>
    <br>
    subroutine create_1d_grid(num_cells, cell_size, dm, ierr)<br>
    <br>
      ! Creates grid consisting of a 1D line of cube cells in 3D.<br>
    <br>
      implicit none<br>
    <br>
    #include <finclude/petsc.h90><br>
    <br>
      PetscInt, intent(in) :: num_cells<br>
      PetscReal, intent(in):: cell_size<br>
      DM, intent(out) :: dm<br>
      PetscErrorCode, intent(out) :: ierr<br>
      ! Locals:<br>
      PetscInt, parameter :: dim = 3<br>
      PetscMPIInt :: rank<br>
      PetscInt  :: num_vertices_per_face, num_vertices_per_cell<br>
      PetscInt  :: num_vertices, total_num_points<br>
      PetscInt  :: i, i0, i1, l<br>
      PetscInt  :: depth<br>
      real(8)   :: x0, dx<br>
      PetscInt, allocatable :: num_points(:)<br>
      PetscInt, allocatable :: cone_size(:)<br>
      PetscInt, allocatable :: cones(:)<br>
      PetscInt, allocatable :: cone_orientations(:)<br>
      PetscScalar, allocatable :: vertex_coords(:)<br>
      DM :: interp_dm<br>
    <br>
      ierr = 0<br>
      call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr); CHKERRQ(ierr)<br>
    <br>
      call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr); CHKERRQ(ierr)<br>
      call PetscObjectSetName(dm, "testplex", ierr); CHKERRQ(ierr)<br>
      call DMPlexSetDimension(dm, dim, ierr); CHKERRQ(ierr)<br>
    <br>
      if (rank == 0) then<br>
    <br>
         num_vertices_per_face = 2 ** (dim-1)<br>
         num_vertices_per_cell = 2 * num_vertices_per_face<br>
         num_vertices = num_vertices_per_face * (num_cells + 1)<br>
         l = num_vertices_per_face * dim<br>
         total_num_points = num_vertices + num_cells ! # points in DAG<br>
    <br>
         allocate(vertex_coords(0:num_vertices*dim-1),
    cones(0:num_cells*num_vertices_per_cell-1))<br>
    <br>
         ! Set up vertex coordinates<br>
         dx = real(cell_size)<br>
    <br>
         do i = 0, num_cells<br>
            i0 = i * l<br>
            i1 = (i+1)*l<br>
            x0 = real(i * cell_size)<br>
            vertex_coords(i0:i1) = \<br>
            [x0, 0.d0, 0.d0,\<br>
            x0, dx, 0.d0,\<br>
            x0, 0.d0, dx,\<br>
            x0, dx, dx]<br>
         end do<br>
    <br>
         ! Set up cells<br>
         do i = 0, num_cells-1<br>
            i0 = i * num_vertices_per_cell<br>
            i1 = (i+1) * num_vertices_per_cell<br>
            cones(i0:i1) = [0, 1, 5, 4, 2, 6, 7, 3] + i *
    num_vertices_per_face <br>
         end do<br>
         cones = cones + num_cells<br>
    <br>
         num_points = [num_vertices, num_cells]<br>
         allocate(cone_size(0: total_num_points-1))<br>
        
    allocate(cone_orientations(0:num_cells*num_vertices_per_cell-1))<br>
         cone_size = 0<br>
         cone_size(0: num_cells-1) = num_vertices_per_cell<br>
         cone_orientations = 0<br>
         depth = 1<br>
    <br>
         call DMPlexCreateFromDAG(dm, depth, num_points, cone_size,
    cones, &<br>
              cone_orientations, vertex_coords, ierr); CHKERRQ(ierr)<br>
    <br>
      else ! rank > 0: create empty DMPlex<br>
    <br>
         num_points = [0, 0, 0]; cone_size = [0]<br>
         cones = [0]; cone_orientations = [0]<br>
         vertex_coords = [0]<br>
         call DMPlexCreateFromDAG(dm, 1, num_points, cone_size, &<br>
              cones, cone_orientations, vertex_coords, ierr);
    CHKERRQ(ierr)<br>
    <br>
      end if<br>
    <br>
      interp_dm = PETSC_NULL_OBJECT<br>
      call DMPlexInterpolate(dm, interp_dm, ierr); CHKERRQ(ierr)<br>
      call DMPlexCopyCoordinates(dm, interp_dm, ierr); CHKERRQ(ierr)<br>
      call DMDestroy(dm, ierr); CHKERRQ(ierr)<br>
      dm = interp_dm<br>
    <br>
      deallocate(num_points, cone_size, cones, cone_orientations,
    vertex_coords)<br>
    <br>
    end subroutine create_1d_grid<br>
    <br>
!------------------------------------------------------------------------<br>
    <br>
    subroutine create_partition_vector(dm, cell_dm, partition, ierr)<br>
    <br>
      ! Returns a DM and vector containing partitioning data<br>
    <br>
      implicit none<br>
    <br>
    #include <finclude/petsc.h90><br>
    <br>
      DM, intent(in) :: dm<br>
      DM, intent(out) :: cell_dm<br>
      Vec, intent(out) :: partition<br>
      PetscErrorCode, intent(out) :: ierr<br>
      ! Locals:<br>
      MPI_Comm :: comm<br>
      PetscMPIInt :: rank<br>
      PetscSF :: sfPoint<br>
      PetscSection :: coord_section<br>
      Vec :: coordinates<br>
      PetscSection :: cell_section<br>
      PetscScalar, pointer :: part(:)<br>
      PetscInt :: cStart, cEnd, c, off<br>
      <br>
      call DMGetCoordinateSection(dm, coord_section, ierr);
    CHKERRQ(ierr)<br>
      call DMGetCoordinatesLocal(dm, coordinates, ierr); CHKERRQ(ierr)<br>
      call DMClone(dm, cell_dm, ierr); CHKERRQ(ierr)<br>
      call DMGetPointSF(dm, sfPoint, ierr); CHKERRQ(ierr)<br>
      call DMSetPointSF(cell_dm, sfPoint, ierr); CHKERRQ(ierr)<br>
      call DMSetCoordinateSection(cell_dm, coord_section, ierr);
    CHKERRQ(ierr)<br>
      call DMSetCoordinatesLocal(cell_dm, coordinates, ierr);
    CHKERRQ(ierr)<br>
      call PetscObjectGetComm(dm, comm, ierr); CHKERRQ(ierr)<br>
      call PetscSectionCreate(comm, cell_section, ierr); CHKERRQ(ierr)<br>
      call DMPlexGetHeightStratum(cell_dm, 0, cStart, cEnd, ierr);
    CHKERRQ(ierr)<br>
      call PetscSectionSetChart(cell_section, cStart, cEnd, ierr);
    CHKERRQ(ierr)<br>
      do c = cStart, cEnd-1<br>
         call PetscSectionSetDof(cell_section, c, 1, ierr);
    CHKERRQ(ierr)<br>
      end do<br>
      call PetscSectionSetUp(cell_section, ierr); CHKERRQ(ierr)<br>
      call DMSetDefaultSection(cell_dm, cell_section, ierr);
    CHKERRQ(ierr)<br>
      call PetscSectionDestroy(cell_section, ierr); CHKERRQ(ierr)<br>
      call DMCreateLocalVector(cell_dm, partition, ierr); CHKERRQ(ierr)<br>
      call PetscObjectSetName(partition, "partition", ierr);
    CHKERRQ(ierr)<br>
      call VecGetArrayF90(partition, part, ierr); CHKERRQ(ierr)<br>
      call MPI_Comm_rank(comm, rank); CHKERRQ(ierr)<br>
      do c = cStart, cEnd-1<br>
         call PetscSectionGetOffset(cell_section, c, off, ierr);
    CHKERRQ(ierr)<br>
         part(off) = rank<br>
      end do<br>
      call VecRestoreArrayF90(partition, part, ierr); CHKERRQ(ierr)<br>
    <br>
    end subroutine create_partition_vector<span class="HOEnZb"><font color="#888888"><br>
    <br>
    <br>
    <br>
    <pre cols="72">-- 
Dr Adrian Croucher
Department of Engineering Science
University of Auckland
New Zealand
tel 64-9-373-7599 ext 84611
</pre>
  </font></span></div>

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