[petsc-users] Accessing Global Vector Data with

[Daniel Mckinnell via petsc-users]

Hi,Attached is a way I came up with to access data in a global vector, is this 
the best way to do this or are there other ways? It would seem intuitive to use 
the global PetscSection and VecGetValuesSection but this doesn't seem to work 
on global vectors.
Instead I have used VecGetValues and VecSetValues, however I also have a 
problem with these when extracting more than one value, I initialise the output 
of VecGetValues as PetscScalar *values; and then call VecGetValues(Other 
stuff... , values). This seems to work some times and not others and I can't 
find any rhyme or reason to it?

Finally I was wondering if there is a good reference code base on Github 
including DMPlex that would be helpful in viewing applications of the DMPlex 
functions?Thanks,Daniel Mckinnell

#include <petscdmplex.h>
#include <petscdmlabel.h>
#include <petscsf.h>
#include <iostream>

PetscErrorCode dm_set_cell_values(DM *dm, Vec *v){
    Vec u;
    PetscSection    globalsection;
    PetscInt        cStart, cEnd;

    DMGetGlobalVector(*dm, &u);
    DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
    DMGetGlobalSection(*dm,&globalsection);
    PetscSectionView(globalsection, PETSC_VIEWER_STDOUT_WORLD);
    for (int i = cStart; i < cEnd; i++){
        PetscInt     csize, coff;
        PetscSectionGetDof(globalsection, i, &csize);
        PetscSectionGetOffset(globalsection,i,&coff);
        PetscInt extractvalues[2] = {coff, coff+1};
        if (coff>=0){
          for (int j=0;j<csize;j++){
            PetscScalar value;
            PetscInt extract=coff+j, num=1;
            VecGetValues(u, num, &extract, &value);
            if (j==1){
              value += i+1;
            }
            VecSetValues(u, num, &extract, &value, INSERT_VALUES);
          }
        }
    }
    *v=u;
}

int main(int argc, char **argv)
{
  DM dm, dmDist = NULL, dmtest;
  PetscSF sf, pointsf;
  Vec v;
  PetscSection section, globalsection;
  PetscInt dim = 2, numFields, numBC, i, cStart, cEnd;
  PetscScalar *u_value;
  PetscInt numComp[3];
  PetscInt numDof[12];
  PetscInt bcField[1];
  IS bcPointIS[1];
  PetscBool interpolate = PETSC_TRUE;
  PetscScalar *u_array;

  try
  {
    PetscInitialize(&argc, &argv, NULL, NULL);
  }
  catch (const std::exception &e)
  {
    return (-1);
  };
  PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL);
  DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, NULL, NULL, NULL, NULL, interpolate, &dm);
  DMPlexDistribute(dm, 0, &sf, &dmDist);
  if (dmDist)
  {
    DMDestroy(&dm);
    dm = dmDist;
  }
  numFields = 3;
  numComp[0] = 1;
  numComp[1] = dim;
  numComp[2] = dim - 1;
  for (i = 0; i < numFields * (dim + 1); ++i)
    numDof[i] = 0;
  numDof[0 * (dim + 1) + 0] = 1;
  numDof[1 * (dim + 1) + dim] = dim;
  numDof[2 * (dim + 1) + dim - 1] = dim - 1;
  numBC = 1;
  bcField[0] = 0;

  DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]);
  DMSetNumFields(dm, numFields);
  DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);
  ISDestroy(&bcPointIS[0]);
  /* Name the Field variables */
  PetscSectionSetFieldName(section, 0, "u");
  PetscSectionSetFieldName(section, 1, "v");
  PetscSectionSetFieldName(section, 2, "w");
  DMSetSection(dm, section);

  DMGetPointSF(dm, &pointsf);
  PetscSectionCreateGlobalSection(section, pointsf, PETSC_TRUE, PETSC_FALSE, &globalsection);
  DMSetGlobalSection(dm, globalsection);
  dm_set_cell_values(&dm,&v);
  PetscSectionDestroy(&globalsection);
  PetscSectionDestroy(&section);
  DMDestroy(&dm);
  PetscFinalize();
  return (0);
}