[petsc-users] Boundary condition IS

[Yann Jobic]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171207/09fcdb70/attachment.html>