<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr">Hi PETSc developers and users,<div><br></div><div>I am trying to setup a test case of DMInterpolationEvaluate. I am referring to the example provided in (<a href="https://bitbucket.org/petsc/petsc/src/master/src/snes/examples/tests/ex2.c" target="_blank">https://bitbucket.org/petsc/petsc/src/master/src/snes/examples/tests/ex2.c</a>) to code it up.</div><div><br></div><div>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:</div><div><br></div><div>Attached is the code which I wrote to test it:</div><div><br></div><div><div>static char help[] = "Test interpolation\n\</div><div> using a parallel unstructured mesh (DMPLEX) to discretize it.\n";</div><div><br></div><div>#include <petsctime.h></div><div>#include <mpi.h></div><div>#include <petsc.h></div><div>#include <petscdmda.h></div><div>#include <petscdmplex.h></div><div>#include <petscds.h></div><div><br></div><div>#undef __FUNCT__</div><div>#define __FUNCT__ "main"</div><div><br></div><div>PetscErrorCode SolveProblem(int i)</div><div>{</div><div><br></div><div>  PetscBool interpolate=0;</div><div>  char filename[]="mesh2.msh";</div><div>  //char filename[]="mesh1";</div><div>  PetscErrorCode ierr;</div><div>  DM dm;</div><div>  DM               distributedMesh = NULL;</div><div>  PetscViewer    viewer;</div><div>  DMLabel   label;</div><div>  PetscBool      simplex=PETSC_TRUE;</div><div>  PetscFE         fe;</div><div>  PetscQuadrature q;</div><div>  PetscDS        prob=NULL;</div><div>  PetscInt id = 1;</div><div>  PetscInt order;</div><div>  PetscPartitioner part;</div><div>  SNES           snes;        /* nonlinear solver */</div><div>  int rank;</div><div>  Vec       fieldvec;</div><div><br></div><div><br></div><div>  MPI_Comm_rank(MPI_COMM_WORLD,&rank);</div><div><br></div><div>  // read mesh from file</div><div>  ierr=DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,filename,interpolate,&dm);CHKERRQ(ierr);</div><div>  </div><div>  //    Distribute mesh over processes</div><div>  ierr = DMPlexGetPartitioner(dm,&part);CHKERRQ(ierr);</div><div>  ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);</div><div>  ierr = DMPlexDistribute(dm, 0, NULL, &distributedMesh);CHKERRQ(ierr);</div><div> </div><div>  if (distributedMesh) {</div><div>   ierr = DMDestroy(&dm);CHKERRQ(ierr);</div><div>    dm  = distributedMesh;</div><div>  }</div><div><br></div><div><br></div><div>    </div><div>  Vec vcoord;</div><div>  DMGetCoordinatesLocal(dm,&vcoord);</div><div>  VecView(vcoord,PETSC_VIEWER_STDOUT_SELF);</div><div><br></div><div><span style="white-space:pre-wrap">        </span>  </div><div>  // interpolation</div><div>  PetscInt            spaceDim, c, Np=8, p;</div><div>  DMInterpolationInfo interpolator;</div><div>  PetscReal          pcoords[16];</div><div>  PetscBool           pointsAllProcs=PETSC_TRUE;</div><div>  Vec                 lu, fieldVals;</div><div>  pcoords[0]=0.0; pcoords[1]=0; </div><div>  pcoords[2]=1; pcoords[3]=0;</div><div>  pcoords[4]=1; pcoords[5]=1;</div><div>  pcoords[6]=0; pcoords[7]=1;</div><div>  pcoords[8]=0.5; pcoords[9]=0.5;</div><div>  pcoords[10]=2; pcoords[11]=0;</div><div>  pcoords[12]=2; pcoords[13]=1;</div><div>  pcoords[14]=1.5; pcoords[15]=0.5;</div><div><br></div><div><span style="white-space:pre-wrap">   </span>      </div><div>  ierr = DMGetCoordinateDim(dm, &spaceDim);CHKERRQ(ierr);</div><div><span style="white-space:pre-wrap">  </span> </div><div><span style="white-space:pre-wrap">       </span>  </div><div>  /* Create interpolator */</div><div>  ierr = DMInterpolationCreate(PETSC_COMM_WORLD, &interpolator);CHKERRQ(ierr);</div><div>  ierr = DMInterpolationSetDim(interpolator, spaceDim);CHKERRQ(ierr);</div><div>  ierr = DMInterpolationAddPoints(interpolator, Np, pcoords);CHKERRQ(ierr);</div><div><br></div><div>  //</div><div><br></div><div>  //VecGetArray(vcoord,&interpolator->points);</div><div><br></div><div>  ierr = DMInterpolationSetUp(interpolator, dm, pointsAllProcs);CHKERRQ(ierr);</div><div>  </div><div>  /* Check locations */</div><div>  for (c = 0; c < interpolator->n; ++c) {</div><div>    ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]Point %D is in Cell %D\n", rank, c, interpolator->cells[c]);CHKERRQ(ierr);</div><div>  }</div><div>  ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);CHKERRQ(ierr);</div><div>  ierr = VecView(interpolator->coords, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</div><div><br></div><div>  /* Setup Discretization */</div><div>  ierr = PetscFECreateDefault(dm, 2, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr);</div><div>  ierr = DMSetField(dm, 0, (PetscObject)fe);CHKERRQ(ierr);</div><div>  //  ierr = DMSetDS(dm);CHKERRQ(ierr);</div><div>  ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);</div><div><br></div><div>  /* Create function */</div><div>  ierr = DMGetLocalVector(dm, &lu);CHKERRQ(ierr);</div><div>  VecSet(lu,1.0);</div><div> </div><div><br></div><div>  /* Check interpolant */</div><div>  ierr = VecCreateSeq(PETSC_COMM_SELF, interpolator->n, &fieldVals);CHKERRQ(ierr);</div><div>  VecZeroEntries(fieldVals);</div><div>  ierr = DMInterpolationSetDof(interpolator, 1);CHKERRQ(ierr);</div><div>  ierr = DMInterpolationEvaluate(interpolator, dm, lu, fieldVals);CHKERRQ(ierr);</div><div>  ierr = VecView(fieldVals, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);</div><div><br></div><div> </div><div>  return ierr;</div><div>}</div><div><br></div><div>int main(int argc, char **argv)</div><div>{</div><div><br></div><div>  PetscInitialize(&argc, &argv, (char *) 0, help);</div><div>  SolveProblem(1);</div><div>   </div><div>  printf("done \n");</div><div>  PetscFinalize();</div><div><br></div><div>  return 0;</div><div>}</div></div><div><br></div><div><br></div><div>The interpolated field at the 8 points defined by pcoords should be all 1s, and instead I get:</div><div><div>Vec Object: 1 MPI processes</div><div>  type: seq</div><div>1.</div><div>0.</div><div>0.</div><div>0.5</div><div>1.</div><div><br></div><div>Vec Object: 1 MPI processes</div><div>  type: seq</div><div>-1.5</div><div>-1.</div><div>2.</div></div><div><br></div><div><br></div><div>Also including the test gmesh file:</div><div><br></div><div><div>$MeshFormat</div><div>2.2 0 8</div><div>$EndMeshFormat</div><div>$PhysicalNames</div><div>2</div><div>2 1 "left"</div><div>2 2 "right"</div><div>$EndPhysicalNames</div><div>$Nodes</div><div>8</div><div>1 0 0 0</div><div>2 1 0 0</div><div>3 1 1 0</div><div>4 0 1 0</div><div>5 2 0 0</div><div>6 2 1 0</div><div>7 0.5 0.5 0</div><div>8 1.5 0.5 0</div><div>$EndNodes</div><div>$Elements</div><div>8</div><div>1 2 2 1 1 1 2 7</div><div>2 2 2 1 1 1 7 4</div><div>3 2 2 1 1 2 3 7</div><div>4 2 2 1 1 3 4 7</div><div>5 2 2 2 2 2 8 3</div><div>6 2 2 2 2 2 5 8</div><div>7 2 2 2 2 3 8 6</div><div>8 2 2 2 2 5 6 8</div><div>$EndElements</div><div>$Periodic</div><div>7</div><div>0 1 4</div><div>1</div><div>1 4</div><div>0 2 6</div><div>1</div><div>2 6</div><div>0 5 4</div><div>1</div><div>5 4</div><div>0 6 1</div><div>1</div><div>6 1</div><div>1 1 3</div><div>0</div><div>1 5 7</div><div>0</div><div>1 6 4</div><div>0</div><div>$EndPeriodic</div></div><div><br></div><div> Thanks,</div><div>Swarnava</div></div></div></div></div></div></div>