[petsc-users] coordinates of vertices of a cell in 3D DMPlex

[Swarnava Ghosh]

Dear Mark and Matt,

Thank you for your suggestions. I tried what you mentioned. It works for
one call of DMPlexVecGetClosure+DMPlexVecRestoreClosure.
However, in my code, I would need to loop over cell values, and this is
when it crashes. A vanilla code snippet is

ierr=DMGetCoordinateDM(pCgdft->dmplexloc,&pCgdft->cdm);CHKERRQ(ierr);

ierr=DMGetCoordinatesLocal(pCgdft->dmplexloc,&pCgdft->VCloc);CHKERRQ(ierr);
    ierr=DMConvert(pCgdft->cdm,DMPLEX,&pCgdft->cdmplexloc);CHKERRQ(ierr);


ierr=DMGetCoordinatesLocal(pCgdft->dmplexloc,&pCgdft->VCloc);CHKERRQ(ierr);

ierr=DMPlexVecGetClosure(pCgdft->cdmplexloc,NULL,pCgdft->VCloc,cell,&closureSize,&coef);CHKERRQ(ierr);

    // first call
    if(rank==0)
      {
    for(i=0;i<closureSize;i++)
     {
       printf("&& cell=0, first call, coef[%d]=%lf \n",i,coef[i]);
     }
      }

ierr=DMPlexVecRestoreClosure(pCgdft->cdmplexloc,NULL,pCgdft->VCloc,cell,&closureSize,&coef);CHKERRQ(ierr);

    // second call

ierr=DMPlexVecGetClosure(pCgdft->cdmplexloc,NULL,pCgdft->VCloc,cell,&closureSize,&coef);CHKERRQ(ierr);

    if(rank==0)
      {
    for(i=0;i<closureSize;i++)
     {
       printf("&& cell=0, second call, coef[%d]=%lf \n",i,coef[i]);
     }
      }

ierr=DMPlexVecRestoreClosure(pCgdft->cdmplexloc,NULL,pCgdft->VCloc,cell,&closureSize,&coef);CHKERRQ(ierr);


the first printf statement gets printed, and then it crashes with the
following error message

[44]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[44]PETSC ERROR: Object is in wrong state
[44]PETSC ERROR: Array was not checked out
[44]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[44]PETSC ERROR: Petsc Release Version 3.8.4, Mar, 24, 2018
[44]PETSC ERROR: ../lib/cgdft on a arch-linux2-c-opt named
hpc-82-22.cm.cluster by swarnava Sun Jul 26 12:26:40 2020
[44]PETSC ERROR: Configure options --prefix=/software/PETSc/3.8.4-intel
--with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90
--with-blaslapack-dir=/software/Intel/2018.1/compilers_and_libraries_2018.1.163/linux/mkl\
/lib/intel64_lin --with-debugging=no --with-shared-libraries=0
--download-metis --download-parmetis --download-superlu_dist
[44]PETSC ERROR: #1 DMRestoreWorkArray() line 1281 in
/groups/hpc-support/install/PETSc/petsc-3.8.4_intel_no-debug/src/dm/interface/dm.c
[44]PETSC ERROR: #2 DMPlexVecRestoreClosure() line 4099 in
/groups/hpc-support/install/PETSc/petsc-3.8.4_intel_no-debug/src/dm/impls/plex/plex.c
[44]PETSC ERROR: #3 CreateLocalDMPLEX() line 749 in
./src/cgdft_localdmplex.cc

Also attached is the error file.

Sincerely,
SG



On Sat, Jul 25, 2020 at 8:15 AM Matthew Knepley <knepley at gmail.com> wrote:

> On Sat, Jul 25, 2020 at 4:10 AM Swarnava Ghosh <swarnava89 at gmail.com>
> wrote:
>
>> Dear Petsc users,
>>
>> I had a trivial question about DMPlex. Suppose I have a 3D mesh of
>> tetrahedrons. I want to find out the 3D coordinates of the vertices of a
>> particular cell. What would be the function to do this?
>>
>
> Lots of good responses. You can see that it is taking some time for
> canonical patterns to emerge as the "right way" to do something. The reason
> we provide
> multiple layers of interface is that user codes rely on different
> abstractions and would like to interact with PETSc using different
> assumptions.
>
> If you just want the coordinates of each vertex in some cell with the
> vertices in a canonical ordering, you can do as Mark suggested, with a
> slight modification:
>
> DM  plex, cdm;
> Vec coordinates;
>
> ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
> ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
> ierr = DMConvert(cdm, DMPLEX, &plex);CHKERRQ(ierr);
> ierr = DMPlexVecGetClosure(plex, NULL, coordinates, cell, NULL,
> &coef);CHKERRQ(ierr);
> ....
> ierr = DMPlexVecRestoreClosure(plex, NULL, coordinates, cell, NULL,
> &coef);CHKERRQ(ierr);
>
> We get a local coordinate vector, because local vectors are guaranteed to
> store everything in the closure of anything in the Plex. Global vectors
> are non-overlapping partitions, suitable for solvers, and might not have
> some of the values.
>
>   Thanks,
>
>      Matt
>
>
>> Thank you,
>> SG
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200726/79cfd510/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: MgVNx3Ny3Nz3eta1.6297.Err_out
Type: application/octet-stream
Size: 140619 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200726/79cfd510/attachment-0001.obj>