[petsc-users] Altering Vertex Coordinates of an Unstructured Mesh

[Daniel Mckinnell via petsc-users]

Hello,I am a student currently doing a summer project so am fairly new to 
Petsc. I am trying to alter vertex coordinates of an unstructured mesh to look 
at mesh movement in DMPlex
after they have been created by DMPlexCreateBoxMesh. I have been altering the 
first example on the 
DMPlex documentation here: 
https://www.mcs.anl.gov/petsc/petsc-current/src/dm/impls/plex/examples/tutorials/ex1.c.html
 
I have written something that vaguely does what I want in parallel, but it has 
a problem as the output from DMPlexGetPointLocal 
seems to be very inconsistent. Below is my code, I was wondering if anyone had 
a better way of doing this? I believeI have a bug but I am not sure where to 
start. Any help and advice/comments very welcome.
Thanks,
Daniel
#include <petscdmplex.h>
#include <Eigen/Core>
#include <iostream>

using namespace Eigen;

int move_vertex(DM dm, int global_index)
{
    PetscInt    n, extract, st, end, dim = 2;
    IS          global_num;
    Vec         u;
    const PetscInt *nindices;
    PetscScalar    u_array;

    DMGetCoordinatesLocal(dm, &u);
    DMPlexGetVertexNumbering(dm, &global_num);

    bool local_point = false;
    ISGetLocalSize(global_num,&n);
    ISGetIndices(global_num,&nindices);
    for (int j=0; j<n; j++){
        if (nindices[j]==global_index){
            local_point = true;
            extract = j;
            break;
        } // Finds the processor that the vertex is assigned to and gets the local index of this vertex 
    }

    if (local_point){
        DMPlexGetPointLocal(dm, extract, &st, &end);
        PetscInt ix[2] = {st,st+1};
        VecGetValues(u,dim,ix,&u_array);
        Map<VectorXd> test(&u_array,2);
        test(0)+=-0.2; // Subtract 0.2 from x coord as a sample test
        VecSetValues(u,dim,ix,&u_array,INSERT_VALUES);
    } // Extracts the coordinates of the vertex, alters it and then writes it back to the Vec

    DMSetCoordinatesLocal(dm, u);
    return(0);
}

int main(int argc, char **argv)
{
DM             dm, dmDist = NULL;
Vec            c, x;
PetscSection   section, sect1;
PetscInt       dim = 2, numFields, numBC, i, ls;
PetscInt       numComp[3];
PetscInt       numDof[12];
PetscInt       bcField[1];
IS             bcPointIS[1];
PetscBool      interpolate = PETSC_TRUE, flag0, flag1, flag2;
PetscScalar    *c_array, v;
PetscScalar check1[18] = {-0.2,0.0, 0.5,0.0, 1.0,0.0, 0.0,0.5, 0.5,0.5, 1.0,0.5, 0.0,1.0, 0.5,1.0, 1.0,1.0};

try {
     PetscInitialize(&argc, &argv, NULL, NULL);
}
catch (const std::exception& e) {return(0);};
PetscOptionsGetInt(NULL,NULL, "-dim", &dim, NULL);
DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_FALSE, NULL, NULL, NULL, NULL, interpolate, &dm);
DMPlexDistribute(dm, 0, NULL, &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;
/* Let u be defined on vertices */
numDof[0*(dim+1)+0]     = 1;
/* Let v be defined on cells */
numDof[1*(dim+1)+dim]   = dim;
/* Let v be defined on faces */
numDof[2*(dim+1)+dim-1] = dim-1;
/* Setup boundary conditions */
numBC = 1;
/* Prescribe a Dirichlet condition on u on the boundary
     Label "marker" is made by the mesh creation routine */
bcField[0] = 0;
DMGetStratumIS(dm, "marker", 1, &bcPointIS[0]);
/* Create a PetscSection with this data layout */
DMSetNumFields(dm, numFields);
DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);

ISDestroy(&bcPointIS[0]);
DMSetSection(dm, section);

move_vertex(dm, 0);
DMGetCoordinates(dm, &c);

VecCreate(PETSC_COMM_WORLD,&x);
VecGetLocalSize(c, &ls);
VecSetSizes(x,ls,PETSC_DECIDE);
VecSetFromOptions(x);
VecSet(x,v);
PetscInt j=0;
for (j=0; j<18; j++) {
v = check1[j];
VecSetValues(x,1,&j,&v,INSERT_VALUES);
}

// Tests I would like it to pass
VecAssemblyBegin(x);
VecAssemblyEnd(x);
VecEqual(x,c, &flag0); // Moving a vertex on the mesh boundary
assert(flag0 == PETSC_TRUE);
DMSetCoordinates(dm, c);

move_vertex(dm, 4);
DMGetCoordinates(dm, &c);
j=8;
v=check1[j]-0.2;
VecSetValues(x,1,&j,&v,INSERT_VALUES);
VecAssemblyBegin(x);
VecAssemblyEnd(x);

VecEqual(x,c, &flag1); // Moving a vertex on the processor boundary
assert(flag1 == PETSC_TRUE);
DMSetCoordinates(dm, c);

move_vertex(dm, 5);
DMGetCoordinates(dm, &c);
j=10;
v=check1[j]-0.2;
VecSetValues(x,1,&j,&v,INSERT_VALUES);
VecAssemblyBegin(x);
VecAssemblyEnd(x);

VecEqual(x,c, &flag2); //Moving a vertex on both the processor boundary and the mesh boundary
assert(flag2 == PETSC_TRUE);
DMSetCoordinates(dm, c);

DMDestroy(&dm);
PetscFinalize();
return(0);
}