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

[Swarnava Ghosh]

Got it! Thanks a lot.

Sincerely,
Swarnava

On Sun, Jul 26, 2020 at 1:23 PM Matthew Knepley <knepley at gmail.com> wrote:

> 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/b83b4b7a/attachment-0001.html>