Dear Matt, Dear Petsc

After my recent pull of petsc master I encounter test fails that use
dmplex meshes.

I attached a simple script to show an issue where DMPlexRestoreCone
mangles the cone relationships in the DM.

The behaviour happens with gcc 7.4.0 aswell as icc 19.0.3.199
@ current petsc/master d3dc048301fdb3e4c282b528e7b28910551cfe13

We recently had a system update and I am not 100% sure if this is some
weird stuff going on at my end or if this is indeed Plex related or if I
should not call Restore in the inner loop.

Anyway, I would be most grateful if you could try to reproduce the error
and hint me towards a solution.

As always, many many thanks!

Fabian

P.S. here the output of
make run

and note the cones, the second time they are outputted...

./plex_prism -show_plex ::ascii_info_detail
DM Object: testplex 1 MPI processes
  type: plex
testplex in 3 dimensions:
Supports:
[0] Max support size: 3
[0]: 1 ----> 0
[0]: 2 ----> 0
[0]: 3 ----> 0
[0]: 4 ----> 0
[0]: 5 ----> 0
[0]: 6 ----> 1
[0]: 6 ----> 4
[0]: 7 ----> 1
[0]: 7 ----> 3
[0]: 8 ----> 1
[0]: 8 ----> 5
[0]: 9 ----> 3
[0]: 9 ----> 4
[0]: 10 ----> 4
[0]: 10 ----> 5
[0]: 11 ----> 3
[0]: 11 ----> 5
[0]: 12 ----> 2
[0]: 12 ----> 4
[0]: 13 ----> 2
[0]: 13 ----> 3
[0]: 14 ----> 2
[0]: 14 ----> 5
[0]: 15 ----> 6
[0]: 15 ----> 7
[0]: 15 ----> 9
[0]: 16 ----> 6
[0]: 16 ----> 8
[0]: 16 ----> 10
[0]: 17 ----> 7
[0]: 17 ----> 8
[0]: 17 ----> 11
[0]: 18 ----> 9
[0]: 18 ----> 12
[0]: 18 ----> 13
[0]: 19 ----> 10
[0]: 19 ----> 12
[0]: 19 ----> 14
[0]: 20 ----> 11
[0]: 20 ----> 13
[0]: 20 ----> 14
Cones:
[0] Max cone size: 5
[0]: 0 <---- 1 (0)
[0]: 0 <---- 2 (0)
[0]: 0 <---- 3 (0)
[0]: 0 <---- 4 (0)
[0]: 0 <---- 5 (0)
[0]: 1 <---- 6 (0)
[0]: 1 <---- 7 (0)
[0]: 1 <---- 8 (0)
[0]: 2 <---- 12 (0)
[0]: 2 <---- 13 (0)
[0]: 2 <---- 14 (0)
[0]: 3 <---- 7 (0)
[0]: 3 <---- 9 (0)
[0]: 3 <---- 11 (0)
[0]: 3 <---- 13 (0)
[0]: 4 <---- 6 (0)
[0]: 4 <---- 9 (0)
[0]: 4 <---- 10 (0)
[0]: 4 <---- 12 (0)
[0]: 5 <---- 8 (0)
[0]: 5 <---- 10 (0)
[0]: 5 <---- 11 (0)
[0]: 5 <---- 14 (0)
[0]: 6 <---- 15 (0)
[0]: 6 <---- 16 (0)
[0]: 7 <---- 15 (0)
[0]: 7 <---- 17 (0)
[0]: 8 <---- 16 (0)
[0]: 8 <---- 17 (0)
[0]: 9 <---- 15 (0)
[0]: 9 <---- 18 (0)
[0]: 10 <---- 16 (0)
[0]: 10 <---- 19 (0)
[0]: 11 <---- 17 (0)
[0]: 11 <---- 20 (0)
[0]: 12 <---- 18 (0)
[0]: 12 <---- 19 (0)
[0]: 13 <---- 18 (0)
[0]: 13 <---- 20 (0)
[0]: 14 <---- 19 (0)
[0]: 14 <---- 20 (0)
coordinates with 1 fields
  field 0 with 3 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  0 offset   0
  (   9) dim  0 offset   0
  (  10) dim  0 offset   0
  (  11) dim  0 offset   0
  (  12) dim  0 offset   0
  (  13) dim  0 offset   0
  (  14) dim  0 offset   0
  (  15) dim  3 offset   0 0. 0. 0.
  (  16) dim  3 offset   3 0. 2. 0.
  (  17) dim  3 offset   6 2. 1. 0.
  (  18) dim  3 offset   9 0. 0. 1.
  (  19) dim  3 offset  12 0. 2. 1.
  (  20) dim  3 offset  15 2. 1. 1.









DM Object: testplex 1 MPI processes
  type: plex
testplex in 3 dimensions:
Supports:
[0] Max support size: 3
[0]: 1 ----> 0
[0]: 2 ----> 0
[0]: 3 ----> 0
[0]: 4 ----> 0
[0]: 5 ----> 0
[0]: 6 ----> 1
[0]: 6 ----> 4
[0]: 7 ----> 1
[0]: 7 ----> 3
[0]: 8 ----> 1
[0]: 8 ----> 5
[0]: 9 ----> 3
[0]: 9 ----> 4
[0]: 10 ----> 4
[0]: 10 ----> 5
[0]: 11 ----> 3
[0]: 11 ----> 5
[0]: 12 ----> 2
[0]: 12 ----> 4
[0]: 13 ----> 2
[0]: 13 ----> 3
[0]: 14 ----> 2
[0]: 14 ----> 5
[0]: 15 ----> 6
[0]: 15 ----> 7
[0]: 15 ----> 9
[0]: 16 ----> 6
[0]: 16 ----> 8
[0]: 16 ----> 10
[0]: 17 ----> 7
[0]: 17 ----> 8
[0]: 17 ----> 11
[0]: 18 ----> 9
[0]: 18 ----> 12
[0]: 18 ----> 13
[0]: 19 ----> 10
[0]: 19 ----> 12
[0]: 19 ----> 14
[0]: 20 ----> 11
[0]: 20 ----> 13
[0]: 20 ----> 14
Cones:
[0] Max cone size: 5
[0]: 0 <---- 0 (0)
[0]: 0 <---- 0 (0)
[0]: 0 <---- 3 (0)
[0]: 0 <---- 4 (0)
[0]: 0 <---- 5 (0)
[0]: 1 <---- 0 (0)
[0]: 1 <---- 0 (0)
[0]: 1 <---- 8 (0)
[0]: 2 <---- 0 (0)
[0]: 2 <---- 0 (0)
[0]: 2 <---- 14 (0)
[0]: 3 <---- 0 (0)
[0]: 3 <---- 0 (0)
[0]: 3 <---- 11 (0)
[0]: 3 <---- 13 (0)
[0]: 4 <---- 0 (0)
[0]: 4 <---- 0 (0)
[0]: 4 <---- 10 (0)
[0]: 4 <---- 12 (0)
[0]: 5 <---- 0 (0)
[0]: 5 <---- 0 (0)
[0]: 5 <---- 11 (0)
[0]: 5 <---- 14 (0)
[0]: 6 <---- 0 (0)
[0]: 6 <---- 0 (0)
[0]: 7 <---- 0 (0)
[0]: 7 <---- 0 (0)
[0]: 8 <---- 0 (0)
[0]: 8 <---- 0 (0)
[0]: 9 <---- 0 (0)
[0]: 9 <---- 0 (0)
[0]: 10 <---- 0 (0)
[0]: 10 <---- 0 (0)
[0]: 11 <---- 0 (0)
[0]: 11 <---- 0 (0)
[0]: 12 <---- 0 (0)
[0]: 12 <---- 0 (0)
[0]: 13 <---- 0 (0)
[0]: 13 <---- 0 (0)
[0]: 14 <---- 0 (0)
[0]: 14 <---- 0 (0)
coordinates with 1 fields
  field 0 with 3 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  0 offset   0
  (   9) dim  0 offset   0
  (  10) dim  0 offset   0
  (  11) dim  0 offset   0
  (  12) dim  0 offset   0
  (  13) dim  0 offset   0
  (  14) dim  0 offset   0
  (  15) dim  3 offset   0 0. 0. 0.
  (  16) dim  3 offset   3 0. 2. 0.
  (  17) dim  3 offset   6 2. 1. 0.
  (  18) dim  3 offset   9 0. 0. 1.
  (  19) dim  3 offset  12 0. 2. 1.
  (  20) dim  3 offset  15 2. 1. 1.

./plex_prism -show_plex ::ascii_info_detail
DM Object: testplex 1 MPI processes
  type: plex
testplex in 3 dimensions:
Supports:
[0] Max support size: 3
[0]: 1 ----> 0
[0]: 2 ----> 0
[0]: 3 ----> 0
[0]: 4 ----> 0
[0]: 5 ----> 0
[0]: 6 ----> 1
[0]: 6 ----> 4
[0]: 7 ----> 1
[0]: 7 ----> 3
[0]: 8 ----> 1
[0]: 8 ----> 5
[0]: 9 ----> 3
[0]: 9 ----> 4
[0]: 10 ----> 4
[0]: 10 ----> 5
[0]: 11 ----> 3
[0]: 11 ----> 5
[0]: 12 ----> 2
[0]: 12 ----> 4
[0]: 13 ----> 2
[0]: 13 ----> 3
[0]: 14 ----> 2
[0]: 14 ----> 5
[0]: 15 ----> 6
[0]: 15 ----> 7
[0]: 15 ----> 9
[0]: 16 ----> 6
[0]: 16 ----> 8
[0]: 16 ----> 10
[0]: 17 ----> 7
[0]: 17 ----> 8
[0]: 17 ----> 11
[0]: 18 ----> 9
[0]: 18 ----> 12
[0]: 18 ----> 13
[0]: 19 ----> 10
[0]: 19 ----> 12
[0]: 19 ----> 14
[0]: 20 ----> 11
[0]: 20 ----> 13
[0]: 20 ----> 14
Cones:
[0] Max cone size: 5
[0]: 0 <---- 1 (0)
[0]: 0 <---- 2 (0)
[0]: 0 <---- 3 (0)
[0]: 0 <---- 4 (0)
[0]: 0 <---- 5 (0)
[0]: 1 <---- 6 (0)
[0]: 1 <---- 7 (0)
[0]: 1 <---- 8 (0)
[0]: 2 <---- 12 (0)
[0]: 2 <---- 13 (0)
[0]: 2 <---- 14 (0)
[0]: 3 <---- 7 (0)
[0]: 3 <---- 9 (0)
[0]: 3 <---- 11 (0)
[0]: 3 <---- 13 (0)
[0]: 4 <---- 6 (0)
[0]: 4 <---- 9 (0)
[0]: 4 <---- 10 (0)
[0]: 4 <---- 12 (0)
[0]: 5 <---- 8 (0)
[0]: 5 <---- 10 (0)
[0]: 5 <---- 11 (0)
[0]: 5 <---- 14 (0)
[0]: 6 <---- 15 (0)
[0]: 6 <---- 16 (0)
[0]: 7 <---- 15 (0)
[0]: 7 <---- 17 (0)
[0]: 8 <---- 16 (0)
[0]: 8 <---- 17 (0)
[0]: 9 <---- 15 (0)
[0]: 9 <---- 18 (0)
[0]: 10 <---- 16 (0)
[0]: 10 <---- 19 (0)
[0]: 11 <---- 17 (0)
[0]: 11 <---- 20 (0)
[0]: 12 <---- 18 (0)
[0]: 12 <---- 19 (0)
[0]: 13 <---- 18 (0)
[0]: 13 <---- 20 (0)
[0]: 14 <---- 19 (0)
[0]: 14 <---- 20 (0)
coordinates with 1 fields
  field 0 with 3 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  0 offset   0
  (   9) dim  0 offset   0
  (  10) dim  0 offset   0
  (  11) dim  0 offset   0
  (  12) dim  0 offset   0
  (  13) dim  0 offset   0
  (  14) dim  0 offset   0
  (  15) dim  3 offset   0 0. 0. 0.
  (  16) dim  3 offset   3 0. 2. 0.
  (  17) dim  3 offset   6 2. 1. 0.
  (  18) dim  3 offset   9 0. 0. 1.
  (  19) dim  3 offset  12 0. 2. 1.
  (  20) dim  3 offset  15 2. 1. 1.
DM Object: testplex 1 MPI processes
  type: plex
testplex in 3 dimensions:
Supports:
[0] Max support size: 3
[0]: 1 ----> 0
[0]: 2 ----> 0
[0]: 3 ----> 0
[0]: 4 ----> 0
[0]: 5 ----> 0
[0]: 6 ----> 1
[0]: 6 ----> 4
[0]: 7 ----> 1
[0]: 7 ----> 3
[0]: 8 ----> 1
[0]: 8 ----> 5
[0]: 9 ----> 3
[0]: 9 ----> 4
[0]: 10 ----> 4
[0]: 10 ----> 5
[0]: 11 ----> 3
[0]: 11 ----> 5
[0]: 12 ----> 2
[0]: 12 ----> 4
[0]: 13 ----> 2
[0]: 13 ----> 3
[0]: 14 ----> 2
[0]: 14 ----> 5
[0]: 15 ----> 6
[0]: 15 ----> 7
[0]: 15 ----> 9
[0]: 16 ----> 6
[0]: 16 ----> 8
[0]: 16 ----> 10
[0]: 17 ----> 7
[0]: 17 ----> 8
[0]: 17 ----> 11
[0]: 18 ----> 9
[0]: 18 ----> 12
[0]: 18 ----> 13
[0]: 19 ----> 10
[0]: 19 ----> 12
[0]: 19 ----> 14
[0]: 20 ----> 11
[0]: 20 ----> 13
[0]: 20 ----> 14
Cones:
[0] Max cone size: 5
[0]: 0 <---- 0 (0)
[0]: 0 <---- 0 (0)
[0]: 0 <---- 3 (0)
[0]: 0 <---- 4 (0)
[0]: 0 <---- 5 (0)
[0]: 1 <---- 0 (0)
[0]: 1 <---- 0 (0)
[0]: 1 <---- 8 (0)
[0]: 2 <---- 0 (0)
[0]: 2 <---- 0 (0)
[0]: 2 <---- 14 (0)
[0]: 3 <---- 0 (0)
[0]: 3 <---- 0 (0)
[0]: 3 <---- 11 (0)
[0]: 3 <---- 13 (0)
[0]: 4 <---- 0 (0)
[0]: 4 <---- 0 (0)
[0]: 4 <---- 10 (0)
[0]: 4 <---- 12 (0)
[0]: 5 <---- 0 (0)
[0]: 5 <---- 0 (0)
[0]: 5 <---- 11 (0)
[0]: 5 <---- 14 (0)
[0]: 6 <---- 0 (0)
[0]: 6 <---- 0 (0)
[0]: 7 <---- 0 (0)
[0]: 7 <---- 0 (0)
[0]: 8 <---- 0 (0)
[0]: 8 <---- 0 (0)
[0]: 9 <---- 0 (0)
[0]: 9 <---- 0 (0)
[0]: 10 <---- 0 (0)
[0]: 10 <---- 0 (0)
[0]: 11 <---- 0 (0)
[0]: 11 <---- 0 (0)
[0]: 12 <---- 0 (0)
[0]: 12 <---- 0 (0)
[0]: 13 <---- 0 (0)
[0]: 13 <---- 0 (0)
[0]: 14 <---- 0 (0)
[0]: 14 <---- 0 (0)
coordinates with 1 fields
  field 0 with 3 components
Process 0:
  (   0) dim  0 offset   0
  (   1) dim  0 offset   0
  (   2) dim  0 offset   0
  (   3) dim  0 offset   0
  (   4) dim  0 offset   0
  (   5) dim  0 offset   0
  (   6) dim  0 offset   0
  (   7) dim  0 offset   0
  (   8) dim  0 offset   0
  (   9) dim  0 offset   0
  (  10) dim  0 offset   0
  (  11) dim  0 offset   0
  (  12) dim  0 offset   0
  (  13) dim  0 offset   0
  (  14) dim  0 offset   0
  (  15) dim  3 offset   0 0. 0. 0.
  (  16) dim  3 offset   3 0. 2. 0.
  (  17) dim  3 offset   6 2. 1. 0.
  (  18) dim  3 offset   9 0. 0. 1.
  (  19) dim  3 offset  12 0. 2. 1.
  (  20) dim  3 offset  15 2. 1. 1.


 Cone           0 :           1           2           3           4           5
 Cone           1 :           6           7           8
 Cone           2 :          12          13          14
 Cone           3 :           7           9          11          13
 Cone           4 :           6           9          10          12
 Cone           5 :           8          10          11          14
 Cone           6 :          15          16
 Cone           7 :          15          17
 Cone           8 :          16          17
 Cone           9 :          15          18
 Cone          10 :          16          19
 Cone          11 :          17          20
 Cone          12 :          18          19
 Cone          13 :          18          20
 Cone          14 :          19          20
 Cone          15 :
 Cone          16 :
 Cone          17 :
 Cone          18 :
 Cone          19 :
 Cone          20 :
include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules

run:: plex_prism
                ./plex_prism -show_plex ::ascii_info_detail

plex_prism:: plex_prism.F90
                ${PETSC_FCOMPILE} -c plex_prism.F90
                ${FLINKER} plex_prism.o -o plex_prism ${PETSC_LIB}

clean::
                rm -rf *.o prism.h5 prism.xmf plex_prism
program main
#include "petsc/finclude/petsc.h"

      use petsc
      implicit none

      PetscErrorCode :: ierr
      PetscInt :: i

      DM :: dm

      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

      call create_plex()
      call PetscObjectViewFromOptions(dm, PETSC_NULL_VEC, "-show_plex", ierr)
      call print_cones()
      call PetscObjectViewFromOptions(dm, PETSC_NULL_VEC, "-show_plex", ierr)

      call DMDestroy(dm, ierr)
      call PetscFinalize(ierr)

      contains
        subroutine print_cones()
          PetscInt :: pstart, pEnd
          PetscInt, pointer :: cone(:)
          call DMPlexGetChart(dm, pStart, pEnd, ierr); CHKERRQ(ierr)
          do i = pStart, pEnd-1
            call DMPlexGetCone(dm, i, cone, ierr); CHKERRQ(ierr)
            print *,'Cone',i,':',cone
            ! Next line is troublesome because it mangles the DM
            call DMPlexRestoreCone(dm, i, cone, ierr); CHKERRQ(ierr)
          enddo
        end subroutine

        subroutine create_plex()
          call DMPlexCreate(PETSC_COMM_WORLD, dm, ierr);CHKERRQ(ierr)
          call PetscObjectSetName(dm, 'testplex', ierr);CHKERRQ(ierr)
          call DMSetDimension(dm, 3, ierr);CHKERRQ(ierr)

          call DMPlexSetChart(dm, 0, 21, ierr); CHKERRQ(ierr)
          ! Preallocation
          ! cell has 5 faces
          call DMPlexSetConeSize(dm, 0, 5, ierr); CHKERRQ(ierr)

          ! Faces have 3 or 4 edges
          do i=1,2
            call DMPlexSetConeSize(dm, i, 3, ierr); CHKERRQ(ierr)
          enddo
          do i=3,5
            call DMPlexSetConeSize(dm, i, 4, ierr); CHKERRQ(ierr)
          enddo

          ! Edges have 2 vertices
          do i=6,14
            call DMPlexSetConeSize(dm, i, 2, ierr); CHKERRQ(ierr)
          enddo

          call DMSetUp(dm, ierr); CHKERRQ(ierr) ! Allocate space for cones

          ! Setup Connections
          ! 0 cell has faces 1(bottom) and 2(top) as well as three faces on sides
          call DMPlexSetCone(dm,  0, [1,2,3,4,5], ierr); CHKERRQ(ierr)

          ! bottom face has 3 edges
          call DMPlexSetCone(dm,  1, [ 6, 7, 8], ierr); CHKERRQ(ierr)
          ! top face has 3 edges
          call DMPlexSetCone(dm,  2, [12,13,14], ierr); CHKERRQ(ierr)

          ! side faces have 4 edges
          call DMPlexSetCone(dm,  3, [ 7, 9,11,13], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm,  4, [ 6, 9,10,12], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm,  5, [ 8,10,11,14], ierr); CHKERRQ(ierr)

          ! and the corresponding vertices
          call DMPlexSetCone(dm,  6, [15,16], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm,  7, [15,17], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm,  8, [16,17], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm,  9, [15,18], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm, 10, [16,19], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm, 11, [17,20], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm, 12, [18,19], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm, 13, [18,20], ierr); CHKERRQ(ierr)
          call DMPlexSetCone(dm, 14, [19,20], ierr); CHKERRQ(ierr)

          call DMPlexSymmetrize(dm, ierr); CHKERRQ(ierr)
          call DMPlexStratify(dm, ierr); CHKERRQ(ierr)

          call set_coords()
        end subroutine

        subroutine set_coords()
          PetscScalar, pointer :: coords(:)
          Vec                  :: coordinates
          PetscInt             :: dimEmbed, coordSize, vStart, vEnd, pStart, pEnd, voff
          PetscSection         :: coordSection

          call DMGetCoordinateDim(dm, dimEmbed, ierr); CHKERRQ(ierr)

          call DMGetCoordinateSection(dm, coordSection, ierr); CHKERRQ(ierr)

          call PetscSectionSetNumFields(coordSection, 1, ierr); CHKERRQ(ierr)
          call PetscSectionSetUp(coordSection, ierr); CHKERRQ(ierr)
          call PetscSectionSetFieldComponents(coordSection, 0, dimEmbed, ierr); CHKERRQ(ierr)

          call DMPlexGetChart(dm, pStart, pEnd, ierr); CHKERRQ(ierr)
          call DMPlexGetDepthStratum (dm, 0, vStart, vEnd, ierr); CHKERRQ(ierr) ! vertices

          call PetscSectionSetChart(coordSection, pStart, pEnd, ierr);CHKERRQ(ierr)

          do i = vStart, vEnd-1
            call PetscSectionSetDof(coordSection, i, dimEmbed, ierr); CHKERRQ(ierr)
            call PetscSectionSetFieldDof(coordSection, i, 0, dimEmbed, ierr); CHKERRQ(ierr)
          enddo

          call PetscSectionSetUp(coordSection, ierr); CHKERRQ(ierr)
          call PetscSectionGetStorageSize(coordSection, coordSize, ierr); CHKERRQ(ierr)

          call VecCreate(PETSC_COMM_SELF, coordinates, ierr); CHKERRQ(ierr)
          call VecSetSizes(coordinates, coordSize, PETSC_DETERMINE, ierr);CHKERRQ(ierr)
          call VecSetBlockSize(coordinates, dimEmbed, ierr);CHKERRQ(ierr)
          call VecSetType(coordinates, VECSTANDARD, ierr);CHKERRQ(ierr)

          call PetscObjectSetName(coordinates, "coordinates", ierr); CHKERRQ(ierr)

          call VecGetArrayF90(coordinates, coords, ierr); CHKERRQ(ierr)

          i = 15; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,0,0]
          i = 16; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,2,0]
          i = 17; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [2,1,0]
          i = 18; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,0,1]
          i = 19; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [0,2,1]
          i = 20; call PetscSectionGetOffset(coordSection, i, voff, ierr); coords(voff+1:voff+3) = [2,1,1]

          call VecRestoreArrayF90(coordinates, coords, ierr); CHKERRQ(ierr)
          call DMSetCoordinatesLocal(dm, coordinates, ierr);CHKERRQ(ierr)
          call PetscObjectViewFromOptions(coordinates, PETSC_NULL_VEC, "-show_plex_coordinates", ierr); CHKERRQ(ierr)
          call VecDestroy(coordinates, ierr);CHKERRQ(ierr)
        end subroutine

      end program

Reply via email to