[petsc-users] DMPlex: Ghost points after DMRefine

[Matthew Knepley]

On Tue, Dec 1, 2015 at 4:23 AM, Morten Nobel-Jørgensen <mono at mek.dtu.dk>
wrote:

> I found a solution to my problem by using the global section instead. I
> still don’t quite understand what my problem ISLocalToGlobalMapping was.
>

Yes, that is the solution I use.


> // New code
>
> bool isGhost(DM dm, PetscInt point){
>     PetscSection globalSection;
>     DMGetDefaultGlobalSection(dm, &globalSection);
>     PetscInt offset;
>     PetscSectionGetOffset(globalSection,point, &offset);
>     return offset < 0;
> }
>
>
I think I ignore nonlocal indices in the l2g mapping rather than making
them negative. This can of course be changed.

   Matt


> // Old code
>
> bool isGhost(DM dm, PetscInt point){
>     const PetscInt* array;
>     ISLocalToGlobalMapping ltogm;
>     DMGetLocalToGlobalMapping(dm,&ltogm);
>     ISLocalToGlobalMappingGetIndices(ltogm, &array);
>     return array[point]<0;
> }
>
>
> From: Morten Nobel-Jørgensen <mono at mek.dtu.dk>
> Date: Monday 30 November 2015 at 17:24
> To: Matthew Knepley <knepley at gmail.com>
>
> Cc: "petsc-users at mcs.anl.gov" <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] DMPlex: Ghost points after DMRefine
>
> Hi Matt
>
> I don’t think the problem is within Petsc - rather somewhere in my code.
> When I dump the DMPlex using DMView (ascii-info–detail) the ghost mapping
> seems to be setup correctly.
>
> Is there a better way to determine if a local point is a ghost point?
>
> The way I iterate the DMPlex is like this:
>
> void iterateDMPlex(DM dm){
>     Vec coordinates;
>     DMGetCoordinatesLocal(dm, &coordinates);
>     PetscSection defaultSection;
>     DMGetDefaultSection(dm, &defaultSection);
>     PetscSection coordSection;
>     DMGetCoordinateSection(dm, &coordSection);
>     PetscScalar *coords;
>     VecGetArray(coordinates, &coords);
>     DM cdm;
>     DMGetCoordinateDM(dm, &cdm);
>
>     // iterate (local) mesh
>     PetscInt cellsFrom, cellsTo;
>     std::string s = "";
>     DMPlexGetHeightStratum(dm, 0, &cellsFrom, &cellsTo);
>     for (PetscInt i=cellsFrom;i<cellsTo;i++) {
>         PetscInt edgesSize;
>         const PetscInt *edgeIndices;
>         DMPlexGetConeSize(dm, i, &edgesSize);
>         DMPlexGetCone(dm, i, &edgeIndices);
>         s = s + "Element: "+std::to_string(i)+"\n";
>
>         for (int edgeId = 0;edgeId <edgesSize;edgeId ++){ // ignore edge orientation
>             PetscInt points = edgeIndices[edgeId];
>             PetscInt edgePoint = edgeIndices[edgeId];
>             s = s + "\tEdge: "+std::to_string(edgePoint)+"\n";
>             PetscInt vertexSize;
>             const PetscInt *vertexIndices;
>             DMPlexGetConeSize(dm, edgePoint, &vertexSize);
>             DMPlexGetCone(dm, edgePoint, &vertexIndices);
>
>             for (int j=0;j<vertexSize;j++){
>                 s = s + "\t\tVertex: "+std::to_string(vertexIndices[j]);
>
>                 s = s + " coords ";
>                 PetscScalar* values;
>                 VecGetValuesSection(coordinates, coordSection,vertexIndices[j],&values);
>
>                 int dim = 2;
>                 for (int k=0;k<dim;k++){
>                     double coordinate = values[k];
>
>                     s = s +std::to_string(coordinate)+" ";
>                 }
>
>                 s = s + (isGhost(cdm, vertexIndices[j])?" ghost":"");
>
>                 s = s + "\n";
>
>             }
>         }
>     }
>
>     VecRestoreArray(coordinates, &coords);
>
>     int rank;
>     MPI_Comm_rank (PETSC_COMM_WORLD, &rank);   // get current process id
>
>     PetscSynchronizedPrintf(PETSC_COMM_WORLD,*”*dmplex iteration rank %d \n %s\n",rank, s.c_str());
>     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
> }
>
>
> From: Matthew Knepley <knepley at gmail.com>
> Date: Monday 30 November 2015 at 14:08
> To: Morten Nobel-Jørgensen <mono at mek.dtu.dk>
> Cc: "petsc-users at mcs.anl.gov" <petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] DMPlex: Ghost points after DMRefine
>
> On Mon, Nov 30, 2015 at 7:01 AM, Morten Nobel-Jørgensen <mono at mek.dtu.dk>
> wrote:
>
>> I have a very simple unstructured mesh composed of two triangles (four
>> vertices) with one shared edge using a DMPlex:
>>
>>  /|\
>> / | \
>> \ | /
>>  \|/
>>
>> After distributing this mesh to two processes, each process owns a
>> triangle. However one process owns tree vertices, while the last vertex is
>> owned by the other process.
>>
>> The problem occurs after uniformly refining the dm. The mesh now looks
>> like this:
>>
>>  /|\
>> /\|/\
>> \/|\/
>>  \|/
>>
>> The new center vertex is now not listed as a ghost vertex but instead
>> exists as two individual points.
>>
>> Is there any way that this new center vertex could be created as a ghost
>> vertex during refinement?
>>
>
> This could be a bug with the l2g mapping. I do not recreate it when
> refining, only the SF defining the mapping.
> Here is an experiment: do not retrieve the mapping until after the
> refinement. Do you get what you want? If so,
> I can easily fix this by destroying the map when I refine.
>
>   Thanks,
>
>     Matt
>
>
>> Kind regards,
>> Morten
>>
>> Ps. Here are some code snippets for getting global point index and test
>> of point is a ghost point:
>>
>> int localToGlobal(DM dm, PetscInt point){
>>     const PetscInt* array;
>>     ISLocalToGlobalMapping ltogm;
>>     DMGetLocalToGlobalMapping(dm,&ltogm);
>>     ISLocalToGlobalMappingGetIndices(ltogm, &array);
>>     PetscInt res = array[point];
>>     if (res < 0){ // if ghost
>>         res = -res +1;
>>     }
>>     return res;
>> }
>>
>> bool isGhost(DM dm, PetscInt point){
>>     const PetscInt* array;
>>     ISLocalToGlobalMapping ltogm;
>>     DMGetLocalToGlobalMapping(dm,&ltogm);
>>     ISLocalToGlobalMappingGetIndices(ltogm, &array);
>>     return array[point]<0;
>> }
>>
>>
>
>
> --
> 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
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151201/de215115/attachment-0001.html>