Re: [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.

> Yes, that is the solution I use.

Thanks – good to hear that I’m on the right track :)

> I think I ignore nonlocal indices in the l2g mapping rather than making them 
> negative. This can of course be changed.

Ok

- Morten

Re: [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, );
PetscInt offset;
PetscSectionGetOffset(globalSection,point, );
return offset < 0;
}


// Old code

bool isGhost(DM dm, PetscInt point){
const PetscInt* array;
ISLocalToGlobalMapping ltogm;
DMGetLocalToGlobalMapping(dm,);
ISLocalToGlobalMappingGetIndices(ltogm, );
return array[point]<0;
}

From: Morten Nobel-Jørgensen <m...@mek.dtu.dk<mailto:m...@mek.dtu.dk>>
Date: Monday 30 November 2015 at 17:24
To: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>
Cc: "petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>" 
<petsc-users@mcs.anl.gov<mailto:petsc-users@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, );
PetscSection defaultSection;
DMGetDefaultSection(dm, );
PetscSection coordSection;
DMGetCoordinateSection(dm, );
PetscScalar *coords;
VecGetArray(coordinates, );
DM cdm;
DMGetCoordinateDM(dm, );

// iterate (local) mesh
PetscInt cellsFrom, cellsTo;
std::string s = "";
DMPlexGetHeightStratum(dm, 0, , );
for (PetscInt i=cellsFrom;i<cellsTo;i++) {
PetscInt edgesSize;
const PetscInt *edgeIndices;
DMPlexGetConeSize(dm, i, );
DMPlexGetCone(dm, i, );
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, );
DMPlexGetCone(dm, edgePoint, );

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],);

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, );

int rank;
MPI_Comm_rank (PETSC_COMM_WORLD, );   // 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 <knep...@gmail.com<mailto:knep...@gmail.com>>
Date: Monday 30 November 2015 at 14:08
To: Morten Nobel-Jørgensen <m...@mek.dtu.dk<mailto:m...@mek.dtu.dk>>
Cc: "petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>" 
<petsc-users@mcs.anl.gov<mailto:petsc-users@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 
<m...@mek.dtu.dk<mailto:m...@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 

Re: [petsc-users] DMPlex: Ghost points after DMRefine

[Morten Nobel-Jørgensen]

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, );
PetscSection defaultSection;
DMGetDefaultSection(dm, );
PetscSection coordSection;
DMGetCoordinateSection(dm, );
PetscScalar *coords;
VecGetArray(coordinates, );
DM cdm;
DMGetCoordinateDM(dm, );

// iterate (local) mesh
PetscInt cellsFrom, cellsTo;
std::string s = "";
DMPlexGetHeightStratum(dm, 0, , );
for (PetscInt i=cellsFrom;i<cellsTo;i++) {
PetscInt edgesSize;
const PetscInt *edgeIndices;
DMPlexGetConeSize(dm, i, );
DMPlexGetCone(dm, i, );
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, );
DMPlexGetCone(dm, edgePoint, );

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],);

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, );

int rank;
MPI_Comm_rank (PETSC_COMM_WORLD, );   // 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 <knep...@gmail.com<mailto:knep...@gmail.com>>
Date: Monday 30 November 2015 at 14:08
To: Morten Nobel-Jørgensen <m...@mek.dtu.dk<mailto:m...@mek.dtu.dk>>
Cc: "petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>" 
<petsc-users@mcs.anl.gov<mailto:petsc-users@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 
<m...@mek.dtu.dk<mailto:m...@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,);
ISLocalToGlobalMappingGetIndices(ltogm, );
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,);
ISLocalToGlobalMappingGetIndices(ltogm, );
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

Re: [petsc-users] DMPlex: Ghost points after DMRefine

[Matthew Knepley]

On Mon, Nov 30, 2015 at 7:01 AM, Morten Nobel-Jørgensen 
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,);
> ISLocalToGlobalMappingGetIndices(ltogm, );
> 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,);
> ISLocalToGlobalMappingGetIndices(ltogm, );
> 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