Le 07/12/2017 à 08:23, Mohammad Hassan Baghaei a écrit :
Hello
I am using DMPlex to construct a fully circular 2D domain for my PDEs.
As I am discretizing the PDEs using finite difference method, I need
to know , to define the layout of PetscSection, how I can find the
boundary condition IS. I know the location of the boundaries in Sieve
chart. However, I noticed that I need to know the Index Set of that
location. Can I use the Sieve point number for Index Set. It seems I
think I am not familiar with IS. Can you help me? Thank for your help.
Amir
Hello,
In snes/examples/tutorials/ex56.c you can find the following code, which
does what you want as i understood :
260: {
261: DMLabel label;
262: IS is;
263: DMCreateLabel(dm, "boundary");
264: DMGetLabel(dm, "boundary", &label);
265: DMPlexMarkBoundaryFaces(dm, label);
266: if (run_type==0) {
267: DMGetStratumIS(dm, "boundary", 1, &is);
268: DMCreateLabel(dm,"Faces");
269: if (is) {
270: PetscInt d, f, Nf;
271: const PetscInt *faces;
272: PetscInt csize;
273: PetscSection cs;
274: Vec coordinates ;
275: DM cdm;
276: ISGetLocalSize(is, &Nf);
277: ISGetIndices(is, &faces);
278: DMGetCoordinatesLocal(dm, &coordinates);
279: DMGetCoordinateDM(dm, &cdm);
280: DMGetDefaultSection(cdm, &cs);
281: /* Check for each boundary face if any component of its
centroid is either 0.0 or 1.0 */
282: for (f = 0; f < Nf; ++f) {
283: PetscReal faceCoord;
284: PetscInt b,v;
285: PetscScalar *coords = NULL;
286: PetscInt Nv;
287: DMPlexVecGetClosure(cdm, cs, coordinates, faces[f],
&csize, &coords);
288: Nv = csize/dim; /* Calculate mean coordinate vector */
289: for (d = 0; d < dim; ++d) {
290: faceCoord = 0.0;
291: for (v = 0; v < Nv; ++v) faceCoord +=
PetscRealPart(coords[v*dim+d]);
292: faceCoord /= Nv;
293: for (b = 0; b < 2; ++b) {
294: if (PetscAbs(faceCoord - b) < PETSC_SMALL) { /*
domain have not been set yet, still [0,1]^3 */
295: DMSetLabelValue(dm, "Faces", faces[f], d*2+b+1);
296: }
297: }
298: }
299: DMPlexVecRestoreClosure(cdm, cs, coordinates, faces[f],
&csize, &coords);
300: }
301: ISRestoreIndices(is, &faces);
302: }
303: ISDestroy(&is);
304: DMGetLabel(dm, "Faces", &label);
305: DMPlexLabelComplete(dm, label);
306: }
307: }
The idea is to mark boundary faces with DMPlexMarkBoundaryFaces, and get
the IS of the labeled faces, and sort them with the position of the
centroid of the face.
Regards,
Yann
