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

[Matthew Knepley]

On Sun, Jul 26, 2020 at 3:42 PM Swarnava Ghosh <swarnava89 at gmail.com> wrote:

> 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
>

Reset the input array to NULL for each invocation of GetClosure().

  Thanks,

     Matt


> 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/>
>>
>

-- 
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/6abc9001/attachment.html>