[petsc-users] DMPlex: Ghost points after DMRefine

[Morten Nobel-Jørgensen]

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


// New code

bool isGhost(DM dm, PetscInt point){
    PetscSection globalSection;
    DMGetDefaultGlobalSection(dm, &globalSection);
    PetscInt offset;
    PetscSectionGetOffset(globalSection,point, &offset);
    return offset < 0;
}


// 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<mailto:mono at mek.dtu.dk>>
Date: Monday 30 November 2015 at 17:24
To: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Cc: "petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>" <petsc-users at mcs.anl.gov<mailto: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<mailto:knepley at gmail.com>>
Date: Monday 30 November 2015 at 14:08
To: Morten Nobel-Jørgensen <mono at mek.dtu.dk<mailto:mono at mek.dtu.dk>>
Cc: "petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>" <petsc-users at mcs.anl.gov<mailto: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<mailto: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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20151201/879b4933/attachment-0001.html>