[petsc-users] Question about DMDAGetElements

[Jed Brown]

Sajid Ali via petsc-users <petsc-users at mcs.anl.gov> writes:

> Hi,
>
> I'm trying to understand this example
> <https://github.com/haplav/petsc-tut-it4i-2018/blob/master/exercises/ex6.c>
> from a tutorial. DM is used as an alternative to setting up the matrix and
> vector separately (as was done in the previous example
> <https://github.com/haplav/petsc-tut-it4i-2018/blob/master/exercises/ex5.c>).
>
>
> Since DMDACreate1d was used to create the mesh, can someone explain this
> logic more clearly :
>
>   /*
>     Gets an array containing the indices (in local coordinates)
>     of all the local elements.
>     nel	- number of local elements
>     nen	- number of one element's nodes
>     e	- the local indices of the elements' vertices
>   */
>   ierr = DMDAGetElements(da, &nel, &nen, &e);CHKERRQ(ierr);
>
>   /*
>      Assemble matrix and RHS
>   */
>   value[0] = 1.0; value[1] = -1.0; value[2] = -1.0; value[3] = 1.0;
>   bvalue[0] = 1.0; bvalue[1] = 1.0;
>   for (i=0; i<nel; i++) {
>     ierr = MatSetValuesLocal(A, 2, e+nen*i, 2, e+nen*i, value,
> ADD_VALUES);CHKERRQ(ierr);
>     //TODO use VecSetValuesLocal() to set the values
>   }
>   ierr = DMDARestoreElements(da, &nel, &nen, &e);CHKERRQ(ierr);
>
>
> The dm was created using
>
>   ierr = MPI_Allreduce(&n, &N, 1, MPI_INT, MPI_SUM,
> PETSC_COMM_WORLD);CHKERRQ(ierr);
>   ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&da);CHKERRQ(ierr);

So each process owns 3 vertices.  The first process owns two elements,
then each subsequent process gets three elements.  Only the first
nel*nen entries in that array are valid.  nen=2 in all cases because
these are 1D linear elements.

Use DMPlex if you want more general finite element mesh handling.

>
> Using n = 3 and I get the following results when I print the values of
> nel,nen and e using this statement
> PetscPrintf(PETSC_COMM_SELF,"nel : %d,nen : %d, e elements :
> %d,%d,%d,%d,%d,%d,%d,%d ,rank
> %d\n",nel,nen,e[0],e[1],e[2],e[3],e[4],e[5],e[6],e[7],rank);
>
> The output comes out as  :
>
> [sajid at xrm exercises]$ mpirun -np 4 ./ex6
> nel : 2,nen : 2, e elements : 0,1,1,2,1936617321,0,3,0 ,rank 0
> nel : 3,nen : 2, e elements : 0,1,1,2,2,3,4,0 ,rank 1
> nel : 3,nen : 2, e elements : 0,1,1,2,2,3,4,0 ,rank 2
> nel : 3,nen : 2, e elements : 0,1,1,2,2,3,3,0 ,rank 3
>
> Does this mean that I'm supposed to read the first 4 elements for rank 0 as
> row0, col0, row1, col1?
>
> And how does differ If I'm creating a vector using DMCreate1d ?
>
>
> Thank You,
> Sajid Ali
> Applied Physics
> Northwestern University