<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<br>
<div class="moz-cite-prefix">On 05/04/14 02:31, Matthew Knepley
wrote:<br>
</div>
<blockquote
cite="mid:CAMYG4G=k3QN1DA6yHSk99FxakrmrWjO653YnJt9e16G5KgUoEQ@mail.gmail.com"
type="cite">
<meta http-equiv="Content-Type" content="text/html;
charset=ISO-8859-1">
<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 class="moz-txt-link-freetext" href="http://http://www.mcs.anl.gov/petsc/documentation/faq.html">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>
<br>
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>
<br>
- 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<br>
<br>
<br>
<br>
<pre class="moz-signature" cols="72">--
Dr Adrian Croucher
Department of Engineering Science
University of Auckland
New Zealand
tel 64-9-373-7599 ext 84611
</pre>
</body>
</html>