Actual source code: meshBoundary.c

  1: #ifdef PETSC_RCS_HEADER
  2: static char vcid[] = "$Id: meshBoundary.c,v 1.19 2000/10/17 13:48:57 knepley Exp $";
  3: #endif

  5: /*
  6:      Defines the interface to the mesh boundary functions
  7: */

 9:  #include src/mesh/meshimpl.h

 11: /*------------------------------------------- Mesh Boundary Creation ------------------------------------------------*/
 12: #undef  __FUNCT__
 14: /*@
 15:   MeshBoundary1DCreateSimple - This function creates the boundary of a simple linear mesh for the mesh generator.

 17:   Input Parameter:
 18: . geomCtx     - The problem geometry information

 20:   Output Parameters in bdCtx:
 21: + numBd       - The number of closed boundaries
 22: . numVertices - The number of vertices
 23: . vertices    - The vertex coordinates
 24: . markers     - The vertex markers
 25: . numSegments - The number of segments
 26: . segments    - The segment endpoints
 27: . segMarkers  - The segment markers
 28: . numHoles    - The number of holes (particles)
 29: - holes       - The hole coordinates

 31:   Level: beginner

 33: .seealso MeshBoundary1DCreate(), MeshBoundary1DDestroy()
 34: .keywords mesh, boundary
 35: @*/
 36: int MeshBoundary1DCreateSimple(MPI_Comm comm, MeshGeometryContext *geomCtx, MeshBoundary2D *bdCtx) {
 37:   double     length       = geomCtx->size[0];
 38:   double     xStart       = geomCtx->start[0];
 39:   PetscTruth periodicMesh = geomCtx->isPeriodic[0];
 40:   double    *vertices;
 41:   int       *markers;
 42:   int        rank;
 43:   int        ierr;

 46:   if (periodicMesh == PETSC_TRUE) {
 47:     bdCtx->numBd       = 0;
 48:     bdCtx->numVertices = 0;
 49:   } else {
 50:     bdCtx->numBd       = 2;
 51:     bdCtx->numVertices = 2;
 52:   }
 53:   MPI_Comm_rank(comm, &rank);
 54:   if (rank != 0) return(0);

 56:   /* Setup nodes */
 57:   PetscMalloc((bdCtx->numVertices)*2 * sizeof(double), &vertices);
 58:   PetscMalloc((bdCtx->numVertices)   * sizeof(int),    &markers);
 59:   if (periodicMesh == PETSC_FALSE) {
 60:     vertices[0] = xStart;
 61:     vertices[1] = xStart + length;
 62:     markers[0]  = BOTTOM_BD;
 63:     markers[1]  = TOP_BD;
 64:   }

 66:   bdCtx->vertices = vertices;
 67:   bdCtx->markers  = markers;
 68:   return(0);
 69: }

 71: #undef  __FUNCT__
 73: /*@
 74:   MeshBoundary1DDestroy - This function frees the memory associated with the boundary of the initial mesh.

 76:   Input Parameters in bdCtx:
 77: + vertices   - The vertex coordinates
 78: . markers    - The vertex markers
 79: . segments   - The segment endpoints
 80: . segMarkers - The segment markers
 81: - holes      - The hole coordinates

 83:   Level: beginner

 85: .seealso MeshBoundary1DCreateSimple, MeshBoundary1DCreate
 86: .keywords mesh, boundary
 87: @*/
 88: int MeshBoundary1DDestroy(MPI_Comm comm, MeshBoundary2D *bdCtx) {
 89:   int rank;

 93:   MPI_Comm_rank(comm, &rank);
 94:   if (rank == 0) {
 95:     PetscFree(bdCtx->vertices);
 96:     PetscFree(bdCtx->markers);
 97:   }
 98:   return(0);
 99: }

101: #undef  __FUNCT__
103: /*@
104:   MeshBoundary2DCreateSimple - This function creates the boundary of a simple square mesh for the mesh generator.

106:   Input Parameter:
107: . geomCtx     - The problem geometry information

109:   Output Parameters in bdCtx:
110: + numBd       - The number of closed boundaries
111: . numVertices - The number of vertices
112: . vertices    - The vertex coordinates
113: . markers     - The vertex markers
114: . numSegments - The number of segments
115: . segments    - The segment endpoints
116: . segMarkers  - The segment markers
117: . numHoles    - The number of holes (particles)
118: - holes       - The hole coordinates

120:   Level: beginner

122: .seealso MeshBoundary2DCreate(), MeshBoundary2DDestroy()
123: .keywords mesh, boundary
124: @*/
125: int MeshBoundary2DCreateSimple(MPI_Comm comm, MeshGeometryContext *geomCtx, MeshBoundary2D *bdCtx)
126: {
127:   double     length       = geomCtx->size[0];
128:   double     width        = geomCtx->size[1];
129:   double     xStart       = geomCtx->start[0];
130:   double     yStart       = geomCtx->start[1];
131:   PetscTruth periodicMesh = geomCtx->isPeriodic[0];
132:   double    *vertices;
133:   int       *markers;
134:   int       *segments;
135:   int       *segMarkers;
136:   int        rank;
137:   int        node, seg;
138:   int        ierr;

141:   if (periodicMesh == PETSC_TRUE) {
142:     bdCtx->numBd       = 2;
143:     bdCtx->numVertices = 12;
144:     bdCtx->numSegments = 20;
145:     bdCtx->numHoles    = 0;
146:   } else {
147:     bdCtx->numBd       = 1;
148:     bdCtx->numVertices = 9;
149:     bdCtx->numSegments = 12;
150:     bdCtx->numHoles    = 0;
151:   }
152:   MPI_Comm_rank(comm, &rank);
153:   if (rank != 0) return(0);

155:   /* Setup nodes */
156:   PetscMalloc((bdCtx->numVertices)*2 * sizeof(double), &vertices);
157:   PetscMalloc((bdCtx->numVertices)   * sizeof(int),    &markers);
158:   if (periodicMesh == PETSC_TRUE) {
159:     vertices[0]  = xStart;                  vertices[1]  = yStart;
160:     vertices[2]  = xStart +     length/4.0; vertices[3]  = yStart;
161:     vertices[4]  = xStart +     length/2.0; vertices[5]  = yStart;
162:     vertices[6]  = xStart + 3.0*length/4.0; vertices[7]  = yStart;
163:     vertices[8]  = xStart +     length;     vertices[9]  = yStart + width;
164:     vertices[10] = xStart + 3.0*length/4.0; vertices[11] = yStart + width;
165:     vertices[12] = xStart +     length/2.0; vertices[13] = yStart + width;
166:     vertices[14] = xStart +     length/4.0; vertices[15] = yStart + width;
167:     vertices[16] = xStart;                  vertices[17] = yStart + width/2.0;
168:     vertices[18] = xStart +     length/4.0; vertices[19] = yStart + width/2.0;
169:     vertices[20] = xStart +     length/2.0; vertices[21] = yStart + width/2.0;
170:     vertices[22] = xStart + 3.0*length/4.0; vertices[23] = yStart + width/2.0;
171:     for(node = 0; node < 4; node++)
172:       markers[node] = BOTTOM_BD;
173:     for(node = 4; node < 8; node++)
174:       markers[node] = TOP_BD;
175:     for(node = 8; node < 12; node++)
176:       markers[node] = 0;
177:   } else {
178:     vertices[0]  = xStart;              vertices[1]  = yStart;
179:     vertices[2]  = xStart + length/2.0; vertices[3]  = yStart;
180:     vertices[4]  = xStart + length;     vertices[5]  = yStart;
181:     vertices[6]  = xStart + length;     vertices[7]  = yStart + width/2.0;
182:     vertices[8]  = xStart + length;     vertices[9]  = yStart + width;
183:     vertices[10] = xStart + length/2.0; vertices[11] = yStart + width;
184:     vertices[12] = xStart;              vertices[13] = yStart + width;
185:     vertices[14] = xStart;              vertices[15] = yStart + width/2.0;
186:     vertices[16] = xStart + length/2.0; vertices[17] = yStart + width/2.0;
187:     for(node = 0; node < 8; node++)
188:       markers[node] = OUTER_BD;
189:     markers[8]      = 0;
190:   }

192:   /* Setup segments */
193:   PetscMalloc((bdCtx->numSegments)*2 * sizeof(int), &segments);
194:   PetscMalloc((bdCtx->numSegments)   * sizeof(int), &segMarkers);
195:   if (periodicMesh == PETSC_TRUE) {
196:     segments[0]  = 0;  segments[1]  = 1;  segments[2]  = 1;  segments[3]  = 2;
197:     segments[4]  = 2;  segments[5]  = 3;  segments[6]  = 3;  segments[7]  = 0;
198:     segments[8]  = 4;  segments[9]  = 5;  segments[10] = 5;  segments[11] = 6;
199:     segments[12] = 6;  segments[13] = 7;  segments[14] = 7;  segments[15] = 4;
200:     segments[16] = 8;  segments[17] = 9;  segments[18] = 9;  segments[19] = 10;
201:     segments[20] = 10; segments[21] = 11; segments[22] = 11; segments[23] = 8;
202:     segments[24] = 0;  segments[25] = 8;  segments[26] = 4;  segments[27] = 8;
203:     segments[28] = 1;  segments[29] = 9;  segments[30] = 7;  segments[31] = 9;
204:     segments[32] = 2;  segments[33] = 10; segments[34] = 6;  segments[35] = 10;
205:     segments[36] = 3;  segments[37] = 11; segments[38] = 5;  segments[39] = 11;
206:     for(seg = 0; seg < 4; seg++)
207:       segMarkers[seg] = BOTTOM_BD;
208:     for(seg = 4; seg < 8; seg++)
209:       segMarkers[seg] = TOP_BD;
210:     for(seg = 8; seg < 20; seg++)
211:       segMarkers[seg] = 0;
212:   } else {
213:     segments[0]  = 0;  segments[1]  = 1; segments[2]  = 1;  segments[3]  = 2;
214:     segments[4]  = 2;  segments[5]  = 3; segments[6]  = 3;  segments[7]  = 4;
215:     segments[8]  = 4;  segments[9]  = 5; segments[10] = 5;  segments[11] = 6;
216:     segments[12] = 6;  segments[13] = 7; segments[14] = 7;  segments[15] = 0;
217:     segments[16] = 1;  segments[17] = 8; segments[18] = 3;  segments[19] = 8;
218:     segments[20] = 5;  segments[21] = 8; segments[22] = 7;  segments[23] = 8;
219:     for(seg = 0; seg < 8; seg++)
220:       segMarkers[seg] = OUTER_BD;
221:     for(seg = 8; seg < 12; seg++)
222:       segMarkers[seg] = 0;
223:   }

225:   bdCtx->vertices   = vertices;
226:   bdCtx->markers    = markers;
227:   bdCtx->segments   = segments;
228:   bdCtx->segMarkers = segMarkers;
229:   bdCtx->holes      = PETSC_NULL;
230:   return(0);
231: }

233: #undef  __FUNCT__
235: /*@
236:   MeshBoundary2DDestroy - This function frees the memory associated with the boundary of the initial mesh.

238:   Input Parameters in bdCtx:
239: + vertices   - The vertex coordinates
240: . markers    - The vertex markers
241: . segments   - The segment endpoints
242: . segMarkers - The segment markers
243: - holes      - The hole coordinates

245:   Level: beginner

247: .seealso MeshBoundary2DCreateSimple, MeshBoundary2DCreate
248: .keywords mesh, boundary
249: @*/
250: int MeshBoundary2DDestroy(MPI_Comm comm, MeshBoundary2D *bdCtx)
251: {
252:   int rank;

256:   MPI_Comm_rank(comm, &rank);
257:   if (rank == 0) {
258:     PetscFree(bdCtx->vertices);
259:     PetscFree(bdCtx->markers);
260:     PetscFree(bdCtx->segments);
261:     PetscFree(bdCtx->segMarkers);
262:     if (bdCtx->numHoles > 0) {
263:       PetscFree(bdCtx->holes);
264:     }
265:   }
266:   return(0);
267: }

269: /*---------------------------------------- Mesh Boundary Visualization ----------------------------------------------*/
270: #undef  __FUNCT__
272: static int MeshBoundary2DView_File(Mesh mesh, MeshBoundary2D *bdCtx, PetscViewer viewer)
273: {
274:   FILE *fd;
275:   int   i;
276:   int   ierr;

279:   PetscViewerASCIIGetPointer(viewer, &fd);
280:   PetscFPrintf(PETSC_COMM_WORLD, fd, "Mesh Boundary Object:n");
281:   if (bdCtx->numBd == 1) {
282:     PetscFPrintf(PETSC_COMM_WORLD, fd, "    %d closed boundaryn", bdCtx->numBd);
283:   } else {
284:     PetscFPrintf(PETSC_COMM_WORLD, fd, "    %d closed boundariesn", bdCtx->numBd);
285:   }
286:   PetscFPrintf(PETSC_COMM_WORLD, fd, "      Boundary vertices: %d segments: %d holes: %dn",
287:                bdCtx->numVertices, bdCtx->numSegments, bdCtx->numHoles);
288:   for(i = 0; i < bdCtx->numVertices; i++) {
289:     PetscFPrintf(PETSC_COMM_WORLD, fd, "    %d %g %g %dn", i, bdCtx->vertices[i*2], bdCtx->vertices[i*2+1], bdCtx->markers[i]);
290:   }
291:   for(i = 0; i < bdCtx->numSegments; i++) {
292:     PetscFPrintf(PETSC_COMM_WORLD, fd, "    %d %d %d %dn", i, bdCtx->segments[i*2], bdCtx->segments[i*2+1], bdCtx->segMarkers[i]);
293:   }
294:   for(i = 0; i < bdCtx->numHoles; i++) {
295:     PetscFPrintf(PETSC_COMM_WORLD, fd, "    %d %g %gn", i, bdCtx->holes[i*2], bdCtx->holes[i*2+1]);
296:   }
297:   return(0);
298: }

300: #undef  __FUNCT__
302: static int MeshBoundary2DView_Draw_Mesh(Mesh mesh, MeshBoundary2D *bdCtx, PetscDraw draw)
303: {
304:   double *vertices   = bdCtx->vertices;
305:   int    *segments   = bdCtx->segments;
306:   int    *segMarkers = bdCtx->segMarkers;
307:   int     color;
308:   int     seg;
309:   int     ierr;

312:   for(seg = 0; seg < bdCtx->numSegments; seg++) {
313:     if (segMarkers[seg]) {
314:       color = PETSC_DRAW_RED;
315:     } else {
316:       color = PETSC_DRAW_BLACK;
317:     }

319:     MeshDrawLine(mesh, draw, vertices[segments[seg*2]*2], vertices[segments[seg*2]*2+1],
320:                         vertices[segments[seg*2+1]*2], vertices[segments[seg*2+1]*2+1], color);
321: 
322:   }
323:   PetscDrawSynchronizedFlush(draw);
324:   return(0);
325: }

327: #undef  __FUNCT__
329: static int MeshBoundary2DView_Draw(Mesh mesh, MeshBoundary2D *bdCtx, PetscViewer v)
330: {
331:   PetscReal       startX, endX;
332:   PetscReal       startY, endY;
333:   PetscReal       sizeX,  sizeY;
334:   double         *vertices;
335:   PetscDraw       draw;
336:   PetscTruth      isnull;
337:   PetscDrawButton button;
338:   char            statusLine[256];
339:   int             pause;
340:   PetscReal       xc, yc, scale;
341:   int             xsize, ysize;
342:   int             vert;
343:   int             ierr;

346:   PetscViewerDrawGetDraw(v, 0, &draw);
347:   PetscDrawIsNull(draw, &isnull);
348:   if (isnull) return(0);
349:   MeshGetBoundingBox(mesh, &startX, &startY, PETSC_NULL, &endX, &endY, PETSC_NULL);
350:   vertices = bdCtx->vertices;
351:   for(vert = 0; vert < bdCtx->numVertices; vert++) {
352:     if (startX > vertices[vert*2])   startX = vertices[vert*2];
353:     if (startY > vertices[vert*2+1]) startY = vertices[vert*2+1];
354:     if (endX   < vertices[vert*2])   endX   = vertices[vert*2];
355:     if (endY   < vertices[vert*2+1]) endY   = vertices[vert*2+1];
356:   }
357:   sizeX = endX - startX;
358:   sizeY = endY - startY;
359:   if (sizeX > sizeY) {
360:     ysize  = 300;
361:     xsize  = ysize * (int) (sizeX/sizeY);
362:   } else {
363:     xsize  = 300;
364:     ysize  = xsize * (int) (sizeY/sizeX);
365:   }
366:   PetscDrawResizeWindow(draw, xsize, ysize);
367:   PetscDrawSynchronizedClear(draw);

369:   ierr  = PetscDrawSetCoordinates(draw, startX-0.02*sizeX, startY-0.02*sizeY, endX+0.02*sizeX, endY+0.02*sizeY);
370:   ierr  = MeshBoundary2DView_Draw_Mesh(mesh, bdCtx, draw);
371:   ierr  = PetscDrawGetPause(draw, &pause);
372:   if (pause >= 0) {
373:     PetscSleep(pause);
374:     return(0);
375:   }

377:   /* Allow the mesh to zoom or shrink */
378:   PetscDrawCheckResizedWindow(draw);
379:   PetscDrawSynchronizedGetMouseButton(draw, &button, &xc, &yc, 0, 0);
380:   while ((button != BUTTON_RIGHT) && (button != BUTTON_CENTER_SHIFT))
381:   {
382:     PetscDrawSynchronizedClear(draw);
383:     statusLine[0] = 0;
384:     switch (button)
385:     {
386:     case BUTTON_LEFT:
387:       scale = 0.5;
388:       break;
389:     case BUTTON_CENTER:
390:       scale = 2.0;
391:       break;
392:     default:
393:       scale = 1.0;
394:       break;
395:     }
396:     if (scale != 1.0) {
397:       startX = scale*(startX + sizeX - xc) + xc - sizeX*scale;
398:       endX   = scale*(endX   - sizeX - xc) + xc + sizeX*scale;
399:       startY = scale*(startY + sizeY - yc) + yc - sizeY*scale;
400:       endY   = scale*(endY   - sizeY - yc) + yc + sizeY*scale;
401:       sizeX *= scale;
402:       sizeY *= scale;
403:     }

405:     PetscDrawSetCoordinates(draw, startX-0.02*sizeX, startY-0.02*sizeY, endX+0.02*sizeX, endY+0.02*sizeY);
406:     PetscDrawString(draw, startX - 0.02*sizeX, startY - 0.02*sizeY, PETSC_DRAW_BLACK, statusLine);
407:     MeshBoundary2DView_Draw_Mesh(mesh, bdCtx, draw);
408:     PetscDrawCheckResizedWindow(draw);
409:     PetscDrawSynchronizedGetMouseButton(draw, &button, &xc, &yc, 0, 0);
410:   }

412:   return(0);
413: }

415: #undef  __FUNCT__
417: int MeshBoundary2DView(Mesh mesh, MeshBoundary2D *bdCtx, PetscViewer viewer)
418: {
419:   PetscTruth isascii, isdraw;
420:   int        ierr;

423:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_ASCII, &isascii);
424:   PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_DRAW,  &isdraw);
425:   if (isascii == PETSC_TRUE) {
426:     MeshBoundary2DView_File(mesh, bdCtx, viewer);
427:   } else if (isdraw == PETSC_TRUE) {
428:     MeshBoundary2DView_Draw(mesh, bdCtx, viewer);
429:   }

431:   return(0);
432: }

434: /*------------------------------------------- Mesh Boundary Functions ------------------------------------------------*/
435: #undef  __FUNCT__
437: int MeshBoundary2DSetBoundingBox(Mesh mesh, MeshBoundary2D *bdCtx)
438: {
439:   MPI_Comm comm;
440:   double   startX, startY;
441:   double   endX,   endY;
442:   double  *vertices;
443:   int      rank;
444:   int      p;
445:   int      ierr;

448:   PetscObjectGetComm((PetscObject) mesh, &comm);
449:   MPI_Comm_rank(comm, &rank);
450:   if (rank == 0) {
451:     vertices        = bdCtx->vertices;
452:     startX = endX = vertices[0];
453:     startY = endY = vertices[1];
454:     for(p = 0; p < bdCtx->numVertices; p++) {
455:       /* Ignore moving boundaries (to accomodate split particles)? This should be done better */
456:       if (bdCtx->markers[p] < 0) continue;
457:       if (startX > vertices[p*2])   startX = vertices[p*2];
458:       if (startY > vertices[p*2+1]) startY = vertices[p*2+1];
459:       if (endX   < vertices[p*2])   endX   = vertices[p*2];
460:       if (endY   < vertices[p*2+1]) endY   = vertices[p*2+1];
461:     }
462:   }
463:   MPI_Bcast(&startX, 1, MPI_DOUBLE, 0, comm);
464:   MPI_Bcast(&endX,   1, MPI_DOUBLE, 0, comm);
465:   MPI_Bcast(&startY, 1, MPI_DOUBLE, 0, comm);
466:   MPI_Bcast(&endY,   1, MPI_DOUBLE, 0, comm);
467:   MeshSetBoundingBox(mesh, startX, startY, 0.0, endX, endY, 0.0);
468:   return(0);
469: }