[petsc-users] Interpolation using DMInterpolationEvaluate

[Swarnava Ghosh]

Hi PETSc developers and users,

I am trying to setup a test case of DMInterpolationEvaluate. I am referring
to the example provided in (
https://bitbucket.org/petsc/petsc/src/master/src/snes/examples/tests/ex2.c)
to code it up.

In my case. I have read a mesh from a file created using Gmesh, and I am
trying to interpolate a constant field of ones using
DMInterpolationEvaluate. Since the field is constant, then the interpolated
values should also be the same. However, I am getting other values:

Attached is the code which I wrote to test it:

static char help[] = "Test interpolation\n\
 using a parallel unstructured mesh (DMPLEX) to discretize it.\n";

#include <petsctime.h>
#include <mpi.h>
#include <petsc.h>
#include <petscdmda.h>
#include <petscdmplex.h>
#include <petscds.h>

#undef __FUNCT__
#define __FUNCT__ "main"

PetscErrorCode SolveProblem(int i)
{

  PetscBool interpolate=0;
  char filename[]="mesh2.msh";
  //char filename[]="mesh1";
  PetscErrorCode ierr;
  DM dm;
  DM               distributedMesh = NULL;
  PetscViewer    viewer;
  DMLabel   label;
  PetscBool      simplex=PETSC_TRUE;
  PetscFE         fe;
  PetscQuadrature q;
  PetscDS        prob=NULL;
  PetscInt id = 1;
  PetscInt order;
  PetscPartitioner part;
  SNES           snes;        /* nonlinear solver */
  int rank;
  Vec       fieldvec;


  MPI_Comm_rank(MPI_COMM_WORLD,&rank);

  // read mesh from file

ierr=DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,filename,interpolate,&dm);CHKERRQ(ierr);

  //    Distribute mesh over processes
  ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr);
  ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
  ierr = DMPlexDistribute(dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);

  if (distributedMesh) {
   ierr = DMDestroy(&dm);CHKERRQ(ierr);
    dm  = distributedMesh;
  }



  Vec vcoord;
  DMGetCoordinatesLocal(dm,&vcoord);
  VecView(vcoord,PETSC_VIEWER_STDOUT_SELF);


  // interpolation
  PetscInt            spaceDim, c, Np=8, p;
  DMInterpolationInfo interpolator;
  PetscReal          pcoords[16];
  PetscBool           pointsAllProcs=PETSC_TRUE;
  Vec                 lu, fieldVals;
  pcoords[0]=0.0; pcoords[1]=0;
  pcoords[2]=1; pcoords[3]=0;
  pcoords[4]=1; pcoords[5]=1;
  pcoords[6]=0; pcoords[7]=1;
  pcoords[8]=0.5; pcoords[9]=0.5;
  pcoords[10]=2; pcoords[11]=0;
  pcoords[12]=2; pcoords[13]=1;
  pcoords[14]=1.5; pcoords[15]=0.5;


  ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr);


  /* Create interpolator */
  ierr = DMInterpolationCreate(PETSC_COMM_WORLD,
&interpolator);CHKERRQ(ierr);
  ierr = DMInterpolationSetDim(interpolator, spaceDim);CHKERRQ(ierr);
  ierr = DMInterpolationAddPoints(interpolator, Np, pcoords);CHKERRQ(ierr);

  //

  //VecGetArray(vcoord,&interpolator->points);

  ierr = DMInterpolationSetUp(interpolator, dm,
pointsAllProcs);CHKERRQ(ierr);

  /* Check locations */
  for (c = 0; c < interpolator->n; ++c) {
    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D is in
Cell %D\n", rank, c, interpolator->cells[c]);CHKERRQ(ierr);
  }
  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr);
  ierr = VecView(interpolator->coords,
PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /* Setup Discretization */
  ierr = PetscFECreateDefault(dm, 2, 1, PETSC_TRUE, NULL, -1,
&fe);CHKERRQ(ierr);
  ierr = DMSetField(dm, 0, (PetscObject)fe);CHKERRQ(ierr);
  //  ierr = DMSetDS(dm);CHKERRQ(ierr);
  ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);

  /* Create function */
  ierr = DMGetLocalVector(dm, &lu);CHKERRQ(ierr);
  VecSet(lu,1.0);


  /* Check interpolant */
  ierr = VecCreateSeq(PETSC_COMM_SELF, interpolator->n,
&fieldVals);CHKERRQ(ierr);
  VecZeroEntries(fieldVals);
  ierr = DMInterpolationSetDof(interpolator, 1);CHKERRQ(ierr);
  ierr = DMInterpolationEvaluate(interpolator, dm, lu,
fieldVals);CHKERRQ(ierr);
  ierr = VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);


  return ierr;
}

int main(int argc, char **argv)
{

  PetscInitialize(&argc, &argv, (char *) 0, help);
  SolveProblem(1);

  printf("done \n");
  PetscFinalize();

  return 0;
}


The interpolated field at the 8 points defined by pcoords should be all 1s,
and instead I get:
Vec Object: 1 MPI processes
  type: seq
1.
0.
0.
0.5
1.

Vec Object: 1 MPI processes
  type: seq
-1.5
-1.
2.


Also including the test gmesh file:

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
2
2 1 "left"
2 2 "right"
$EndPhysicalNames
$Nodes
8
1 0 0 0
2 1 0 0
3 1 1 0
4 0 1 0
5 2 0 0
6 2 1 0
7 0.5 0.5 0
8 1.5 0.5 0
$EndNodes
$Elements
8
1 2 2 1 1 1 2 7
2 2 2 1 1 1 7 4
3 2 2 1 1 2 3 7
4 2 2 1 1 3 4 7
5 2 2 2 2 2 8 3
6 2 2 2 2 2 5 8
7 2 2 2 2 3 8 6
8 2 2 2 2 5 6 8
$EndElements
$Periodic
7
0 1 4
1
1 4
0 2 6
1
2 6
0 5 4
1
5 4
0 6 1
1
6 1
1 1 3
0
1 5 7
0
1 6 4
0
$EndPeriodic

 Thanks,
Swarnava
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190523/c8b55e96/attachment.html>