Actual source code: tri2d.c

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

  5: /* Implements 2d triangular grids */
  6: #include "src/mesh/impls/triangular/2d/2dimpl.h"         /*I "mesh.h" I*/
  7: #include "tri2dView.h"
  8: #include "tri2dCSR.h"
  9: #ifdef PETSC_HAVE_TRIANGLE
 10: #include "tri2d_Triangle.h"
 11: #endif

 13: extern int PartitionGhostNodeIndex_Private(Partition, int, int *);
 14: EXTERN_C_BEGIN
 15: int MeshSerialize_Triangular_2D(MPI_Comm, Mesh *, PetscViewer, PetscTruth);
 16: EXTERN_C_END

 18: /*--------------------------------------------- Public Functions ----------------------------------------------------*/
 19: #undef  __FUNCT__
 21: /*@C
 22:   MeshTriangular2DCalcAreas - Calculate and store the area of each element

 24:   Collective on Mesh

 26:   Input Parameters:
 27: + mesh   - The mesh
 28: - create - The flag which creates the storage if it does not already exist

 30:   Note:
 31:   If areas have not been calculated before, and 'create' is PETSC_FALSE,
 32:   then no action taken.

 34:   Level: beginner

 36: .keywords: mesh, element, area
 37: .seealso: MeshTriangular2DCalcAspectRatios()
 38: @*/
 39: int MeshTriangular2DCalcAreas(Mesh mesh, PetscTruth create)
 40: {
 41:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
 42:   int              numFaces   = mesh->numFaces;
 43:   int              numCorners = mesh->numCorners;
 44:   int             *faces      = tri->faces;
 45:   double          *nodes      = tri->nodes;
 46:   double           x21, x31, y21, y31;
 47:   int              elem;
 48:   int              ierr;

 51:   if (tri->areas || create == PETSC_TRUE)
 52:   {
 53:     if (tri->areas) {
 54:       PetscFree(tri->areas);
 55:     }
 56:     PetscMalloc(numFaces * sizeof(double), &tri->areas);
 57:     for(elem = 0; elem < numFaces; elem++)
 58:     {
 59:       if (mesh->isPeriodic == PETSC_TRUE) {
 60:         x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2]);
 61:         y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]);
 62:         x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2]);
 63:         y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]);
 64:       } else {
 65:         x21 = nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2];
 66:         y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1];
 67:         x31 = nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2];
 68:         y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1];
 69:       }
 70:       tri->areas[elem] = PetscAbsReal(x21*y31 - x31*y21)/2.0;
 71:     }
 72:     PetscLogFlops(numFaces*8);
 73:   }
 74:   return(0);
 75: }

 77: #undef  __FUNCT__
 79: /*@C
 80:   MeshTriangular2DCalcAspectRatios - Calculate and store the aspect ratio of each element

 82:   Collective on Mesh

 84:   Input Parameters:
 85: + mesh   - The mesh
 86: - create - The flag which creates the storage if it does not already exist

 88:   Note:
 89:   If aspect ratios have not been calculated before, and 'create' is PETSC_FALSE,
 90:   then no action taken.

 92:   Level: beginner

 94: .keywords: mesh, element, aspect ratio
 95: .seealso: MeshTriangular2DCalcAreas()
 96: @*/
 97: int MeshTriangular2DCalcAspectRatios(Mesh mesh, PetscTruth create)
 98: {
 99:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
100:   int              numFaces   = mesh->numFaces;
101:   int              numCorners = mesh->numCorners;
102:   int             *faces      = tri->faces;
103:   double          *nodes      = tri->nodes;
104:   double           x21, x31, x32, y21, y31, y32;
105:   double           area;
106:   int              elem;
107:   int              ierr;

110:   if (tri->aspectRatios || create == PETSC_TRUE)
111:   {
112:     if (tri->aspectRatios) {
113:       PetscFree(tri->aspectRatios);
114:     }
115:     PetscMalloc(numFaces * sizeof(double), &tri->aspectRatios);
116:     for(elem = 0; elem < numFaces; elem++)
117:     {
118:       if (mesh->isPeriodic == PETSC_TRUE) {
119:         x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2]);
120:         y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]);
121:         x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2]);
122:         y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]);
123:         x32 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners+1]*2]);
124:         y32 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]);
125:       } else {
126:         x21 = nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2];
127:         y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1];
128:         x31 = nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2];
129:         y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1];
130:         x32 = nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners+1]*2];
131:         y32 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1];
132:       }
133:       area = PetscAbsReal(x21*y31 - x31*y21);
134:       tri->aspectRatios[elem] = PetscMax(x32*x32 + y32*y32, PetscMax(x21*x21 + y21*y21, x31*x31 + y31*y31)) / area;
135:     }
136:     PetscLogFlops(numFaces*19);
137:   }
138:   return(0);
139: }

141: /*--------------------------------------------- Private Functions ---------------------------------------------------*/
142: #undef  __FUNCT__
144: static int MeshDestroy_Triangular_2D(Mesh mesh)
145: {
146:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
147:   int              ierr;

150:   if (tri->nodes != PETSC_NULL) {
151:     PetscFree(tri->nodes);
152:     PetscFree(tri->markers);
153:     PetscFree(tri->degrees);
154:   }
155:   if (tri->nodesOld != PETSC_NULL) {
156:     PetscFree(tri->nodesOld);
157:   }
158:   if (tri->edges != PETSC_NULL) {
159:     PetscFree(tri->edges);
160:     PetscFree(tri->edgemarkers);
161:   }
162:   if (tri->faces != PETSC_NULL) {
163:     PetscFree(tri->faces);
164:   }
165:   if (tri->neighbors != PETSC_NULL) {
166:     PetscFree(tri->neighbors);
167:   }
168:   if (tri->bdNodes != PETSC_NULL) {
169:     PetscFree(tri->bdNodes);
170:     PetscFree(tri->bdEdges);
171:     PetscFree(tri->bdMarkers);
172:     PetscFree(tri->bdBegin);
173:     PetscFree(tri->bdEdgeBegin);
174:   }
175: #if 0
176:   if (tri->bdCtx != PETSC_NULL) {
177:     MeshBoundaryDestroy(tri->bdCtx);
178:   }
179:   if (tri->bdCtxNew != PETSC_NULL) {
180:     MeshBoundaryDestroy(tri->bdCtxNew);
181:   }
182: #endif
183:   if (tri->areas != PETSC_NULL) {
184:     PetscFree(tri->areas);
185:   }
186:   if (tri->aspectRatios != PETSC_NULL) {
187:     PetscFree(tri->aspectRatios);
188:   }
189:   PetscFree(tri);
190:   if (mesh->support != PETSC_NULL) {
191:     PetscFree(mesh->support);
192:   }
193:   PartitionDestroy(mesh->part);
194:   if (mesh->nodeOrdering != PETSC_NULL) {
195:     AODestroy(mesh->nodeOrdering);
196:   }
197:   if (mesh->coarseMap != PETSC_NULL) {
198:     AODestroy(mesh->coarseMap);
199:   }

201:   return(0);
202: }

204: #undef  __FUNCT__
206: int MeshDebug_Triangular_2D(Mesh mesh, PetscTruth distributed)
207: {
208:   Mesh_Triangular         *tri = (Mesh_Triangular *) mesh->data;
209:   Partition                p   = mesh->part;
210:   Partition_Triangular_2D *q   = PETSC_NULL;
211:   int                      numBd       = mesh->numBd;
212:   int                      numElements = mesh->numFaces;
213:   int                      numCorners  = mesh->numCorners;
214:   int                      numNodes    = mesh->numNodes;
215:   int                     *elements    = tri->faces;
216:   int                     *markers     = tri->markers;
217:   int                     *bdMarkers   = tri->bdMarkers;
218:   int                      degree;
219:   int                     *support;
220:   int                     *degrees;
221:   int                     *supports;
222:   int                     *sizes;
223:   int                      numProcs, rank;
224:   PetscTruth               isMidnode;
225:   int                      proc, bd, elem, nElem, sElem, sElem2, corner, node, newNode, size;
226:   int                      ierr;

229:   MPI_Comm_size(mesh->comm, &numProcs);
230:   MPI_Comm_rank(mesh->comm, &rank);
231:   if (p != PETSC_NULL) {
232:     q = (Partition_Triangular_2D *) p->data;
233:   }

235:   /* Check basic sizes */
236:   PetscMalloc(numProcs * sizeof(int), &sizes);
237:   MPI_Allgather(&mesh->numBd, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
238:   for(proc = 0, 0; proc < numProcs-1; proc++) {
239:     if (sizes[proc] != sizes[proc+1]) {
240:       PetscPrintf(mesh->comm, "The number of boundaries is different on proc %d(%d) and proc %d(%d)n",
241:                   proc, sizes[proc], proc+1, sizes[proc+1]);
242:       1;
243:     }
244:   }
245: 
246: #if 0
247:   /* I made vertices distributed in linear meshes */
248:   MPI_Allgather(&tri->numVertices, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
249:   for(proc = 0, 0; proc < numProcs-1; proc++) {
250:     if (sizes[proc] != sizes[proc+1]) {
251:       PetscPrintf(mesh->comm, "The number of vertices is different on proc %d(%d) and proc %d(%d)n",
252:                   proc, sizes[proc], proc+1, sizes[proc+1]);
253:       1;
254:     }
255:   }
256: 
257: #endif
258:   MPI_Allgather(&mesh->numCorners, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
259:   for(proc = 0, 0; proc < numProcs-1; proc++) {
260:     if (sizes[proc] != sizes[proc+1]) {
261:       PetscPrintf(mesh->comm, "The number of face corners is different on proc %d(%d) and proc %d(%d)n",
262:                   proc, sizes[proc], proc+1, sizes[proc+1]);
263:       1;
264:     }
265:   }
266: 
267:   MPI_Allgather(&mesh->numBdEdges, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
268:   for(proc = 0, 0; proc < numProcs-1; proc++) {
269:     if (sizes[proc] != sizes[proc+1]) {
270:       PetscPrintf(mesh->comm, "The number of boundary edges is different on proc %d(%d) and proc %d(%d)n",
271:                   proc, sizes[proc], proc+1, sizes[proc+1]);
272:       1;
273:     }
274:   }
275: 
276:   if (distributed == PETSC_FALSE) {
277:     MPI_Allgather(&mesh->numNodes, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
278:     for(proc = 0, 0; proc < numProcs-1; proc++) {
279:       if (sizes[proc] != sizes[proc+1]) {
280:         PetscPrintf(mesh->comm, "The number of nodes is different on proc %d(%d) and proc %d(%d)n",
281:                     proc, sizes[proc], proc+1, sizes[proc+1]);
282:         1;
283:       }
284:     }
285: 
286:     MPI_Allgather(&mesh->numBdNodes, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
287:     for(proc = 0, 0; proc < numProcs-1; proc++) {
288:       if (sizes[proc] != sizes[proc+1]) {
289:         PetscPrintf(mesh->comm, "The number of boundary nodes is different on proc %d(%d) and proc %d(%d)n",
290:                     proc, sizes[proc], proc+1, sizes[proc+1]);
291:         1;
292:       }
293:     }
294: 
295:     MPI_Allgather(&mesh->numEdges, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
296:     for(proc = 0, 0; proc < numProcs-1; proc++) {
297:       if (sizes[proc] != sizes[proc+1]) {
298:         PetscPrintf(mesh->comm, "The number of edges is different on proc %d(%d) and proc %d(%d)n",
299:                     proc, sizes[proc], proc+1, sizes[proc+1]);
300:         1;
301:       }
302:     }
303: 
304:     MPI_Allgather(&mesh->numFaces, 1, MPI_INT, sizes, 1, MPI_INT, mesh->comm);
305:     for(proc = 0, 0; proc < numProcs-1; proc++) {
306:       if (sizes[proc] != sizes[proc+1]) {
307:         PetscPrintf(mesh->comm, "The number of faces is different on proc %d(%d) and proc %d(%d)n",
308:                     proc, sizes[proc], proc+1, sizes[proc+1]);
309:         1;
310:       }
311:     }
312: 
313:   } else {
314:     MPI_Allreduce(&mesh->numNodes, &size, 1, MPI_INT, MPI_SUM, mesh->comm);
315:     if (size != q->numNodes) {
316:       SETERRQ2(PETSC_ERR_PLIB, "The total number of nodes %d should be %dn", size, q->numNodes);
317:     }
318: #if 0
319:     MPI_Allreduce(&tri->numBdNodes, &size, 1, MPI_INT, MPI_SUM, mesh->comm);
320: #else
321:     size = mesh->numBdNodes;
322: #endif
323:     if (size != q->numBdNodes) {
324:       SETERRQ2(PETSC_ERR_PLIB, "The total number of boundary nodes %d should be %dn", size, q->numBdNodes);
325:     }
326:     MPI_Allreduce(&mesh->numEdges, &size, 1, MPI_INT, MPI_SUM, mesh->comm);
327:     if (size != q->numEdges) {
328:       SETERRQ2(PETSC_ERR_PLIB, "The total number of edges %d should be %dn", size, q->numEdges);
329:     }
330:     MPI_Allreduce(&mesh->numFaces, &size, 1, MPI_INT, MPI_SUM, mesh->comm);
331:     if (size != mesh->part->numElements) {
332:       SETERRQ2(PETSC_ERR_PLIB, "The total number of faces %d should be %dn", size, mesh->part->numElements);
333:     }
334:   }
335:   PetscFree(sizes);

337:   if (distributed == PETSC_FALSE) {
338:     /* Check markers */
339:     for(node = 0; node < numNodes; node++) {
340:       if (!markers[node]) continue;

342:       for(bd = 0; bd < numBd; bd++)
343:         if(bdMarkers[bd] == markers[node])
344:           break;
345:       if (bd == numBd) {
346:         SETERRQ2(PETSC_ERR_PLIB, "The marker %d for node %d is invalidn", markers[node], node);
347:       }
348:     }
349:     /* Check mesh connectivity */
350:     for(node = 0; node < numNodes; node++) {
351:       MeshGetNodeSupport(mesh, node, 0, &degree, &support);

353:       /* Check node degree */
354:       PetscMalloc(numProcs * sizeof(int), &degrees);
355:       MPI_Allgather(&degree, 1, MPI_INT, degrees, 1, MPI_INT, mesh->comm);
356:       for(proc = 0, 0; proc < numProcs-1; proc++)
357:         if (degrees[proc] != degrees[proc+1]) {
358:           PetscPrintf(mesh->comm, "Degree of node %d is different on proc %d(%d) and proc %d(%d)n",
359:                       node, proc, degrees[proc], proc+1, degrees[proc+1]);
360:           PetscPrintf(mesh->comm, "Node Information:n");
361:           PetscSequentialPhaseBegin(mesh->comm, 1);
362:             PetscPrintf(PETSC_COMM_SELF, "  Processor %dn", rank);
363:             if ((rank == proc) || (rank == proc+1))
364:               for(sElem = 0; sElem < degree; sElem++)
365:                 PetscPrintf(PETSC_COMM_SELF, "  Support    : %dn", support[sElem]);
366:             PetscPrintf(PETSC_COMM_SELF, "  Marker     : %dn", tri->markers[node]);
367:             PetscPrintf(PETSC_COMM_SELF, "  Location   : (%g,%g)n", tri->nodes[node*2], tri->nodes[node*2+1]);
368:             PetscPrintf(PETSC_COMM_SELF, "  In Elements:");
369:             for(elem = 0, isMidnode = PETSC_FALSE; elem < numElements; elem++)
370:               for(corner = 0; corner < numCorners; corner++)
371:                 if (elements[elem*numCorners+corner] == node) {
372:                   PetscPrintf(PETSC_COMM_SELF, " %d", elem);
373:                   if (corner >= 3)
374:                     isMidnode = PETSC_TRUE;
375:                 }
376:             PetscPrintf(PETSC_COMM_SELF, "n");
377:             if (isMidnode == PETSC_TRUE)
378:               PetscPrintf(PETSC_COMM_SELF, "  Midnoden");
379:             else
380:               PetscPrintf(PETSC_COMM_SELF, "  Vertexn");
381:           PetscSequentialPhaseEnd(mesh->comm, 1);
382:           1;
383:         }
384: 
385:       PetscFree(degrees);

387:       /* Check support elements */
388:       PetscMalloc(degree*numProcs * sizeof(int), &supports);
389:       MPI_Allgather(support, degree, MPI_INT, supports, degree, MPI_INT, mesh->comm);
390:       for(sElem = 0, 0; sElem < degree; sElem++) {
391:         nElem = supports[sElem];
392:         for(proc = 1; proc < numProcs; proc++) {
393:           for(sElem2 = 0; sElem2 < degree; sElem2++)
394:             if (supports[proc*degree+sElem2] == nElem)
395:               break;
396:           if (sElem2 == degree) {
397:             PetscPrintf(mesh->comm, "Support element %d of node %d on proc %d(%d) is not present on proc %dn",
398:                         sElem, node, 0, supports[sElem], proc);
399:             PetscPrintf(mesh->comm, "  Support of node %d on proc %d:n   ", node, proc);
400:             for(sElem2 = 0; sElem2 < degree; sElem2++)
401:               PetscPrintf(mesh->comm, " %d", supports[proc*degree+sElem2]);
402:             PetscPrintf(mesh->comm, "n");
403:             1;
404:           }
405:         }
406:       }
407: 
408:       PetscFree(supports);

410:       /* Check that node only appears inside elements in the support */
411:       for(elem = 0, 0; elem < numElements; elem++)
412:         for(corner = 0; corner < numCorners; corner++) {
413:           newNode = elements[elem*numCorners+corner];
414:           if (node == newNode) {
415:             for(sElem = 0; sElem < degree; sElem++)
416:               if (support[sElem] == elem)
417:                 break;
418:             if (sElem == degree) {
419:               PetscPrintf(mesh->comm, "Node %d found in element %d which is not present in the supportn",
420:                           node, elem);
421:               1;
422:             }
423:           }
424:         }
425: 

427:       MeshRestoreNodeSupport(mesh, node, 0, &degree, &support);
428:     }
429:   } else {
430:     /* Check markers */
431:     for(node = 0; node < q->numOverlapNodes; node++) {
432:       if (!markers[node]) continue;

434:       for(bd = 0; bd < numBd; bd++)
435:         if(bdMarkers[bd] == markers[node])
436:           break;
437:       if (bd == numBd) {
438:         SETERRQ2(PETSC_ERR_PLIB, "The marker %d for node %d is invalidn", markers[node], node);
439:       }
440:     }
441:     /* Check mesh connectivity */
442:     for(node = 0; node < q->numLocNodes; node++) {
443:       MeshGetNodeSupport(mesh, node, 0, &degree, &support);

445:       /* Check that node only appears inside elements in the support */
446:       for(elem = 0, 0; elem < p->numOverlapElements; elem++) {
447:         for(corner = 0; corner < numCorners; corner++) {
448:           newNode = elements[elem*numCorners+corner];
449:           if (node == newNode) {
450:             for(sElem = 0; sElem < degree; sElem++) {
451:               if (support[sElem] == elem) break;
452:             }
453:             if (sElem == degree) {
454:               PetscPrintf(PETSC_COMM_SELF, "[%d]Node %d found in element %d which is not present in the supportn",
455:                           p->rank, node, elem);
456:               1;
457:             }
458:           }
459:         }
460:       }
461: 

463:       MeshRestoreNodeSupport(mesh, node, 0, &degree, &support);
464:     }
465:   }
466:   return(0);
467: }

469: #undef  __FUNCT__
471: int MeshSetupSupport_Triangular_2D(Mesh mesh)
472: {
473:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
474:   int              edge, node;
475:   int              ierr;

478:   /* Calculate maximum degree of vertices */
479:   if (mesh->numNodes > 0) {
480:     PetscMalloc(mesh->numNodes * sizeof(int), &tri->degrees);
481:   }
482:   PetscMemzero(tri->degrees, mesh->numNodes * sizeof(int));
483:   for(edge = 0; edge < mesh->numEdges*2; edge++) {
484:     tri->degrees[tri->edges[edge]]++;
485:   }
486:   for(node = 0, mesh->maxDegree = 0; node < mesh->numNodes; node++) {
487:     mesh->maxDegree = PetscMax(mesh->maxDegree, tri->degrees[node]);
488:   }
489:   if (mesh->maxDegree > 0) {
490:     PetscMalloc(mesh->maxDegree * sizeof(int), &mesh->support);
491:   }
492:   mesh->supportTaken = PETSC_FALSE;
493:   return(0);
494: }

496: #undef  __FUNCT__
498: int MeshCheckBoundary_Triangular_2D(Mesh mesh) {
499:   MeshBoundary2D *bdCtx         = mesh->bdCtx;
500:   int            *markers       = bdCtx->markers;
501:   int            *segments      = bdCtx->segments;
502:   int            *segMarkers    = bdCtx->segMarkers;
503:   int             numBd         = 0;
504:   int             numBdVertices = 0;
505:   int             numBdSegments = 0;
506:   int            *bdMarkers;
507:   int            *bdBegin;
508:   int            *bdSegmentBegin;
509:   int             bd, vertex, segment, marker, rank;
510:   int             ierr;

514:   MPI_Comm_rank(mesh->comm, &rank);
515:   if (rank != 0) return(0);
516:   PetscMalloc(bdCtx->numBd     * sizeof(int), &bdMarkers);
517:   PetscMalloc((bdCtx->numBd+1) * sizeof(int), &bdBegin);
518:   PetscMalloc((bdCtx->numBd+1) * sizeof(int), &bdSegmentBegin);
519:   PetscMemzero(bdBegin,        (bdCtx->numBd+1) * sizeof(int));
520:   PetscMemzero(bdSegmentBegin, (bdCtx->numBd+1) * sizeof(int));
521:   for(vertex = 0; vertex < bdCtx->numVertices; vertex++) {
522:     if (markers[vertex] == 0) continue;
523:     numBdVertices++;
524:     /* Check for new marker */
525:     for(bd = 0; bd < numBd; bd++) {
526:       if (markers[vertex] == bdMarkers[bd]) break;
527:     }
528:     if (bd == numBd) {
529:       if (numBd >= bdCtx->numBd) SETERRQ1(PETSC_ERR_ARG_CORRUPT, "More markers present than declared: %d", bdCtx->numBd);
530:       /* Insert new marker */
531:       for(bd = 0; bd < numBd; bd++) {
532:         if (markers[vertex] < bdMarkers[bd]) break;
533:       }
534:       if (bd < numBd) {
535:         PetscMemmove(&bdMarkers[bd+1], &bdMarkers[bd], (numBd-bd) * sizeof(int));
536:         PetscMemmove(&bdBegin[bd+2],   &bdBegin[bd+1], (numBd-bd) * sizeof(int));
537:         bdBegin[bd+1]        = 0;
538:         bdSegmentBegin[bd+1] = 0;
539:       }
540:       bdMarkers[bd] = markers[vertex];
541:       numBd++;
542:       bdBegin[bd+1]++;
543:     } else {
544:       bdBegin[bd+1]++;
545:     }
546:   }
547:   for(segment = 0; segment < bdCtx->numSegments; segment++) {
548:     if (segMarkers[segment] == 0) continue;
549:     marker = segMarkers[segment];
550:     for(bd = 0; bd < numBd; bd++) {
551:       if (marker == bdMarkers[bd]) break;
552:     }
553:     if (bd == numBd) SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Invalid segment marker %d on segment %d", marker, segment);
554:     numBdSegments++;
555:     bdSegmentBegin[bd+1]++;
556:     if (markers[segments[segment*2]] != marker) {
557:       SETERRQ4(PETSC_ERR_PLIB, "Marker %d on left vertex %d does not match marker %d on its segment %d",
558:                markers[segments[segment*2]], segments[segment*2], marker, segment);
559:     }
560:     if (markers[segments[segment*2+1]] != marker) {
561:       SETERRQ4(PETSC_ERR_PLIB, "Marker %d on right vertex %d does not match marker %d on its segment %d",
562:                markers[segments[segment*2+1]], segments[segment*2+1], marker, segment);
563:     }
564:   }

566:   /* Do prefix sums to get position offsets */
567:   for(bd = 2; bd <= numBd; bd++) {
568:     bdBegin[bd]        = bdBegin[bd-1]        + bdBegin[bd];
569:     bdSegmentBegin[bd] = bdSegmentBegin[bd-1] + bdSegmentBegin[bd];
570:   }

572:   if (numBd != bdCtx->numBd) {
573:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundaries %d should be %d", numBd, bdCtx->numBd);
574:   }
575:   if (bdBegin[numBd] != numBdVertices) {
576:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary vertices %d should be %d", bdBegin[numBd], numBdVertices);
577:   }
578:   if (bdSegmentBegin[numBd] != numBdSegments) {
579:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary segments %d should be %d", bdSegmentBegin[numBd], numBdSegments);
580:   }
581:   PetscFree(bdMarkers);
582:   PetscFree(bdBegin);
583:   PetscFree(bdSegmentBegin);
584:   return(0);
585: }

587: #undef  __FUNCT__
589: int MeshSetupBoundary_Triangular_2D(Mesh mesh)
590: {
591:   Mesh_Triangular *tri         = (Mesh_Triangular *) mesh->data;
592:   int              numBd       = mesh->numBd;
593:   int              numNodes    = mesh->numNodes;
594:   int             *markers     = tri->markers;
595:   int              numEdges    = mesh->numEdges;
596:   int             *edges       = tri->edges;
597:   int             *edgemarkers = tri->edgemarkers;
598:   int              numFaces    = mesh->numFaces;
599:   int              numCorners  = mesh->numCorners;
600:   int             *faces       = tri->faces;
601:   int              newNumBd;      /* Current number of different boundary markers */
602:   int              numBdEdges;    /* Current offset into bdEdges[] */
603:   int             *bdNodes;       /* Node numbers for boundary nodes ordered by boundary */
604:   int             *bdEdges;       /* Node numbers for boundary edges ordered by boundary */
605:   int             *bdBeginOff;    /* Current offset into the bdNodes or bdEdges array */
606:   int             *bdSeen;        /* Flags for boundaries already processed */
607:   PetscTruth       closed;        /* Indicates whether a boundary has been closed */
608:   int              degree;
609:   int             *support;
610:   int              rank, bd, marker, node, nextNode, midnode, bdNode, nextBdNode, midBdNode, edge, bdEdge;
611:   int              elem, cElem, supElem, corner, supCorner, tmp;
612:   int              ierr;

615:   MPI_Comm_rank(mesh->comm, &rank);
616:   PetscLogEventBegin(MESH_SetupBoundary, mesh, 0, 0, 0);
617:   tri->areas        = PETSC_NULL;
618:   tri->aspectRatios = PETSC_NULL;
619:   mesh->numBdNodes   = 0;
620:   mesh->numBdEdges   = 0;
621:   if (numBd == 0) {
622:     tri->bdMarkers   = PETSC_NULL;
623:     tri->bdBegin     = PETSC_NULL;
624:     tri->bdEdgeBegin = PETSC_NULL;
625:     tri->bdNodes     = PETSC_NULL;
626:     tri->bdEdges     = PETSC_NULL;
627:     return(0);
628:   }
629:   PetscMalloc(numBd     * sizeof(int), &tri->bdMarkers);
630:   PetscMalloc((numBd+1) * sizeof(int), &tri->bdBegin);
631:   PetscMalloc((numBd+1) * sizeof(int), &tri->bdEdgeBegin);
632:   PetscMalloc((numBd+1) * sizeof(int), &bdBeginOff);
633:   PetscLogObjectMemory(mesh, (numBd + (numBd+1)*2) * sizeof(int));
634:   PetscMemzero(tri->bdMarkers,   numBd     * sizeof(int));
635:   PetscMemzero(tri->bdBegin,     (numBd+1) * sizeof(int));
636:   PetscMemzero(tri->bdEdgeBegin, (numBd+1) * sizeof(int));
637:   PetscMemzero(bdBeginOff,       (numBd+1) * sizeof(int));

639:   for(node = 0, newNumBd = 0; node < numNodes; node++) {
640:     /* Get number of boundary nodes and markers */
641:     if (markers[node]) {
642:       mesh->numBdNodes++;
643:       /* Check for new marker */
644:       for(bd = 0; bd < newNumBd; bd++) {
645:         if (markers[node] == tri->bdMarkers[bd]) break;
646:       }
647:       if (bd == newNumBd) {
648:         /* Insert new marker */
649:         for(bd = 0; bd < newNumBd; bd++) {
650:           if (markers[node] < tri->bdMarkers[bd]) break;
651:         }
652:         if (bd < newNumBd) {
653:           PetscMemmove(&tri->bdMarkers[bd+1], &tri->bdMarkers[bd], (newNumBd-bd) * sizeof(int));
654:           PetscMemmove(&tri->bdBegin[bd+2],   &tri->bdBegin[bd+1], (newNumBd-bd) * sizeof(int));
655:           tri->bdBegin[bd+1] = 0;
656:         }
657:         tri->bdMarkers[bd] = markers[node];
658:         newNumBd++;
659:         tri->bdBegin[bd+1]++;
660:       } else {
661:                                 tri->bdBegin[bd+1]++;
662:       }
663:     }
664:   }

666:   /* Do prefix sums to get position offsets */
667:   for(bd = 2; bd <= numBd; bd++) {
668:     tri->bdBegin[bd] = tri->bdBegin[bd-1] + tri->bdBegin[bd];
669:   }

671:   /* This is shameful -- I will write the code to look over edges soon */
672:   if (numCorners == 3) {
673:     PetscMemcpy(tri->bdEdgeBegin, tri->bdBegin, (numBd+1) * sizeof(int));
674:   } else if (numCorners == 6) {
675:     /* We know that there are an even number of nodes in every boundary */
676:     for(bd = 0; bd <= numBd; bd++) {
677:       if (tri->bdBegin[bd]%2) {
678:         SETERRQ2(PETSC_ERR_PLIB, "There are %d nodes on boundary %d (not divisible by two)",
679:                  tri->bdBegin[bd]-tri->bdBegin[bd-1], tri->bdMarkers[bd]);
680:       } else {
681:         tri->bdEdgeBegin[bd] = tri->bdBegin[bd]/2;
682:       }
683:     }
684:   } else {
685:     SETERRQ1(PETSC_ERR_SUP, "Number of local nodes %d not supported", numCorners);
686:   }

688:   /* Get number of boundary edges */
689:   for(edge = 0; edge < numEdges; edge++) {
690:     marker = edgemarkers[edge];
691:     if (marker) {
692:       mesh->numBdEdges++;
693:       MeshGetMidnodeFromEdge(mesh, edge, &midnode);
694:       if (markers[edges[edge*2]] != marker) {
695:         if ((markers[edges[edge*2+1]] == marker) &&
696:             ((markers[midnode] == markers[edges[edge*2]]) || (markers[midnode] == markers[edges[edge*2+1]]))) {
697:           /* I assume here that the generator mistakenly included an edge between two boundaries */
698:           PetscPrintf(PETSC_COMM_SELF, "[%d]Removing edge %d between boundaries %d and %d from boundaryn",
699:                       rank, edge, markers[edges[edge*2]], markers[edges[edge*2+1]]);
700:           edgemarkers[edge] = 0;
701:           markers[midnode]  = 0;
702:           continue;
703:         } else {
704:           SETERRQ5(PETSC_ERR_PLIB, "Marker %d on left node %d does not match marker %d on its edge %d, right marker is %d",
705:                    markers[edges[edge*2]], edges[edge*2], marker, edge, markers[edges[edge*2+1]]);
706:         }
707:       }
708:       if (markers[edges[edge*2+1]] != marker) {
709:         if ((markers[edges[edge*2]] == marker) &&
710:             ((markers[midnode] == markers[edges[edge*2]]) || (markers[midnode] == markers[edges[edge*2+1]]))) {
711:           /* I assume here that the generator mistakenly included an edge between two boundaries */
712:           PetscPrintf(PETSC_COMM_SELF, "[%d]Removing edge %d between boundaries %d and %d from boundaryn",
713:                       rank, edge, markers[edges[edge*2]], markers[edges[edge*2+1]]);
714:           edgemarkers[edge] = 0;
715:           markers[midnode]  = 0;
716:           continue;
717:         } else {
718:           SETERRQ5(PETSC_ERR_PLIB, "Marker %d on right node %d does not match marker %d on its edge %d, left marker is %d",
719:                    markers[edges[edge*2+1]], edges[edge*2+1], marker, edge, markers[edges[edge*2]]);
720:         }
721:       }
722:       if (markers[midnode] != marker) {
723:         SETERRQ5(PETSC_ERR_PLIB, "Marker %d on midnode %d does not match marker %d on its edge %d, left marker is %d",
724:                  markers[midnode], midnode, marker, edge, markers[edges[edge*2]]);
725:       }
726:     }
727:   }

729:   /* Check boundary information consistency */
730:   if (newNumBd != numBd) {
731:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundaries %d should be %d", newNumBd, numBd);
732:   }
733:   if (tri->bdBegin[numBd] != mesh->numBdNodes) {
734:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary nodes %d should be %d", tri->bdBegin[numBd], mesh->numBdNodes);
735:   }
736:   if (tri->bdEdgeBegin[numBd] != mesh->numBdEdges) {
737:     SETERRQ2(PETSC_ERR_PLIB, "Invalid number of boundary edges %d should be %d", tri->bdEdgeBegin[numBd], mesh->numBdEdges);
738:   }

740:   PetscMalloc(mesh->numBdNodes * sizeof(int), &bdNodes);
741:   PetscMalloc(mesh->numBdEdges * sizeof(int), &bdEdges);
742:   PetscLogObjectMemory(mesh, (mesh->numBdNodes + mesh->numBdEdges) * sizeof(int));

744:   /* Split nodes by marker */
745:   PetscMemcpy(bdBeginOff, tri->bdBegin, (numBd+1) * sizeof(int));
746:   for(node = 0; node < numNodes; node++) {
747:     for(bd = 0; bd < numBd; bd++) {
748:       if (markers[node] == tri->bdMarkers[bd]) {
749: #ifdef MESH_TRACE_BOUNDARY_CREATE
750:         PetscPrintf(PETSC_COMM_WORLD, "bd: %d bdNode[%d] = %dn", bd, bdBeginOff[bd], node);
751: #endif
752:         bdNodes[bdBeginOff[bd]++] = node;
753:       }
754:     }
755:   }
756:   for(bd = 0; bd < numBd; bd++) {
757:     if (tri->bdBegin[bd+1] != bdBeginOff[bd])
758:       SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node marker information");
759:   }

761:   /* Split edges by marker */
762:   PetscMemcpy(bdBeginOff, tri->bdEdgeBegin, (numBd+1) * sizeof(int));
763:   for(edge = 0; edge < numEdges; edge++) {
764:     for(bd = 0; bd < numBd; bd++) {
765:       if (edgemarkers[edge] == tri->bdMarkers[bd]) {
766: #ifdef MESH_TRACE_BOUNDARY_CREATE
767:         PetscPrintf(PETSC_COMM_WORLD, "bd: %d bdEdge[%d] = %dn", bd, bdBeginOff[bd], edge);
768: #endif
769:         bdEdges[bdBeginOff[bd]++] = edge;
770:       }
771:     }
772:   }
773:   for(bd = 0; bd < numBd; bd++) {
774:     if (tri->bdEdgeBegin[bd+1] != bdBeginOff[bd])
775:       SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge marker information");
776:   }
777:   PetscFree(bdBeginOff);

779:   /* Order boundary counter-clockwise */
780:   PetscMalloc(numBd * sizeof(int), &bdSeen);
781:   PetscMemzero(bdSeen, numBd * sizeof(int));
782:   /* Loop over elements */
783:   for(elem = 0; elem < numFaces; elem++) {
784:     /* Find an element with a node on the boundary */
785:     for(corner = 0; corner < 3; corner++) {
786:       if (markers[faces[elem*numCorners+corner]]) {
787:         /* Find boundary index */
788:         cElem = elem;
789:         node  = faces[elem*numCorners+corner];
790:         for(bd = 0; bd < numBd; bd++)
791:           if (markers[node] == tri->bdMarkers[bd])
792:             break;
793:         if (bd == numBd) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary marker");
794:         /* Check for processed boundaries */
795:         if (bdSeen[bd])
796:           continue;
797:         numBdEdges = tri->bdEdgeBegin[bd];

799:         /* Start building this boundary */
800: #ifdef MESH_TRACE_BOUNDARY_CREATE
801:         PetscPrintf(mesh->comm, "Starting boundary %d beginning at %d ending before %dn",
802:                     bd, tri->bdBegin[bd], tri->bdBegin[bd+1]);
803: #endif
804:         /* Find initial node */
805:         for(bdNode = tri->bdBegin[bd]; bdNode < tri->bdBegin[bd+1]; bdNode++) {
806:           if (bdNodes[bdNode] == node) {
807:               /* Move this node to the head of the list */
808: #ifdef MESH_TRACE_BOUNDARY_CREATE
809:             PetscPrintf(mesh->comm, "   moving node %d from %d to %dn", bdNodes[bdNode], bdNode, tri->bdBegin[bd]);
810: #endif
811:               tmp                       = bdNodes[tri->bdBegin[bd]];
812:               bdNodes[tri->bdBegin[bd]] = bdNodes[bdNode];
813:               bdNodes[bdNode]           = tmp;
814:               break;
815:           }
816:         }

818:         /* Order edges counterclockwise around a boundary */
819:         /* I do not currently check the orientation of the constructed boundary */
820:         for(bdNode = tri->bdBegin[bd], closed = PETSC_FALSE; bdNode < tri->bdBegin[bd+1]; bdNode++) {
821: #ifdef MESH_TRACE_BOUNDARY_CREATE
822:           PetscPrintf(mesh->comm, "n  At boundary node %d x:%lf y: %lfn",
823:                       bdNode, tri->nodes[bdNodes[bdNode]*2], tri->nodes[bdNodes[bdNode]*2+1]);
824: #ifdef MESH_TRACE_BOUNDARY_CREATE_DETAIL
825:           for(node = tri->bdBegin[bd]; node < tri->bdBegin[bd+1]; node++)
826:             PetscPrintf(mesh->comm, "    bdNode[%d]: %dn", node, bdNodes[node]);
827: #endif
828: #endif

830:           /* Find a neighbor of the point -- Could maybe do better than linear search */
831:           for(bdEdge = numBdEdges; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++) {
832:             edge = bdEdges[bdEdge];
833:             if ((edgemarkers[edge] == tri->bdMarkers[bd]) &&
834:                 (((edges[edge*2]   == bdNodes[bdNode]) && (markers[edges[edge*2+1]] == tri->bdMarkers[bd])) ||
835:                  ((edges[edge*2+1] == bdNodes[bdNode]) && (markers[edges[edge*2]]   == tri->bdMarkers[bd]))))
836:             {
837:               /* Get neighboring node number */
838:               if (edges[edge*2] == bdNodes[bdNode]) {
839:                 node = edges[edge*2+1];
840:               } else {
841:                 node = edges[edge*2];
842:               }

844:               /* Find neighboring node in bdNodes[] */
845:               for(nextBdNode = tri->bdBegin[bd]; nextBdNode < tri->bdBegin[bd+1]; nextBdNode++)
846:                 if (bdNodes[nextBdNode] == node) break;
847: #ifdef MESH_TRACE_BOUNDARY_CREATE
848:               PetscPrintf(mesh->comm, "   found connection along edge %d to bd node %d = node %dn",
849:                           edge, nextBdNode, node);
850: #endif

852:               if (nextBdNode > bdNode) {
853:                 /* Insert midnode */
854:                 if (numCorners == 6) {
855:                   /* Move this node next to the connected one */
856: #ifdef MESH_TRACE_BOUNDARY_CREATE
857:                   PetscPrintf(mesh->comm, "   moving node %d from %d to %dn", bdNodes[nextBdNode], nextBdNode, bdNode+2);
858: #endif
859:                   tmp                 = bdNodes[bdNode+2];
860:                   bdNodes[bdNode+2]   = bdNodes[nextBdNode];
861:                   bdNodes[nextBdNode] = tmp;
862:                   nextBdNode = bdNode+2;

864:                   /* Walk around the node looking for a boundary edge and reset containing element */
865:                   node     = bdNodes[bdNode];
866:                   nextNode = bdNodes[nextBdNode];
867:                   ierr     = MeshGetNodeSupport(mesh, node, cElem, &degree, &support);
868:                   for(supElem = 0, midnode = -1; supElem < degree; supElem++) {
869:                     for(supCorner = 0; supCorner < 3; supCorner++) {
870:                       cElem     = support[supElem];
871:                       if ((faces[cElem*numCorners+supCorner]         == node) &&
872:                           (faces[cElem*numCorners+((supCorner+1)%3)] == nextNode))
873:                       {
874:                         midnode = faces[cElem*numCorners+((((supCorner*2+1)%3)*2)%3 + 3)];
875:                         supElem = degree;
876:                         break;
877:                       }
878:                       else if ((faces[cElem*numCorners+supCorner]         == node) &&
879:                                (faces[cElem*numCorners+((supCorner+2)%3)] == nextNode))
880:                       {
881:                         midnode = faces[cElem*numCorners+((((supCorner*2+2)%3)*2)%3 + 3)];
882:                         supElem = degree;
883:                         break;
884:                       }
885:                     }
886:                   }
887:                   ierr     = MeshRestoreNodeSupport(mesh, node, cElem, &degree, &support);
888:                   /* Find midnode in bdNodes[] */
889:                   for(midBdNode = tri->bdBegin[bd]; midBdNode < tri->bdBegin[bd+1]; midBdNode++)
890:                     if (bdNodes[midBdNode] == midnode)
891:                       break;
892:                   if (midBdNode == tri->bdBegin[bd+1]) SETERRQ(PETSC_ERR_PLIB, "Unable to locate midnode on boundary");
893: #ifdef MESH_TRACE_BOUNDARY_CREATE
894:                   PetscPrintf(mesh->comm, "   moving midnode %d in elem %d from %d to %dn",
895:                               midnode, cElem, midBdNode, bdNode+1);
896: #endif
897:                   tmp                = bdNodes[bdNode+1];
898:                   bdNodes[bdNode+1]  = bdNodes[midBdNode];
899:                   bdNodes[midBdNode] = tmp;
900:                   bdNode++;
901:                 } else {
902:                   /* numCorners == 3 */
903:                   /* Move this node next to the connected one */
904: #ifdef MESH_TRACE_BOUNDARY_CREATE
905:                   PetscPrintf(mesh->comm, "   moving node %d from %d to %dn", bdNodes[nextBdNode], nextBdNode, bdNode+1);
906: #endif
907:                   tmp                 = bdNodes[bdNode+1];
908:                   bdNodes[bdNode+1]   = bdNodes[nextBdNode];
909:                   bdNodes[nextBdNode] = tmp;
910:                 }

912:                 /* Reorder bdEdges[] */
913: #ifdef MESH_TRACE_BOUNDARY_CREATE
914:                 PetscPrintf(mesh->comm, "   moving edge %d from %d to %dn", edge, bdEdge, numBdEdges);
915: #endif
916:                 tmp                 = bdEdges[numBdEdges];
917:                 bdEdges[numBdEdges] = bdEdges[bdEdge];
918:                 bdEdges[bdEdge]     = tmp;
919:                 numBdEdges++;
920:                 break;
921:               }
922:               /* Check that first and last nodes are connected (closed boundary) */
923:               else if ((nextBdNode == tri->bdBegin[bd]) && (bdNode != (numCorners == 6 ? nextBdNode+2 : nextBdNode+1)))
924:               {
925:                 if (bdEdges[numBdEdges++] != edge) SETERRQ(PETSC_ERR_PLIB, "Invalid closing edge");
926:                 closed = PETSC_TRUE;
927:                 /* End the loop since only a midnode is left */
928:                 bdNode = tri->bdBegin[bd+1];
929: #ifdef MESH_TRACE_BOUNDARY_CREATE
930:                 PetscPrintf(mesh->comm, "   adding edge %d total: %dn", edge, numBdEdges);
931: #endif
932:                 break;
933:               }
934:             }
935:           }
936:           if (bdEdge == tri->bdEdgeBegin[bd+1]) SETERRQ(PETSC_ERR_PLIB, "No connection");
937:         }

939:         /* Check for closed boundary */
940:         if (closed == PETSC_FALSE) {
941:           PetscPrintf(PETSC_COMM_SELF, "Boundary %d with marker %d is not closedn", bd, tri->bdMarkers[bd]);
942:         }

944:         /* Check size of boundary */
945:         if (tri->bdBegin[bd+1] != numBdEdges*(numCorners == 3 ? 1 : 2)) {
946:           PetscPrintf(PETSC_COMM_SELF, "On boundary %d with marker %d, edge count %d x %d should equal %d (was %d)n",
947:                       bd, tri->bdMarkers[bd], numBdEdges, (numCorners == 3 ? 1 : 2), tri->bdBegin[bd+1], tri->bdBegin[bd]);
948:           SETERRQ(PETSC_ERR_PLIB, "Invalid boundary edge marker information");
949:         }
950: #ifdef MESH_TRACE_BOUNDARY_CREATE
951:         /* Check boundary nodes */
952:         for(bdNode = tri->bdBegin[bd]; bdNode < tri->bdBegin[bd+1]; bdNode++)
953:           PetscPrintf(mesh->comm, "bd: %4d bdNode: %4d node: %5dn",
954:                       bd, bdNode, bdNodes[bdNode], bdEdges[bdNode], edges[bdEdges[bdNode]*2], edges[bdEdges[bdNode]*2+1]);
955:         /* Check boundary edges */
956:         for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++)
957:           PetscPrintf(mesh->comm, "bd: %4d bdEdge: %4d edge: %5d start: %5d end: %5dn",
958:                       bd, bdEdge, bdEdges[bdEdge], edges[bdEdges[bdEdge]*2], edges[bdEdges[bdEdge]*2+1]);
959: #endif
960:         bdSeen[bd] = 1;
961:       }
962:     }
963:   }
964:   PetscFree(bdSeen);

966:   /* Set fields */
967:   tri->bdNodes = bdNodes;
968:   tri->bdEdges = bdEdges;
969:   /* Diagnostics */
970: #ifdef MESH_TRACE_BOUNDARY_REORDER
971:   int minNode, minBdNode, minEdge, minBdEdge, dir, nextBdEdge;

973:   PetscMalloc(mesh->numBdNodes * sizeof(int), &bdNodes);
974:   PetscMalloc(mesh->numBdEdges * sizeof(int), &bdEdges);
975:   for(bd = 0; bd < numBd; bd++) {
976:     /* Find smallest node */
977:     for(bdNode = tri->bdBegin[bd], minNode = tri->numNodes, minBdNode = -1; bdNode < tri->bdBegin[bd+1]; bdNode++)
978:       if (tri->bdNodes[bdNode] < minNode) {
979:         minNode   = tri->bdNodes[bdNode];
980:         minBdNode = bdNode;
981:       }
982:     /* Proceed in the direction of the smaller node */
983:     if (minBdNode == tri->bdBegin[bd]) {
984:       if (tri->bdNodes[minBdNode+1] < tri->bdNodes[tri->bdBegin[bd+1]-1])
985:         dir = 1;
986:       else
987:         dir = 0;
988:     } else if (minBdNode == tri->bdBegin[bd+1]-1) {
989:       if (tri->bdNodes[tri->bdBegin[bd]] < tri->bdNodes[minBdNode-1])
990:         dir = 0;
991:       else
992:         dir = 1;
993:     } else if (tri->bdNodes[minBdNode+1] < tri->bdNodes[minBdNode-1]) {
994:       dir = 1;
995:     } else {
996:       dir = 0;
997:     }
998:     if (dir)
999:     {
1000:       for(nextBdNode = tri->bdBegin[bd], bdNode = minBdNode; bdNode < tri->bdBegin[bd+1]; nextBdNode++, bdNode++)
1001:         bdNodes[nextBdNode] = tri->bdNodes[bdNode];
1002:       for(bdNode = tri->bdBegin[bd]; bdNode < minBdNode; nextBdNode++, bdNode++)
1003:         bdNodes[nextBdNode] = tri->bdNodes[bdNode];
1004:     }
1005:     else
1006:     {
1007:       for(nextBdNode = tri->bdBegin[bd], bdNode = minBdNode; bdNode >= tri->bdBegin[bd]; nextBdNode++, bdNode--)
1008:         bdNodes[nextBdNode] = tri->bdNodes[bdNode];
1009:       for(bdNode = tri->bdBegin[bd+1]-1; bdNode > minBdNode; nextBdNode++, bdNode--)
1010:         bdNodes[nextBdNode] = tri->bdNodes[bdNode];
1011:     }
1012:   }
1013:   for(bd = 0; bd < numBd; bd++) {
1014:     /* Find smallest edge */
1015:     for(bdEdge = tri->bdEdgeBegin[bd], minEdge = tri->numEdges, minBdEdge = -1; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++)
1016:       if (tri->bdEdges[bdEdge] < minEdge) {
1017:         minEdge   = tri->bdEdges[bdEdge];
1018:         minBdEdge = bdEdge;
1019:       }
1020:     /* Proceed in the direction of the smaller edge */
1021:     if (minBdEdge == tri->bdEdgeBegin[bd]) {
1022:       if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[tri->bdEdgeBegin[bd+1]-1])
1023:         dir = 1;
1024:       else
1025:         dir = 0;
1026:     } else if (minBdEdge == tri->bdEdgeBegin[bd+1]-1) {
1027:       if (tri->bdEdges[tri->bdEdgeBegin[bd]] < tri->bdEdges[minBdEdge-1])
1028:         dir = 0;
1029:       else
1030:         dir = 1;
1031:     } else if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[minBdEdge-1]) {
1032:       dir = 1;
1033:     } else {
1034:       dir = 0;
1035:     }
1036:     if (dir)
1037:     {
1038:       for(nextBdEdge = tri->bdEdgeBegin[bd], bdEdge = minBdEdge; bdEdge < tri->bdEdgeBegin[bd+1]; nextBdEdge++, bdEdge++)
1039:         bdEdges[nextBdEdge] = tri->bdEdges[bdEdge];
1040:       for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < minBdEdge; nextBdEdge++, bdEdge++)
1041:         bdEdges[nextBdEdge] = tri->bdEdges[bdEdge];
1042:     }
1043:     else
1044:     {
1045:       for(nextBdEdge = tri->bdEdgeBegin[bd], bdEdge = minBdEdge; bdEdge >= tri->bdEdgeBegin[bd]; nextBdEdge++, bdEdge--)
1046:         bdEdges[nextBdEdge] = tri->bdEdges[bdEdge];
1047:       for(bdEdge = tri->bdEdgeBegin[bd+1]-1; bdEdge > minBdEdge; nextBdEdge++, bdEdge--)
1048:         bdEdges[nextBdEdge] = tri->bdEdges[bdEdge];
1049:     }
1050:   }
1051:   PetscMemcpy(tri->bdNodes, bdNodes, mesh->numBdNodes * sizeof(int));
1052:   PetscMemcpy(tri->bdEdges, bdEdges, mesh->numBdEdges * sizeof(int));
1053:   PetscFree(bdNodes);
1054:   PetscFree(bdEdges);
1055: #endif
1056: #ifdef MESH_TRACE_BOUNDARY
1057:   int minNode, minBdNode, minEdge, minBdEdge, dir;

1059:   for(bd = 0; bd < numBd; bd++) {
1060:     /* Find smallest node */
1061:     for(bdNode = tri->bdBegin[bd], minNode = tri->numNodes, minBdNode = -1; bdNode < tri->bdBegin[bd+1]; bdNode++)
1062:       if (tri->bdNodes[bdNode] < minNode) {
1063:         minNode   = tri->bdNodes[bdNode];
1064:         minBdNode = bdNode;
1065:       }
1066:     /* Proceed in the direction of the smaller node */
1067:     if (minBdNode == tri->bdBegin[bd]) {
1068:       if (tri->bdNodes[minBdNode+1] < tri->bdNodes[tri->bdBegin[bd+1]-1])
1069:         dir = 1;
1070:       else
1071:         dir = 0;
1072:     } else if (minBdNode == tri->bdBegin[bd+1]-1) {
1073:       if (tri->bdNodes[tri->bdBegin[bd]] < tri->bdNodes[minBdNode-1])
1074:         dir = 0;
1075:       else
1076:         dir = 1;
1077:     } else if (tri->bdNodes[minBdNode+1] < tri->bdNodes[minBdNode-1]) {
1078:       dir = 1;
1079:     } else {
1080:       dir = 0;
1081:     }
1082:     if (dir)
1083:     {
1084:       for(bdNode = minBdNode; bdNode < tri->bdBegin[bd+1]; bdNode++)
1085:         PetscPrintf(mesh->comm, "bd: %4d node: %5dn", bd, tri->bdNodes[bdNode]);
1086:       for(bdNode = tri->bdBegin[bd]; bdNode < minBdNode; bdNode++)
1087:         PetscPrintf(mesh->comm, "bd: %4d node: %5dn", bd, tri->bdNodes[bdNode]);
1088:     }
1089:     else
1090:     {
1091:       for(bdNode = minBdNode; bdNode >= tri->bdBegin[bd]; bdNode--)
1092:         PetscPrintf(mesh->comm, "bd: %4d node: %5dn", bd, tri->bdNodes[bdNode]);
1093:       for(bdNode = tri->bdBegin[bd+1]-1; bdNode > minBdNode; bdNode--)
1094:         PetscPrintf(mesh->comm, "bd: %4d node: %5dn", bd, tri->bdNodes[bdNode]);
1095:     }
1096:   }
1097:   for(bd = 0; bd < numBd; bd++) {
1098:     /* Find smallest edge */
1099:     for(bdEdge = tri->bdEdgeBegin[bd], minEdge = tri->numEdges, minBdEdge = -1; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++)
1100:       if (tri->bdEdges[bdEdge] < minEdge) {
1101:         minEdge   = tri->bdEdges[bdEdge];
1102:         minBdEdge = bdEdge;
1103:       }
1104:     /* Proceed in the direction of the smaller edge */
1105:     if (minBdEdge == tri->bdEdgeBegin[bd]) {
1106:       if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[tri->bdEdgeBegin[bd+1]-1])
1107:         dir = 1;
1108:       else
1109:         dir = 0;
1110:     } else if (minBdEdge == tri->bdEdgeBegin[bd+1]-1) {
1111:       if (tri->bdEdges[tri->bdEdgeBegin[bd]] < tri->bdEdges[minBdEdge-1])
1112:         dir = 0;
1113:       else
1114:         dir = 1;
1115:     } else if (tri->bdEdges[minBdEdge+1] < tri->bdEdges[minBdEdge-1]) {
1116:       dir = 1;
1117:     } else {
1118:       dir = 0;
1119:     }
1120:     if (dir)
1121:     {
1122:       for(bdEdge = minBdEdge; bdEdge < tri->bdEdgeBegin[bd+1]; bdEdge++)
1123:         PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5dn",
1124:                     bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]);
1125:       for(bdEdge = tri->bdEdgeBegin[bd]; bdEdge < minBdEdge; bdEdge++)
1126:         PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5dn",
1127:                     bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]);
1128:     }
1129:     else
1130:     {
1131:       for(bdEdge = minBdEdge; bdEdge >= tri->bdEdgeBegin[bd]; bdEdge--)
1132:         PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5dn",
1133:                     bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]);
1134:       for(bdEdge = tri->bdEdgeBegin[bd+1]-1; bdEdge > minBdEdge; bdEdge--)
1135:         PetscPrintf(mesh->comm, "bd: %4d edge: %5d start: %5d end: %5dn",
1136:                     bd, tri->bdEdges[bdEdge], edges[tri->bdEdges[bdEdge]*2], edges[tri->bdEdges[bdEdge]*2+1]);
1137:     }
1138:   }
1139: #endif

1141:   PetscLogEventEnd(MESH_SetupBoundary, mesh, 0, 0, 0);
1142:   return(0);
1143: }

1145: #undef  __FUNCT__
1147: /*
1148:   MeshAssemble_Private - This function assembles the mesh entirely on the first processor.

1150:   Collective on Mesh

1152:   Input Parameters:
1153: . mesh    - The mesh being refined

1155:   Output Parameter:
1156: + nodes   - The node coordinates
1157: . markers - The node markers
1158: . faces   - The nodes for each element
1159: - edges   - The nodes for each edge

1161:   Level: developer

1163: .keywords mesh, assemble
1164: .seealso MeshInitRefineInput_Triangle()
1165: */
1166: int MeshAssemble_Private(Mesh mesh, double **nodes, int **markers, int **faces, int **edges)
1167: {
1168:   Mesh_Triangular         *tri        = (Mesh_Triangular *) mesh->data;
1169:   Partition                p          = mesh->part;
1170:   Partition_Triangular_2D *q          = (Partition_Triangular_2D *) p->data;
1171:   int                      numNodes   = q->numNodes;
1172:   int                      numFaces   = p->numElements;
1173:   int                      numEdges   = q->numEdges;
1174:   int                      numCorners = mesh->numCorners;
1175:   int                      numProcs   = p->numProcs;
1176:   int                      rank       = p->rank;
1177:   int                     *locFaces   = PETSC_NULL;
1178:   int                     *locEdges   = PETSC_NULL;
1179:   int                     *numRecvNodes, *numRecvMarkers, *numRecvFaces, *numRecvEdges;
1180:   int                     *nodeOffsets,  *markerOffsets,  *faceOffsets,  *edgeOffsets;
1181:   int                      proc, elem, edge, size;
1182:   int                      ierr;

1185:   /* Allocate global arrays */
1186:   if (rank == 0) {
1187:     PetscMalloc(numNodes*2          * sizeof(double), nodes);
1188:     PetscMalloc(numNodes            * sizeof(int),    markers);
1189:     PetscMalloc(numFaces*numCorners * sizeof(int),    faces);
1190:     PetscMalloc(numEdges*2          * sizeof(int),    edges);
1191:   }

1193:   if (numProcs > 1) {
1194:     if (mesh->numFaces > 0) {
1195:       PetscMalloc(mesh->numFaces*numCorners * sizeof(int), &locFaces);
1196:     }
1197:     if (mesh->numEdges > 0) {
1198:       PetscMalloc(mesh->numEdges*2          * sizeof(int), &locEdges);
1199:     }

1201:     /* Calculate offsets */
1202:     PetscMalloc(numProcs * sizeof(int), &numRecvNodes);
1203:     PetscMalloc(numProcs * sizeof(int), &numRecvMarkers);
1204:     PetscMalloc(numProcs * sizeof(int), &numRecvEdges);
1205:     PetscMalloc(numProcs * sizeof(int), &numRecvFaces);
1206:     PetscMalloc(numProcs * sizeof(int), &nodeOffsets);
1207:     PetscMalloc(numProcs * sizeof(int), &markerOffsets);
1208:     PetscMalloc(numProcs * sizeof(int), &edgeOffsets);
1209:     PetscMalloc(numProcs * sizeof(int), &faceOffsets);
1210:     for(proc = 0; proc < numProcs; proc++) {
1211:       numRecvNodes[proc]   = (q->firstNode[proc+1]    - q->firstNode[proc])*2;
1212:       numRecvMarkers[proc] = (q->firstNode[proc+1]    - q->firstNode[proc]);
1213:       numRecvEdges[proc]   = (q->firstEdge[proc+1]    - q->firstEdge[proc])*2;
1214:       numRecvFaces[proc]   = (p->firstElement[proc+1] - p->firstElement[proc])*numCorners;
1215:     }
1216:     nodeOffsets[0]   = 0;
1217:     markerOffsets[0] = 0;
1218:     edgeOffsets[0]   = 0;
1219:     faceOffsets[0]   = 0;
1220:     for(proc = 1; proc < numProcs; proc++) {
1221:       nodeOffsets[proc]   = numRecvNodes[proc-1]   + nodeOffsets[proc-1];
1222:       markerOffsets[proc] = numRecvMarkers[proc-1] + markerOffsets[proc-1];
1223:       edgeOffsets[proc]   = numRecvEdges[proc-1]   + edgeOffsets[proc-1];
1224:       faceOffsets[proc]   = numRecvFaces[proc-1]   + faceOffsets[proc-1];
1225:     }

1227:     /* Local to global node number conversion */
1228:     for(elem = 0; elem < mesh->numFaces*numCorners; elem++) {
1229:       PartitionLocalToGlobalNodeIndex(p, tri->faces[elem], &locFaces[elem]);
1230:     }
1231:     for(edge = 0; edge < mesh->numEdges*2; edge++) {
1232:       PartitionLocalToGlobalNodeIndex(p, tri->edges[edge], &locEdges[edge]);
1233:     }

1235:     /* Collect global arrays */
1236:     size = mesh->numNodes*2;
1237:     MPI_Gatherv(tri->nodes,   size, MPI_DOUBLE, *nodes,   numRecvNodes,   nodeOffsets,   MPI_DOUBLE, 0, p->comm);
1238: 
1239:     size = mesh->numNodes;
1240:     MPI_Gatherv(tri->markers, size, MPI_INT,    *markers, numRecvMarkers, markerOffsets, MPI_INT,    0, p->comm);
1241: 
1242:     size = mesh->numEdges*2;
1243:     MPI_Gatherv(locEdges,     size, MPI_INT,    *edges,   numRecvEdges,   edgeOffsets,   MPI_INT,    0, p->comm);
1244: 
1245:     size = mesh->numFaces*numCorners;
1246:     MPI_Gatherv(locFaces,     size, MPI_INT,    *faces,   numRecvFaces,   faceOffsets,   MPI_INT,    0, p->comm);
1247: 

1249:     /* Cleanup */
1250:     if (locFaces != PETSC_NULL) {
1251:       PetscFree(locFaces);
1252:     }
1253:     if (locEdges != PETSC_NULL) {
1254:       PetscFree(locEdges);
1255:     }
1256:     PetscFree(numRecvNodes);
1257:     PetscFree(numRecvMarkers);
1258:     PetscFree(numRecvEdges);
1259:     PetscFree(numRecvFaces);
1260:     PetscFree(nodeOffsets);
1261:     PetscFree(markerOffsets);
1262:     PetscFree(edgeOffsets);
1263:     PetscFree(faceOffsets);
1264:   } else {
1265:     /* Uniprocessor case */
1266:     PetscMemcpy(*nodes,   tri->nodes,   numNodes*2          * sizeof(double));
1267:     PetscMemcpy(*markers, tri->markers, numNodes            * sizeof(int));
1268:     PetscMemcpy(*faces,   tri->faces,   numFaces*numCorners * sizeof(int));
1269:     PetscMemcpy(*edges,   tri->edges,   numEdges*2          * sizeof(int));
1270:   }
1271:   return(0);
1272: }

1274: #undef  __FUNCT__
1276: /*
1277:   MeshCoarsen_Triangular_2D - Coarsens a two dimensional unstructured mesh using area constraints

1279:   Collective on Mesh

1281:   Input Parameters:
1282: + mesh    - The mesh begin coarsened
1283: - area    - A function which gives an area constraint when evaluated inside an element

1285:   Output Parameters:
1286: . newmesh - The coarse mesh

1288:   Note:
1289:   If PETSC_NULL is given as the 'area' argument, the mesh consisting only of vertices is returned.

1291:   Level: developer

1293: .keywords unstructured grid
1294: .seealso GridCoarsen(), MeshRefine(), MeshReform()
1295: */
1296: static int MeshCoarsen_Triangular_2D(Mesh mesh, PointFunction area, Mesh *newmesh)
1297: {
1298:   ParameterDict            dict;
1299:   Mesh                     m;
1300:   Mesh_Triangular         *tri;
1301:   Partition                p;
1302:   Partition_Triangular_2D *q;
1303:   Mesh_Triangular         *triOld = (Mesh_Triangular *) mesh->data;
1304:   Partition                pOld   = mesh->part;
1305:   Partition_Triangular_2D *qOld   = (Partition_Triangular_2D *) pOld->data;
1306:   int                    (*refine)(Mesh, PointFunction, Mesh);
1307:   int                     *newNodes;      /* newNodes[fine node #] = coarse node # or -1 */
1308:   int                     *FineMapping;   /* FineMapping[n]   = node n in the fine   mesh numbering */
1309:   int                     *CoarseMapping; /* CoarseMapping[n] = node n in the coarse mesh numbering */
1310:   PetscTruth               hasIndex;
1311:   int                     *nodePerm, *temp;
1312:   int                      numGhostElements, numGhostNodes;
1313:   int                      proc, bd, elem, edge, corner, node, newNode, gNode, ghostNode, bdNode, count;
1314:   int                      ierr;

1317:   if (area == PETSC_NULL) {
1318:     if (mesh->numCorners == 3) {
1319:       /* Return the mesh and increment the reference count instead of copying */
1320:       PetscObjectReference((PetscObject) mesh);
1321:       *newmesh = mesh;
1322:       return(0);
1323:     } else if (mesh->numCorners != 6) {
1324:       SETERRQ1(PETSC_ERR_SUP, "Number of local nodes %d not supported", mesh->numCorners);
1325:     }

1327:     /* We would like to preserve the partition of vertices, although it may be somewhat
1328:        unbalanced, since this makes mapping from the coarse mesh back to the original
1329:        much easier.
1330:     */
1331:     MeshCreate(mesh->comm, &m);
1332:     ParameterDictCreate(mesh->comm, &dict);
1333:     ParameterDictSetObject(dict, "bdCtx", mesh->bdCtx);
1334:     ParameterDictSetInteger(dict, "dim", 2);
1335:     ParameterDictSetInteger(dict, "numLocNodes", mesh->numCorners);
1336:     PetscObjectSetParameterDict((PetscObject) m, dict);
1337:     ParameterDictDestroy(dict);
1338:     PetscNew(Mesh_Triangular, &tri);
1339:     PetscLogObjectMemory(m, sizeof(Mesh_Triangular));
1340:     PetscMemcpy(m->ops, mesh->ops, sizeof(struct _MeshOps));
1341:     PetscStrallocpy(mesh->type_name,      &m->type_name);
1342:     PetscStrallocpy(mesh->serialize_name, &m->serialize_name);
1343:     m->data             = (void *) tri;
1344:     m->dim              = 2;
1345:     m->startX           = mesh->startX;
1346:     m->endX             = mesh->endX;
1347:     m->sizeX            = mesh->sizeX;
1348:     m->startY           = mesh->startY;
1349:     m->endY             = mesh->endY;
1350:     m->sizeY            = mesh->sizeY;
1351:     m->isPeriodic       = mesh->isPeriodic;
1352:     m->isPeriodicDim[0] = mesh->isPeriodicDim[0];
1353:     m->isPeriodicDim[1] = mesh->isPeriodicDim[1];
1354:     m->isPeriodicDim[2] = mesh->isPeriodicDim[2];
1355:     m->nodeOrdering     = PETSC_NULL;
1356:     m->partitioned      = 1;
1357:     m->highlightElement = mesh->highlightElement;
1358:     m->usr              = mesh->usr;

1360:     /* Copy function list */
1361:     PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", (void (**)(void)) &refine);
1362:     if (refine != PETSC_NULL) {
1363: #ifdef PETSC_HAVE_TRIANGLE
1364:       PetscObjectComposeFunction((PetscObject) m, "MeshTriangular2D_Refine_C", "MeshRefine_Triangle",
1365:                                         (void (*)(void)) MeshRefine_Triangle);
1366: #else
1367:       /* The query needs to return the name also */
1368:       PetscObjectComposeFunction((PetscObject) m, "MeshTriangular2D_Refine_C", "", (void (*)(void)) refine);
1369: #endif
1370: 
1371:     }

1373:     /* Create the partition object */
1374:     PartitionCreate(m, &p);
1375:     PetscNew(Partition_Triangular_2D, &q);
1376:     PetscLogObjectMemory(p, sizeof(Partition_Triangular_2D));
1377:     PetscMemcpy(p->ops, pOld->ops, sizeof(struct _PartitionOps));
1378:     PetscObjectChangeTypeName((PetscObject) p, pOld->type_name);
1379:     PetscObjectChangeSerializeName((PetscObject) p, pOld->serialize_name);
1380:     p->data = (void *) q;
1381:     m->part = p;
1382:     PetscLogObjectParent(m, p);

1384:     /* Construct element partition */
1385:     p->numProcs            = pOld->numProcs;
1386:     p->rank                = pOld->rank;
1387:     p->ordering            = PETSC_NULL;
1388:     p->numLocElements      = pOld->numLocElements;
1389:     p->numElements         = pOld->numElements;
1390:     p->numOverlapElements  = pOld->numOverlapElements;
1391:     PetscMalloc((p->numProcs+1) * sizeof(int), &p->firstElement);
1392:     PetscMemcpy(p->firstElement, pOld->firstElement, (p->numProcs+1) * sizeof(int));
1393:     p->ghostElements       = PETSC_NULL;
1394:     p->ghostElementProcs   = PETSC_NULL;
1395:     numGhostElements       = pOld->numOverlapElements - pOld->numLocElements;
1396:     if (numGhostElements > 0) {
1397:       PetscMalloc(numGhostElements * sizeof(int), &p->ghostElements);
1398:       PetscMalloc(numGhostElements * sizeof(int), &p->ghostElementProcs);
1399:       PetscLogObjectMemory(p, numGhostElements*2 * sizeof(int));
1400:       PetscMemcpy(p->ghostElements,     pOld->ghostElements,     numGhostElements * sizeof(int));
1401:       PetscMemcpy(p->ghostElementProcs, pOld->ghostElementProcs, numGhostElements * sizeof(int));
1402:     }
1403:     q->nodeOrdering = PETSC_NULL;
1404:     q->edgeOrdering = PETSC_NULL;
1405:     q->numLocEdges  = qOld->numLocEdges;
1406:     q->numEdges     = qOld->numEdges;
1407:     PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstEdge);
1408:     PetscMemcpy(q->firstEdge,    qOld->firstEdge,    (p->numProcs+1) * sizeof(int));
1409:     PetscLogObjectMemory(p, (p->numProcs+1)*2 * sizeof(int));

1411:     /* We must handle ghost members since we manually construct the partition */
1412:     m->numBd       = mesh->numBd;
1413:     m->numFaces    = mesh->numFaces;
1414:     m->numCorners  = 3;
1415:     m->numEdges    = mesh->numEdges;
1416:     PetscMalloc(p->numOverlapElements*m->numCorners   * sizeof(int), &tri->faces);
1417:     PetscMalloc(p->numOverlapElements*3               * sizeof(int), &tri->neighbors);
1418:     PetscMalloc(m->numEdges*2                         * sizeof(int), &tri->edges);
1419:     PetscMalloc(m->numEdges                           * sizeof(int), &tri->edgemarkers);
1420:     PetscMalloc(qOld->numOverlapNodes                 * sizeof(int), &newNodes);
1421:     PetscLogObjectMemory(m, (p->numOverlapElements*(m->numCorners + 3) + m->numEdges*3) * sizeof(int));
1422:     for(node = 0; node < qOld->numOverlapNodes; node++)
1423:       newNodes[node] = -1;
1424:     /* Renumber the vertices and construct the face list */
1425:     for(elem = 0, newNode = 0, numGhostNodes = 0; elem < p->numOverlapElements; elem++) {
1426:       for(corner = 0; corner < 3; corner++) {
1427:         node = triOld->faces[elem*mesh->numCorners+corner];
1428:         if (newNodes[node] == -1) {
1429:           if (node >= mesh->numNodes) {
1430:             /* Mark them as ghost nodes */
1431:             numGhostNodes++;
1432:             newNodes[node] = -2;
1433:           } else {
1434:             newNodes[node] = newNode++;
1435:           }
1436:         }
1437:         tri->faces[elem*m->numCorners+corner] = node;
1438:       }
1439:     }
1440:     PetscMemcpy(tri->neighbors,   triOld->neighbors,   p->numOverlapElements*3 * sizeof(int));
1441:     PetscMemcpy(tri->edges,       triOld->edges,       m->numEdges*2           * sizeof(int));
1442:     PetscMemcpy(tri->edgemarkers, triOld->edgemarkers, m->numEdges             * sizeof(int));

1444:     /* We may have extra ghost nodes from the edges */
1445:     for(edge = 0; edge < q->numLocEdges; edge++) {
1446:       for(corner = 0; corner < 2; corner++) {
1447:         node = tri->edges[edge*2+corner];
1448:         if (newNodes[node] == -1) {
1449:           if (node >= mesh->numNodes) {
1450:             /* Mark them as ghost nodes */
1451:             numGhostNodes++;
1452:             newNodes[node] = -2;
1453:           }
1454:         }
1455:       }
1456:     }

1458:     /* Compute vertex sizes */
1459:     m->numNodes      = newNode;
1460:     m->numVertices   = m->numNodes;
1461:     MPI_Allreduce(&m->numNodes, &q->numNodes, 1, MPI_INT, MPI_SUM, m->comm);
1462:     q->numLocNodes     = newNode;
1463:     q->numOverlapNodes = newNode + numGhostNodes;

1465:     /* Compute vertex partition */
1466:     PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstNode);
1467:     PetscLogObjectMemory(p, (p->numProcs+1) * sizeof(int));
1468:     MPI_Allgather(&q->numLocNodes, 1, MPI_INT, &q->firstNode[1], 1, MPI_INT, p->comm);
1469:     for(proc = 1, q->firstNode[0] = 0; proc <= p->numProcs; proc++) {
1470:       q->firstNode[proc] += q->firstNode[proc-1];
1471:     }
1472:     q->ghostNodes       = PETSC_NULL;
1473:     q->ghostNodeProcs   = PETSC_NULL;
1474:     if (numGhostNodes > 0) {
1475:       PetscMalloc(numGhostNodes * sizeof(int), &q->ghostNodes);
1476:       PetscMalloc(numGhostNodes * sizeof(int), &q->ghostNodeProcs);
1477:       PetscLogObjectMemory(p, numGhostNodes*2 * sizeof(int));
1478:     }

1480:     /* Create the mapping from coarse nodes to vertices in the original mesh */
1481:     PetscMalloc(m->numNodes * sizeof(int), &FineMapping);
1482:     PetscMalloc(m->numNodes * sizeof(int), &CoarseMapping);
1483:     for(node = 0, count = 0; node < mesh->numNodes; node++) {
1484:       newNode = newNodes[node];
1485:       if (newNode >= 0) {
1486:         /* These are all interior nodes */
1487:         PartitionLocalToGlobalNodeIndex(pOld, node,    &FineMapping[count]);
1488:         PartitionLocalToGlobalNodeIndex(p,    newNode, &CoarseMapping[count]);
1489:         count++;
1490:       }
1491:     }
1492:     if (count != m->numNodes) SETERRQ2(PETSC_ERR_PLIB, "Invalid coarse map size %d should be %d", count, m->numNodes);
1493:     AOCreateMapping(m->comm, m->numNodes, FineMapping, CoarseMapping, &m->coarseMap);
1494:     PetscFree(FineMapping);
1495:     PetscFree(CoarseMapping);

1497:     /* Setup ghost nodes */
1498:     for(ghostNode = mesh->numNodes, count = 0; ghostNode < qOld->numOverlapNodes; ghostNode++) {
1499:       if (newNodes[ghostNode] == -2) {
1500:         PartitionLocalToGlobalNodeIndex(pOld, ghostNode, &gNode);
1501:         AOApplicationToPetsc(m->coarseMap, 1, &gNode);
1502:         q->ghostNodes[count]     = gNode;
1503:         q->ghostNodeProcs[count] = qOld->ghostNodeProcs[ghostNode-mesh->numNodes];
1504:         count++;
1505:       }
1506:     }
1507:     if (count != numGhostNodes) SETERRQ2(PETSC_ERR_PLIB, "Invalid number of ghost nodes %d should be %d", count, numGhostNodes);
1508:     /* Cleanup */
1509:     PetscFree(newNodes);

1511:     /* Resort ghost nodes */
1512:     if (numGhostNodes > 0) {
1513:       PetscMalloc(numGhostNodes*2 * sizeof(int), &nodePerm);
1514:       temp = nodePerm + numGhostNodes;
1515:       for(node = 0; node < numGhostNodes; node++)
1516:         nodePerm[node] = node;
1517:       PetscSortIntWithPermutation(numGhostNodes, q->ghostNodes, nodePerm);
1518:       for(node = 0; node < numGhostNodes; node++)
1519:         temp[node] = q->ghostNodes[nodePerm[node]];
1520:       for(node = 0; node < numGhostNodes; node++)
1521:         q->ghostNodes[node] = temp[node];
1522:       for(node = 0; node < numGhostNodes; node++)
1523:         temp[node] = q->ghostNodeProcs[nodePerm[node]];
1524:       for(node = 0; node < numGhostNodes; node++)
1525:         q->ghostNodeProcs[node] = temp[node];
1526:       PetscFree(nodePerm);
1527:     }

1529:     /* Copy vertex information */
1530:     PetscMalloc(q->numOverlapNodes*2 * sizeof(double), &tri->nodes);
1531:     PetscMalloc(q->numOverlapNodes   * sizeof(int),    &tri->markers);
1532:     PetscLogObjectMemory(m, q->numOverlapNodes*2 * sizeof(double));
1533:     PetscLogObjectMemory(m, q->numOverlapNodes   * sizeof(int));
1534:     for(newNode = 0; newNode < q->numOverlapNodes; newNode++) {
1535:       PartitionLocalToGlobalNodeIndex(p, newNode, &gNode);
1536:       AOPetscToApplication(m->coarseMap, 1, &gNode);
1537:       PartitionGlobalToLocalNodeIndex(pOld, gNode, &node);
1538:       tri->nodes[newNode*2]   = triOld->nodes[node*2];
1539:       tri->nodes[newNode*2+1] = triOld->nodes[node*2+1];
1540:       tri->markers[newNode]   = triOld->markers[node];
1541:     }

1543:     /* Renumber nodes */
1544:     for(elem = 0; elem < p->numOverlapElements*m->numCorners; elem++) {
1545:       PartitionLocalToGlobalNodeIndex(pOld, tri->faces[elem], &tri->faces[elem]);
1546:     }
1547:     for(edge = 0; edge < m->numEdges*2; edge++) {
1548:       PartitionLocalToGlobalNodeIndex(pOld, tri->edges[edge], &tri->edges[edge]);
1549:     }
1550:     AOApplicationToPetsc(m->coarseMap, p->numOverlapElements*m->numCorners, tri->faces);
1551:     AOApplicationToPetsc(m->coarseMap, m->numEdges*2,                       tri->edges);
1552:     for(elem = 0; elem < p->numOverlapElements*m->numCorners; elem++) {
1553:       PartitionGlobalToLocalNodeIndex(p, tri->faces[elem], &tri->faces[elem]);
1554:     }
1555:     for(edge = 0; edge < m->numEdges*2; edge++) {
1556:       PartitionGlobalToLocalNodeIndex(p, tri->edges[edge], &tri->edges[edge]);
1557:     }

1559:     /* Construct derived and boundary information */
1560:     tri->areas        = PETSC_NULL;
1561:     tri->aspectRatios = PETSC_NULL;
1562:     m->numBdEdges     = mesh->numBdEdges;
1563:     PetscMalloc(m->numBd      * sizeof(int), &tri->bdMarkers);
1564:     PetscMalloc((m->numBd+1)  * sizeof(int), &tri->bdBegin);
1565:     PetscMalloc((m->numBd+1)  * sizeof(int), &tri->bdEdgeBegin);
1566:     PetscMalloc(m->numBdEdges * sizeof(int), &tri->bdEdges);
1567:     PetscLogObjectMemory(m, (m->numBd + (m->numBd+1)*2 + m->numBdEdges) * sizeof(int));
1568:     PetscMemcpy(tri->bdMarkers,   triOld->bdMarkers,    m->numBd      * sizeof(int));
1569:     PetscMemcpy(tri->bdEdgeBegin, triOld->bdEdgeBegin, (m->numBd+1)   * sizeof(int));
1570:     PetscMemcpy(tri->bdEdges,     triOld->bdEdges,      m->numBdEdges * sizeof(int));
1571:     tri->bdBegin[0]   = 0;
1572:     for(bd = 0; bd < m->numBd; bd++) {
1573:       tri->bdBegin[bd+1] = tri->bdBegin[bd];
1574:       for(bdNode = triOld->bdBegin[bd]; bdNode < triOld->bdBegin[bd+1]; bdNode++) {
1575:         node = triOld->bdNodes[bdNode];
1576:         AOMappingHasApplicationIndex(m->coarseMap, node, &hasIndex);
1577:         if (hasIndex == PETSC_TRUE)
1578:           tri->bdBegin[bd+1]++;
1579:       }
1580:     }
1581:     m->numBdNodes = tri->bdBegin[m->numBd];
1582:     /* Assume closed boundaries */
1583:     if (m->numBdNodes != m->numBdEdges) SETERRQ(PETSC_ERR_PLIB, "Invalid mesh boundary");
1584: #ifdef PETSC_USE_BOPT_g
1585:     MPI_Scan(&m->numBdNodes, &node, 1, MPI_INT, MPI_SUM, m->comm);
1586:     if (node != m->numBdNodes*(p->rank+1)) SETERRQ(PETSC_ERR_PLIB, "Invalid mesh boundary");
1587: #endif
1588:     PetscMalloc(m->numBdNodes * sizeof(int), &tri->bdNodes);
1589:     PetscLogObjectMemory(m, m->numBdNodes * sizeof(int));
1590:     for(bd = 0, count = 0; bd < m->numBd; bd++) {
1591:       for(bdNode = triOld->bdBegin[bd]; bdNode < triOld->bdBegin[bd+1]; bdNode++) {
1592:         node = triOld->bdNodes[bdNode];
1593:         AOMappingHasApplicationIndex(m->coarseMap, node, &hasIndex);
1594:         if (hasIndex == PETSC_TRUE) {
1595:           AOApplicationToPetsc(m->coarseMap, 1, &node);
1596:           tri->bdNodes[count++] = node;
1597:         }
1598:       }
1599:     }
1600:     if (count != m->numBdNodes) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node partition");

1602:     /* Partition boundary nodes */
1603:     PetscMalloc((p->numProcs+1) * sizeof(int), &q->firstBdNode);
1604:     PetscLogObjectMemory(p, (p->numProcs+1) * sizeof(int));
1605:     q->numLocBdNodes = 0;
1606:     for(bdNode = 0; bdNode < m->numBdNodes; bdNode++) {
1607:       newNode = tri->bdNodes[bdNode];
1608:       if ((newNode >= q->firstNode[p->rank]) && (newNode < q->firstNode[p->rank+1]))
1609:         q->numLocBdNodes++;
1610:     }
1611:     MPI_Allgather(&q->numLocBdNodes, 1, MPI_INT, &q->firstBdNode[1], 1, MPI_INT, p->comm);
1612:     for(proc = 1, q->firstBdNode[0] = 0; proc <= p->numProcs; proc++) {
1613:       q->firstBdNode[proc] += q->firstBdNode[proc-1];
1614:     }
1615:     q->numBdNodes    = q->firstBdNode[p->numProcs];

1617:     /* Process ghost boundary nodes */
1618:     q->numOverlapBdNodes = q->numLocBdNodes;
1619:     for(node = 0; node < numGhostNodes; node++) {
1620:       if (tri->markers[m->numNodes+node] != 0)
1621:         q->numOverlapBdNodes++;
1622:     }
1623:     q->ghostBdNodes      = PETSC_NULL;
1624:     if (q->numOverlapBdNodes > q->numLocBdNodes) {
1625:       PetscMalloc((q->numOverlapBdNodes - q->numLocBdNodes) * sizeof(int), &q->ghostBdNodes);
1626:       for(node = 0, bdNode = 0; node < numGhostNodes; node++) {
1627:         if (tri->markers[m->numNodes+node] != 0)
1628:           q->ghostBdNodes[bdNode++] = node;
1629:       }
1630:       if (bdNode != q->numOverlapBdNodes - q->numLocBdNodes) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary node partition");
1631:     }

1633:     /* Copy holes from previous mesh */
1634:     m->numHoles = mesh->numHoles;
1635:     m->holes    = PETSC_NULL;
1636:     if (m->numHoles > 0) {
1637:       PetscMalloc(m->numHoles*2 * sizeof(double), &m->holes);
1638:       PetscMemcpy(m->holes, mesh->holes, m->numHoles*2 * sizeof(double));
1639:       PetscLogObjectMemory(m, m->numHoles*2 * sizeof(double));
1640:     }

1642:     /* Calculate maximum degree of vertices */
1643:     m->maxDegree = mesh->maxDegree;
1644:     PetscMalloc(m->numNodes * sizeof(int), &tri->degrees);
1645:     PetscMalloc(m->maxDegree  * sizeof(int), &m->support);
1646:     for(newNode = 0; newNode < m->numNodes; newNode++) {
1647:       PartitionLocalToGlobalNodeIndex(p, newNode, &gNode);
1648:       AOPetscToApplication(m->coarseMap, 1, &gNode);
1649:       PartitionGlobalToLocalNodeIndex(pOld, gNode, &node);
1650:       tri->degrees[newNode] = triOld->degrees[node];
1651:     }
1652:     m->supportTaken = PETSC_FALSE;

1654: #ifdef PETSC_USE_BOPT_g
1655:     /* Check mesh integrity */
1656:     MeshDebug_Triangular_2D(m, PETSC_TRUE);
1657: #endif

1659:     /* Initialize save space */
1660:     tri->nodesOld = PETSC_NULL;

1662:   } else {
1663:     SETERRQ(PETSC_ERR_SUP, "No coarsening strategies currently supported");
1664:   }

1666:   PetscObjectSetName((PetscObject) m, "Coarse Mesh");
1667:   *newmesh = m;
1668:   return(0);
1669: }

1671: #undef  __FUNCT__
1673: /*
1674:   MeshRefine_Triangular_2D - Refines a two dimensional unstructured mesh using area constraints

1676:   Collective on Mesh

1678:   Input Parameters:
1679: + mesh    - The mesh begin refined
1680: - area    - A function which gives an area constraint when evaluated inside an element

1682:   Output Parameter:
1683: . newmesh - The refined mesh

1685:   Level: developer

1687: .keywords unstructured grid
1688: .seealso GridRefine(), MeshCoarsen(), MeshReform()
1689: */
1690: static int MeshRefine_Triangular_2D(Mesh mesh, PointFunction area, Mesh *newmesh)
1691: {
1692:   ParameterDict dict;
1693:   Mesh          m;
1694:   int         (*f)(Mesh, PointFunction, Mesh);
1695:   int           ierr;

1698:   MeshCreate(mesh->comm, &m);
1699:   ParameterDictCreate(mesh->comm, &dict);
1700:   ParameterDictSetObject(dict, "bdCtx", mesh->bdCtx);
1701:   ParameterDictSetInteger(dict, "dim", 2);
1702:   ParameterDictSetInteger(dict, "numLocNodes", mesh->numCorners);
1703:   PetscObjectSetParameterDict((PetscObject) m, dict);
1704:   ParameterDictDestroy(dict);

1706:   PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", (void (**)(void)) &f);
1707:   if (f == PETSC_NULL) SETERRQ(PETSC_ERR_SUP, "Mesh has no refinement function");
1708:   (*f)(mesh, area, m);
1709:   PetscObjectSetName((PetscObject) m, "Refined Mesh");
1710:   *newmesh = m;
1711:   return(0);
1712: }

1714: #undef  __FUNCT__
1716: static int MeshResetNodes_Triangular_2D(Mesh mesh, PetscTruth resetBd)
1717: {
1718:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1719:   int             *elements   = tri->faces;
1720:   int             *markers    = tri->markers;
1721:   double          *nodes      = tri->nodes;
1722:   int              numCorners = mesh->numCorners;
1723:   int              elem, corner, node1, node2, midNode, midCorner;

1726:   if (numCorners == 3)
1727:     return(0);

1729:   if (mesh->isPeriodic == PETSC_TRUE) {
1730:     for(elem = 0; elem < mesh->part->numOverlapElements; elem++) {
1731:       for(corner = 0; corner < 3; corner++) {
1732:         node1 = elements[elem*numCorners+corner];
1733:         node2 = elements[elem*numCorners+((corner+1)%3)];
1734:         if ((markers[node1] != markers[node2]) || (markers[node1] == 0) || (resetBd == PETSC_TRUE)) {
1735:           midCorner          = (corner+2)%3 + 3;
1736:           midNode            = elements[elem*numCorners+midCorner];
1737:           nodes[midNode*2]   = MeshPeriodicX(mesh, 0.5*(nodes[node1*2]   +
1738:                                                         MeshPeriodicRelativeX(mesh, nodes[node2*2],   nodes[node1*2])));
1739:           nodes[midNode*2+1] = MeshPeriodicY(mesh, 0.5*(nodes[node1*2+1] +
1740:                                                         MeshPeriodicRelativeY(mesh, nodes[node2*2+1], nodes[node1*2+1])));
1741:         }
1742:       }
1743:     }
1744:   } else {
1745:     for(elem = 0; elem < mesh->part->numOverlapElements; elem++) {
1746:       for(corner = 0; corner < 3; corner++) {
1747:         node1 = elements[elem*numCorners+corner];
1748:         node2 = elements[elem*numCorners+((corner+1)%3)];
1749:         if ((markers[node1] != markers[node2]) || (markers[node1] == 0) || (resetBd == PETSC_TRUE)) {
1750:           midCorner          = (corner+2)%3 + 3;
1751:           midNode            = elements[elem*numCorners+midCorner];
1752:           nodes[midNode*2]   = 0.5*(nodes[node1*2]   + nodes[node2*2]);
1753:           nodes[midNode*2+1] = 0.5*(nodes[node1*2+1] + nodes[node2*2+1]);
1754:         }
1755:       }
1756:     }
1757:   }
1758:   return(0);
1759: }

1761: #undef  __FUNCT__
1763: int MeshPartition_Triangular_2D(Mesh mesh)
1764: {

1768:   PartitionCreateTriangular2D(mesh, &mesh->part);
1769:   mesh->partitioned = 1;
1770:   PetscFunctionReturn(ierr);
1771: }

1773: #undef  __FUNCT__
1775: int MeshUpdateBoundingBox_Triangular_2D(Mesh mesh)
1776: {
1777:   Mesh_Triangular *tri   = (Mesh_Triangular *) mesh->data;
1778:   double          *nodes = tri->nodes;
1779:   Partition        part;
1780:   int              dim, numOverlapNodes;
1781:   int              node;
1782:   int              ierr;

1784:   /* Calculate the local bounding box */
1786:   MeshGetPartition(mesh, &part);
1787:   MeshGetDimension(mesh, &dim);
1788:   PartitionGetNumOverlapNodes(part, &numOverlapNodes);
1789:   if (numOverlapNodes > 0) {
1790:     mesh->locStartX = mesh->locEndX = nodes[0];
1791:     mesh->locStartY = mesh->locEndY = nodes[1];
1792:   } else {
1793:     mesh->locStartX = mesh->locEndX = 0.0;
1794:     mesh->locStartY = mesh->locEndY = 0.0;
1795:   }
1796:   for(node = 0; node < numOverlapNodes; node++) {
1797:     if (nodes[node*dim]   < mesh->locStartX) mesh->locStartX = nodes[node*dim];
1798:     if (nodes[node*dim]   > mesh->locEndX)   mesh->locEndX   = nodes[node*dim];
1799:     if (nodes[node*dim+1] < mesh->locStartY) mesh->locStartY = nodes[node*dim+1];
1800:     if (nodes[node*dim+1] > mesh->locEndY)   mesh->locEndY   = nodes[node*dim+1];
1801:   }
1802:   mesh->locSizeX = mesh->locEndX - mesh->locStartX;
1803:   mesh->locSizeY = mesh->locEndY - mesh->locStartY;
1804:   return(0);
1805: }

1807: #undef  __FUNCT__
1809: int MeshGetElemNeighbor_Triangular_2D(Mesh mesh, int elem, int corner, int *neighbor)
1810: {
1811:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1812: #ifdef PETSC_USE_BOPT_g
1813:   int              numOverlapElements;
1814:   int              ierr;
1815: #endif

1818: #ifdef PETSC_USE_BOPT_g
1819:   PartitionGetNumOverlapElements(mesh->part, &numOverlapElements);
1820:   if ((elem < 0)   || (elem >= numOverlapElements)) {
1821:     SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid element %d should be in [0,%d)", elem, numOverlapElements);
1822:   }
1823:   if ((corner < 0) || (corner >= mesh->numCorners)) {
1824:     SETERRQ2(PETSC_ERR_ARG_WRONG, "Invalid corner %d should be in [0,%d)", corner, mesh->numCorners);
1825:   }
1826: #endif
1827:   *neighbor = tri->neighbors[elem*3+corner];
1828:   return(0);
1829: }

1831: #undef  __FUNCT__
1833: int MeshGetNodeCoords_Triangular_2D(Mesh mesh, int node, double *x, double *y, double *z)
1834: {
1835:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1836:   int              dim = mesh->dim;
1837: #ifdef PETSC_USE_BOPT_g
1838:   int              numOverlapNodes;
1839:   int              ierr;
1840: #endif

1843: #ifdef PETSC_USE_BOPT_g
1844:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
1845:   if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified");
1846: #endif
1847:   if (x != PETSC_NULL) *x = tri->nodes[node*dim];
1848:   if (y != PETSC_NULL) *y = tri->nodes[node*dim+1];
1849:   if (z != PETSC_NULL) *z = 0.0;
1850:   return(0);
1851: }

1853: #undef  __FUNCT__
1855: int MeshSetNodeCoords_Triangular_2D(Mesh mesh, int node, double x, double y, double z)
1856: {
1857:   Mesh_Triangular *tri = (Mesh_Triangular *)         mesh->data;
1858:   int              dim = mesh->dim;
1859: #ifdef PETSC_USE_BOPT_g
1860:   int              numOverlapNodes;
1861:   int              ierr;
1862: #endif

1865: #ifdef PETSC_USE_BOPT_g
1866:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
1867:   if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified");
1868: #endif
1869:   tri->nodes[node*dim]   = MeshPeriodicX(mesh, x);
1870:   tri->nodes[node*dim+1] = MeshPeriodicY(mesh, y);
1871:   return(0);
1872: }

1874: #undef  __FUNCT__
1876: int MeshGetNodeCoordsSaved_Triangular_2D(Mesh mesh, int node, double *x, double *y, double *z)
1877: {
1878:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1879:   int              dim = mesh->dim;
1880: #ifdef PETSC_USE_BOPT_g
1881:   int              numOverlapNodes;
1882:   int              ierr;
1883: #endif

1886: #ifdef PETSC_USE_BOPT_g
1887:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
1888:   if ((node < 0) || (node >= numOverlapNodes)) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node specified");
1889: #endif
1890:   if (tri->nodesOld != PETSC_NULL) {
1891:     if (x != PETSC_NULL) *x = tri->nodesOld[node*dim];
1892:     if (y != PETSC_NULL) *y = tri->nodesOld[node*dim+1];
1893:     if (z != PETSC_NULL) *z = 0.0;
1894:   } else {
1895:     if (x != PETSC_NULL) *x = tri->nodes[node*dim];
1896:     if (y != PETSC_NULL) *y = tri->nodes[node*dim+1];
1897:     if (z != PETSC_NULL) *z = 0.0;
1898:   }
1899:   return(0);
1900: }

1902: #undef  __FUNCT__
1904: int MeshGetNodeBoundary_Triangular_2D(Mesh mesh, int node, int *bd)
1905: {
1906:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1907: #ifdef PETSC_USE_BOPT_g
1908:   int              numOverlapNodes;
1909:   int              ierr;
1910: #endif

1913: #ifdef PETSC_USE_BOPT_g
1914:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
1915:   if ((node < 0) || (node >= numOverlapNodes)) {
1916:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", node, numOverlapNodes);
1917:   }
1918: #endif
1919:   *bd = tri->markers[node];
1920:   return(0);
1921: }

1923: #undef  __FUNCT__
1925: int MeshGetNodeFromElement_Triangular_2D(Mesh mesh, int elem, int corner, int *node)
1926: {
1927:   Mesh_Triangular *tri        = (Mesh_Triangular *) mesh->data;
1928:   int              numCorners = mesh->numCorners;

1931:   *node = tri->faces[elem*numCorners+corner];
1932:   return(0);
1933: }

1935: #undef  __FUNCT__
1937: int MeshGetNodeFromEdge_Triangular_2D(Mesh mesh, int edge, int corner, int *node)
1938: {
1939:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;

1942:   *node = tri->edges[edge*2+corner];
1943:   return(0);
1944: }

1946: #undef  __FUNCT__
1948: int MeshNodeIsVertex_Triangular_2D(Mesh mesh, int node, PetscTruth *isVertex)
1949: {
1950:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1951: #ifdef PETSC_USE_BOPT_g
1952:   int              numOverlapNodes;
1953:   int              ierr;
1954: #endif

1957: #ifdef PETSC_USE_BOPT_g
1958:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
1959:   if ((node < 0) || (node >= numOverlapNodes)) {
1960:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Invalid node %d should be in [0,%d)", node, numOverlapNodes);
1961:   }
1962: #endif
1963:   if (tri->degrees[node] == 0)
1964:     *isVertex = PETSC_FALSE;
1965:   else
1966:     *isVertex = PETSC_TRUE;
1967:   return(0);
1968: }

1970: #undef  __FUNCT__
1972: int MeshGetBoundarySize_Triangular_2D(Mesh mesh, int boundary, int *size)
1973: {
1974:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1975:   int              b; /* Canonical boundary number */
1976:   int              ierr;

1979:   /* Find canonical boundary number */
1980:   MeshGetBoundaryIndex(mesh, boundary, &b);
1981:   *size = tri->bdBegin[b+1] - tri->bdBegin[b];
1982:   return(0);
1983: }

1985: #undef  __FUNCT__
1987: int MeshGetBoundaryIndex_Triangular_2D(Mesh mesh, int boundary, int *index)
1988: {
1989:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1990:   int              b; /* Canonical boundary number */

1993:   /* Find canonical boundary number */
1994:   for(b = 0; b < mesh->numBd; b++)
1995:     if (tri->bdMarkers[b] == boundary) break;
1996:   if (b == mesh->numBd) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid boundary marker %d", boundary);
1997:   *index = b;
1998:   return(0);
1999: }

2001: #undef  __FUNCT__
2003: int MeshGetBoundaryStart_Triangular_2D(Mesh mesh, int boundary, PetscTruth ghost, int *node)
2004: {
2005:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
2006:   int              b; /* Canonical boundary number */
2007:   int              ierr;

2010:   /* Find canonical boundary number */
2011:   MeshGetBoundaryIndex(mesh, boundary, &b);
2012:   if (mesh->activeBd != -1) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE, "Already iterating over boundary %d", mesh->activeBd);
2013:   /* Find first boundary node */
2014:   mesh->activeBd     = b;
2015:   mesh->activeBdOld  = b;
2016:   mesh->activeBdNode = tri->bdBegin[b] - 1;
2017:   MeshGetBoundaryNext(mesh, boundary, ghost, node);
2018:   return(0);
2019: }

2021: #undef  __FUNCT__
2023: int MeshGetBoundaryNext_Triangular_2D(Mesh mesh, int boundary, PetscTruth ghost, int *node)
2024: {
2025:   Mesh_Triangular         *tri      = (Mesh_Triangular *) mesh->data;
2026:   Partition                p        = mesh->part;
2027:   Partition_Triangular_2D *q        = (Partition_Triangular_2D *) p->data;
2028:   int                      numNodes = mesh->numNodes;
2029:   int                      b;        /* Canonical boundary number */
2030:   int                      tempNode; /* Canonical local node number */
2031:   int                      ierr;

2034:   /* Find canonical boundary number */
2035:   MeshGetBoundaryIndex(mesh, boundary, &b);
2036:   /* Find next boundary node */
2037:   if (mesh->activeBd != b) SETERRQ1(PETSC_ERR_ARG_WRONG, "Boundary %d not active", boundary);
2038:   do {
2039:     mesh->activeBdNode++;
2040:     if (mesh->activeBdNode == tri->bdBegin[b+1]) {
2041:       mesh->activeBd = mesh->activeBdNode = -1;
2042:       *node = -1;
2043:       return(0);
2044:     }
2045:     tempNode = tri->bdNodes[mesh->activeBdNode] - q->firstNode[p->rank];
2046:     /* Translate ghost nodes */
2047:     if (ghost == PETSC_TRUE) {
2048:       if ((tempNode < 0) || (tempNode >= numNodes)) {
2049:         PartitionGhostNodeIndex_Private(p, tri->bdNodes[mesh->activeBdNode], &tempNode);
2050:         if (ierr == 0) {
2051:           tempNode += numNodes;
2052:           break;
2053:         }
2054:       }
2055:     }
2056:   }
2057:   while((tempNode < 0) || (tempNode >= numNodes));
2058:   *node = tempNode;
2059:   return(0);
2060: }

2062: #undef  __FUNCT__
2064: int MeshGetActiveBoundary_Triangular_2D(Mesh mesh, int *boundary)
2065: {
2066:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;

2069:   if (mesh->activeBdOld >= 0) {
2070:     *boundary = tri->bdMarkers[mesh->activeBdOld];
2071:   } else {
2072:     *boundary = 0;
2073:   }
2074:   return(0);
2075: }

2077: #undef  __FUNCT__
2079: int MeshGetNodeSupport_Triangular_2D(Mesh mesh, int node, int startElem, int *degree, int **support)
2080: {
2081:   Mesh_Triangular *tri                = (Mesh_Triangular *) mesh->data;
2082:   int              numOverlapElements = mesh->numFaces;
2083:   int              numCorners         = mesh->numCorners;
2084:   int             *elements           = tri->faces;
2085:   int             *neighbors          = tri->neighbors;
2086:   int              deg                = -1;
2087:   PetscTruth       found;
2088:   int              firstElem;
2089:   int              prevElem, nextElem;
2090:   int              elem, neighbor, corner;

2093:   /* First check suggested element */
2094:   firstElem = -1;
2095:   for(corner = 0; corner < numCorners; corner++) {
2096:     if (elements[startElem*numCorners+corner] == node) {
2097:       firstElem = startElem;
2098:       if (corner >= 3) {
2099:         deg = 1;
2100:       } else {
2101:         deg = 0;
2102:       }
2103:       break;
2104:     }
2105:   }

2107:   /* Locate an element containing the node */
2108:   if (mesh->part != PETSC_NULL) {
2109:     numOverlapElements = mesh->part->numOverlapElements;
2110:   }
2111:   for(elem = 0; (elem < numOverlapElements) && (firstElem < 0); elem++) {
2112:     for(corner = 0; corner < numCorners; corner++) {
2113:       if (elements[elem*numCorners+corner] == node) {
2114:         firstElem = elem;
2115:         if (corner >= 3) {
2116:           deg = 1;
2117:         } else {
2118:           deg = 0;
2119:         }
2120:         break;
2121:       }
2122:     }
2123:   }
2124:   if (firstElem < 0) {
2125:     SETERRQ1(PETSC_ERR_ARG_WRONG, "Node %d not found in mesh", node);
2126:   }

2128:   /* Check for midnode */
2129:   if (deg != 0) {
2130:     mesh->support[0] = firstElem;
2131:     mesh->support[1] = -1;
2132:     /* Find other element containing node */
2133:     for(neighbor = 0; neighbor < 3; neighbor++) {
2134:       elem = neighbors[firstElem*3+neighbor];
2135:       if (elem < 0)
2136:         continue;

2138:       for(corner = 3; corner < numCorners; corner++) {
2139:         if (elements[elem*numCorners+corner] == node) {
2140:           mesh->support[1] = elem;
2141:           deg++;
2142:           break;
2143:         }
2144:       }
2145:       if (corner < numCorners) break;
2146:     }
2147: #ifdef PETSC_USE_BOPT_g
2148:     if (mesh->support[1] == -1) {
2149:       if (tri->markers[node] == 0) SETERRQ1(PETSC_ERR_PLIB, "Interior node %d with incomplete support", node);
2150:     }
2151: #endif
2152:     *support = mesh->support;
2153:     *degree  = deg;
2154:     return(0);
2155:   }

2157:   /* Walk around node */
2158:   prevElem = -1;
2159:   nextElem = -1;
2160:   found    = PETSC_TRUE;
2161:   while((nextElem != firstElem) && (found == PETSC_TRUE)) {
2162:     /* First time through */
2163:     if (prevElem == -1)
2164:       nextElem = firstElem;

2166:     /* Add element to the list */
2167:     mesh->support[deg++] = nextElem;

2169:     /* Look for a neighbor containing the node */
2170:     found  = PETSC_FALSE;
2171:     for(neighbor = 0; neighbor < 3; neighbor++) {
2172:       /* Reject the element we just came from (we could reject the opposite face also) */
2173:       elem = neighbors[nextElem*3+neighbor];
2174:       if ((elem == prevElem) || (elem < 0)) continue;

2176:       for(corner = 0; corner < 3; corner++) {
2177:         if (elements[elem*numCorners+corner] == node) {
2178:           prevElem = nextElem;
2179:           nextElem = elem;
2180:           found    = PETSC_TRUE;
2181:           break;
2182:         }
2183:       }
2184:       if (corner < 3) break;
2185:     }

2187:     if (found == PETSC_FALSE) {
2188: #ifdef PETSC_USE_BOPT_g
2189:       /* Make sure that this is a boundary node */
2190:       if (tri->markers[node] == 0) SETERRQ(PETSC_ERR_PLIB, "Interior node with incomplete support");
2191: #endif

2193:       /* Also walk the other way for a boundary node --
2194:            We could be nice and reorder the elements afterwards to all be adjacent */
2195:       if (deg > 1) {
2196:         prevElem = mesh->support[1];
2197:         nextElem = mesh->support[0];
2198:         found    = PETSC_TRUE;
2199:         while(found == PETSC_TRUE) {
2200:           /* Look for a neighbor containing the node */
2201:           found    = PETSC_FALSE;
2202:           for(neighbor = 0; neighbor < 3; neighbor++) {
2203:             /* Reject the element we just came from (we could reject the opposite face also) */
2204:             elem = neighbors[nextElem*3+neighbor];
2205:             if ((elem == prevElem) || (elem < 0)) continue;

2207:             for(corner = 0; corner < 3; corner++) {
2208:               if (elements[elem*numCorners+corner] == node) {
2209:                 prevElem = nextElem;
2210:                 nextElem = elem;
2211:                 found    = PETSC_TRUE;

2213:                 /* Add element to the list */
2214:                 mesh->support[deg++] = nextElem;
2215:                 break;
2216:               }
2217:             }
2218:             if (corner < 3) break;
2219:           }
2220:         }
2221:       }
2222:     }

2224: #ifdef PETSC_USE_BOPT_g
2225:     /* Check for overflow */
2226:     if (deg > mesh->maxDegree) SETERRQ1(PETSC_ERR_MEM, "Maximum degree %d exceeded", mesh->maxDegree);
2227: #endif
2228:   }

2230:   *support = mesh->support;
2231:   *degree  = deg;
2232:   return(0);
2233: }

2235: #if 0
2236: #undef  __FUNCT__
2238: int MeshLocatePoint_Triangular_2D_Directed(Mesh mesh, double xx, double yy, double zz, int *containingElem)
2239: {
2240:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
2241:   double           start_x, start_y; /* Coordinates of a vertex */
2242:   double           v_x, v_y;         /* Coordinates of an edge vector */
2243:   double           p_x, p_y;         /* Coordinates of the vector from a vertex (start) to a point */
2244:   double           cross;            /* The magnitude of the cross product of two vectors */
2245:   PetscTruth       inside[3];        /* Are we in the correct halfplane? */
2246:   PetscTruth       found;            /* We have found the element containing a point */
2247:   int              nextElem;         /* The next element to search */
2248:   int              prevElem;         /* The last element searched */
2249:   int              numCorners  = mesh->numCorners;
2250:   int              numElements = mesh->numFaces;
2251:   double          *nodes       = tri->nodes;
2252:   int             *elements    = tri->faces;
2253:   int             *neighbors   = tri->neighbors;
2254:   double           x           = xx;
2255:   double           y           = yy;
2256:   double           coords[6];
2257:   int              neighbor;
2258:   int              elem, edge;

2261:   *containingElem = -1;
2262:   /* Quick bounding rectangle check */
2263:   if ((mesh->isPeriodic == PETSC_FALSE) &&
2264:       ((mesh->locStartX > x) || (mesh->locEndX < x) || (mesh->locStartY > y) || (mesh->locEndY < y)))
2265:     return(0);
2266:   /* Simple search */
2267:   for(elem = 0, found = PETSC_FALSE; elem < numElements; elem++) {
2268:     /* A point lies within a triangle is it is on the inner side of every edge --
2269:          We take the cross product of the edge oriented counterclockwise with
2270:          the vector from the edge start to the point:

2272:          2
2273:          |
2274:          | 
2275:          |  
2276:          0---1----> v_0
2277:           
2278:             theta
2279:             
2280:              P

2282:          If

2284:            |v_i cross P| equiv |v_i| |P| sintheta > 0

2286:          then the point lies in the correct half-plane. This has the advantage of
2287:          only using elementary arithmetic. Robust predicates could be used to catch
2288:          points near a boundary.
2289:     */
2290:     if (mesh->isPeriodic == PETSC_TRUE) {
2291:       coords[0] = nodes[elements[elem*numCorners]*2];
2292:       coords[1] = nodes[elements[elem*numCorners]*2+1];
2293:       coords[2] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+1]*2],   nodes[elements[elem*numCorners]*2]);
2294:       coords[3] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+1]*2+1], nodes[elements[elem*numCorners]*2+1]);
2295:       coords[4] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+2]*2],   nodes[elements[elem*numCorners]*2]);
2296:       coords[5] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+2]*2+1], nodes[elements[elem*numCorners]*2+1]);
2297:       x         = MeshPeriodicRelativeX(mesh, xx, nodes[elements[elem*numCorners]*2]);
2298:       y         = MeshPeriodicRelativeY(mesh, yy, nodes[elements[elem*numCorners]*2+1]);
2299:     }
2300:     for(edge = 0; edge < 3; edge++) {
2301:       /* First edge 0--1 */
2302:       if (mesh->isPeriodic == PETSC_TRUE) {
2303:         start_x = coords[edge*2];
2304:         start_y = coords[edge*2+1];
2305:         v_x     = coords[((edge+1)%3)*2]   - start_x;
2306:         v_y     = coords[((edge+1)%3)*2+1] - start_y;
2307:       } else {
2308:         start_x = nodes[elements[elem*numCorners+edge]*2];
2309:         start_y = nodes[elements[elem*numCorners+edge]*2+1];
2310:         v_x     = nodes[elements[elem*numCorners+((edge+1)%3)]*2]   - start_x;
2311:         v_y     = nodes[elements[elem*numCorners+((edge+1)%3)]*2+1] - start_y;
2312:       }
2313:       p_x       = x - start_x;
2314:       p_y       = y - start_y;
2315:       cross     = v_x*p_y - v_y*p_x;

2317:       /* We will use this routine to find integration points which lie on the
2318:          boundary so we must give a little slack to the test. I do not think
2319:          this will bias our results. */
2320:       if ((cross >= 0) || (PetscAbsScalar(cross) < PETSC_MACHINE_EPSILON)) {
2321:         inside[edge] = PETSC_TRUE;
2322:       } else {
2323:         inside[edge] = PETSC_FALSE;
2324:       }
2325:     }

2327:     /* This information could be used for a directed search --

2329:           4 T|
2330:             F|
2331:             F|
2332:              |
2333:              |
2334:               2
2335:           T   |        T
2336:         5 T   |        F 3
2337:           F   |T       T
2338:           ----0---1------
2339:           F   | F      F
2340:         6 T   | T 1    F 2
2341:           F   | T      T

2343:        and clearly all F's means that the area of the triangle is 0.
2344:        The regions are ordered so that all the Hamming distances are
2345:        1 (Gray code).
2346:     */
2347:     if (inside[2] == PETSC_TRUE) {
2348:       if (inside[1] == PETSC_TRUE) {
2349:         if (inside[0] == PETSC_TRUE) {
2350:           /* Region 0 -- TTT */
2351:           found = PETSC_TRUE;
2352:         } else {
2353:           /* Region 1 -- FTT */
2354:           nextElem = neighbors[elem*3+2];
2355:         }
2356:       } else {
2357:         if (inside[0] == PETSC_TRUE) {
2358:           /* Region 3 -- TFT */
2359:           nextElem = neighbors[elem*3];
2360:         } else {
2361:           /* Region 2 -- FFT */
2362:           if (neighbors[elem*3] >= 0) {
2363:             nextElem = neighbors[elem*3];
2364:           } else {
2365:             nextElem = neighbors[elem*3+2];
2366:           }
2367:         }
2368:       }
2369:     } else {
2370:       if (inside[1] == PETSC_TRUE) {
2371:         if (inside[0] == PETSC_TRUE) {
2372:           /* Region 5 -- TTF */
2373:           nextElem = neighbors[elem*3+1];
2374:         } else {
2375:           /* Region 6 -- FTF */
2376:           if (neighbors[elem*3+1] >= 0) {
2377:             nextElem = neighbors[elem*3+1];
2378:           } else {
2379:             nextElem = neighbors[elem*3+2];
2380:           }
2381:         }
2382:       } else {
2383:         if (inside[0] == PETSC_TRUE) {
2384:           /* Region 4 -- TFF */
2385:           if (neighbors[elem*3] >= 0) {
2386:             nextElem = neighbors[elem*3];
2387:           } else {
2388:             nextElem = neighbors[elem*3+1];
2389:           }
2390:         } else {
2391:           /* Region 7 -- FFF */
2392:           PetscPrintf(PETSC_COMM_SELF, "Element %d: (%g, %g)-(%g, %g)-(%g, %g) and (%g, %g)n", elem,
2393:                       nodes[elements[elem*numCorners]*2],   nodes[elements[elem*numCorners]*2+1],
2394:                       nodes[elements[elem*numCorners+1]*2], nodes[elements[elem*numCorners+1]*2+1],
2395:                       nodes[elements[elem*numCorners+2]*2], nodes[elements[elem*numCorners+2]*2+1], x, y);
2396:           SETERRQ1(PETSC_ERR_MESH_NULL_ELEM, "Element %d has no area", elem);
2397:         }
2398:       }
2399:     }

2401:     /* We found the point */
2402:     if (found) {
2403:       *containingElem = elem;
2404:       break;
2405:     }

2407:     /* Choose next direction */
2408:     if (nextElem < 0) {
2409:       /* Choose randomly */
2410:       for(neighbor = 0; neighbor < 3; neighbor++) {
2411:         tempElem = neighbors[elem*3+neighbor];
2412:         if ((tempElem >= 0) && (tempElem != prevElem)) {
2413:           nextElem = tempElem;
2414:           break;
2415:         }
2416:       }
2417:       if (nextElem < 0) {
2418:         SETERRQ1(PETSC_ERR_MAX_ITER, "Element search hit dead end in element %d", elem);
2419:       }
2420:     }

2422:     prevElem = elem;
2423:   }

2425:   return(0);
2426: }
2427: #endif

2429: #undef  __FUNCT__
2431: int MeshLocatePoint_Triangular_2D(Mesh mesh, double xx, double yy, double zz, int *containingElem) {
2432:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
2433:   double           start_x, start_y; /* Coordinates of a vertex */
2434:   double           v_x, v_y;         /* Coordinates of an edge vector */
2435:   double           p_x, p_y;         /* Coordinates of the vector from a vertex (start) to a point */
2436:   double           cross;            /* The magnitude of the cross product of two vectors */
2437:   PetscTruth       inside[3];        /* Are we in the correct halfplane? */
2438:   PetscTruth       found;            /* We have found the element containing a point */
2439:   int              numCorners  = mesh->numCorners;
2440:   int              numElements = mesh->numFaces;
2441:   double          *nodes       = tri->nodes;
2442:   int             *markers     = tri->markers;
2443:   int             *elements    = tri->faces;
2444:   int             *neighbors   = tri->neighbors;
2445:   double           x           = xx;
2446:   double           y           = yy;
2447:   double           coords[6];
2448:   double           pNormSq;
2449:   int              elem, edge;

2452:   *containingElem = -1;
2453:   /* Quick bounding rectangle check */
2454:   if ((mesh->isPeriodic == PETSC_FALSE) &&
2455:       ((mesh->locStartX > x) || (mesh->locEndX < x) || (mesh->locStartY > y) || (mesh->locEndY < y)))
2456:     return(0);
2457:   /* Simple search */
2458:   for(elem = 0, found = PETSC_FALSE; elem < numElements; elem++) {
2459:     /* A point lies within a triangle is it is on the inner side of every edge --
2460:          We take the cross product of the edge oriented counterclockwise with
2461:          the vector from the edge start to the point:

2463:          2
2464:          |
2465:          | 
2466:          |  
2467:          0---1----> v_0
2468:           
2469:             theta
2470:             
2471:              P

2473:          If

2475:            |v_i cross P| equiv |v_i| |P| sintheta > 0

2477:          then the point lies in the correct half-plane. This has the advantage of
2478:          only using elementary arithmetic. Robust predicates could be used to catch
2479:          points near a boundary.
2480:     */
2481:     if (mesh->isPeriodic == PETSC_TRUE) {
2482:       coords[0] = nodes[elements[elem*numCorners]*2];
2483:       coords[1] = nodes[elements[elem*numCorners]*2+1];
2484:       coords[2] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+1]*2],   nodes[elements[elem*numCorners]*2]);
2485:       coords[3] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+1]*2+1], nodes[elements[elem*numCorners]*2+1]);
2486:       coords[4] = MeshPeriodicRelativeX(mesh, nodes[elements[elem*numCorners+2]*2],   nodes[elements[elem*numCorners]*2]);
2487:       coords[5] = MeshPeriodicRelativeY(mesh, nodes[elements[elem*numCorners+2]*2+1], nodes[elements[elem*numCorners]*2+1]);
2488:       x         = MeshPeriodicRelativeX(mesh, xx, nodes[elements[elem*numCorners]*2]);
2489:       y         = MeshPeriodicRelativeY(mesh, yy, nodes[elements[elem*numCorners]*2+1]);
2490:     }
2491:     for(edge = 0; edge < 3; edge++) {
2492:       /* First edge 0--1 */
2493:       if (mesh->isPeriodic == PETSC_TRUE) {
2494:         start_x = coords[edge*2];
2495:         start_y = coords[edge*2+1];
2496:         v_x     = coords[((edge+1)%3)*2]   - start_x;
2497:         v_y     = coords[((edge+1)%3)*2+1] - start_y;
2498:       } else {
2499:         start_x = nodes[elements[elem*numCorners+edge]*2];
2500:         start_y = nodes[elements[elem*numCorners+edge]*2+1];
2501:         v_x     = nodes[elements[elem*numCorners+((edge+1)%3)]*2]   - start_x;
2502:         v_y     = nodes[elements[elem*numCorners+((edge+1)%3)]*2+1] - start_y;
2503:       }
2504:       p_x       = x - start_x;
2505:       p_y       = y - start_y;
2506:       cross     = v_x*p_y - v_y*p_x;
2507:       pNormSq   = p_x*p_x+p_y*p_y;

2509:       /* Check for identical points, being looser at the boundary */
2510:       if (pNormSq < PETSC_MACHINE_EPSILON) {
2511:         *containingElem = elem;
2512:         return(0);
2513:       } else if (markers[elements[elem*numCorners+edge]] != 0) {
2514:         if (pNormSq < 0.001*PetscMax(v_x*v_x, v_y*v_y)) {
2515:           *containingElem = elem;
2516:           return(0);
2517:         }
2518:       }

2520:       if (cross >= 0) {
2521:           inside[edge] = PETSC_TRUE;
2522:       } else if ((neighbors[elem*3+((edge+2)%3)] < 0) && (cross*cross < 0.087155743*(v_x*v_x+v_y*v_y)*pNormSq)) {
2523:         /* We are losing precision here (I think) for large problems (since points are turning up inside particles)
2524:           and therefore I am putting a 5 degree limit on this check (sin 5 = 0.087155743), which is only active for
2525:           boundary edges. */
2526:         inside[edge] = PETSC_TRUE;
2527:       } else {
2528:         inside[edge] = PETSC_FALSE;
2529:       }
2530:     }

2532:     /* This information could be used for a directed search --

2534:           4 T|
2535:             F|
2536:             F|
2537:              |
2538:              |
2539:               2
2540:           T   |        T
2541:         5 T   |        F 3
2542:           F   |T       T
2543:           ----0---1------
2544:           F   | F      F
2545:         6 T   | T 1    F 2
2546:           F   | T      T

2548:        and clearly all F's means that the area of the triangle is 0.
2549:        The regions are ordered so that all the Hamming distances are
2550:        1 (Gray code).
2551:     */
2552:     if ((inside[2] == PETSC_TRUE) && (inside[1] == PETSC_TRUE) && (inside[0] == PETSC_TRUE)) {
2553:       double minX = coords[0];
2554:       double maxX = coords[0];
2555:       double minY = coords[1];
2556:       double maxY = coords[1];
2557:       double offX, offY;

2559:       if (coords[2] < minX) minX = coords[2];
2560:       if (coords[2] > maxX) maxX = coords[2];
2561:       if (coords[3] < minY) minY = coords[3];
2562:       if (coords[3] > maxY) maxY = coords[3];
2563:       if (coords[4] < minX) minX = coords[4];
2564:       if (coords[4] > maxX) maxX = coords[4];
2565:       if (coords[5] < minY) minY = coords[5];
2566:       if (coords[5] > maxY) maxY = coords[5];
2567:       offX = 0.01*(maxX - minX);
2568:       offY = 0.01*(maxY - minY);

2570:       if ((x >= minX-offX) && (x <= maxX+offX) && (y >= minY-offY) && (y <= maxY+offY)) {
2571:         /* Region 0 -- TTT */
2572:         found = PETSC_TRUE;
2573:       }
2574:     } else if ((inside[2] == PETSC_FALSE) && (inside[1] == PETSC_FALSE) && (inside[0] == PETSC_FALSE)) {
2575:       /* Region 7 -- FFF */
2576:       PetscPrintf(PETSC_COMM_SELF, "Element %d: (%g, %g)-(%g, %g)-(%g, %g) and (%g, %g)n", elem, coords[0], coords[1], coords[2],
2577:                   coords[3], coords[4], coords[5], x, y);
2578:       SETERRQ1(PETSC_ERR_MESH_NULL_ELEM, "Element %d has no area or is inverted", elem);
2579:     }

2581:     /* We found the point */
2582:     if (found) {
2583:       *containingElem = elem;
2584:       break;
2585:     }
2586:   }

2588:   return(0);
2589: }

2591: #undef  __FUNCT__
2593: int MeshNearestNode_Triangular_2D(Mesh mesh, double x, double y, double z, PetscTruth outsideDomain, int *node)
2594: {
2595:   Mesh_Triangular         *tri        = (Mesh_Triangular *) mesh->data;
2596:   Partition                p          = mesh->part;
2597:   Partition_Triangular_2D *q          = (Partition_Triangular_2D *) p->data;
2598:   int                      numCorners = mesh->numCorners;
2599:   int                     *elements   = tri->faces;
2600:   int                      numNodes   = mesh->numNodes;
2601:   double                  *nodes      = tri->nodes;
2602:   int                     *firstNode  = q->firstNode;
2603:   int                      rank       = p->rank;
2604:   int                      elem;
2605:   double                   minDist, dist, allMinDist;
2606:   int                      minNode, globalMinNode, allMinNode;
2607:   int                      corner, n;
2608:   int                      ierr;

2611:   if (outsideDomain == PETSC_FALSE) {
2612:     MeshLocatePoint(mesh, x, y, z, &elem);
2613:     if (elem >= 0) {
2614:       minDist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners]*2] - x)) +
2615:         PetscSqr(MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners]*2+1] - y));
2616:       minNode = elements[elem*numCorners];
2617:       for(corner = 1; corner < numCorners; corner++) {
2618:         dist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[elements[elem*numCorners+corner]*2] - x)) +
2619:           PetscSqr(MeshPeriodicDiffY(mesh, nodes[elements[elem*numCorners+corner]*2+1] - y));
2620:         if (dist < minDist) {
2621:           minDist = dist;
2622:           minNode = elements[elem*numCorners+corner];
2623:         }
2624:       }
2625:       PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
2626:     } else {
2627:       minNode       = -1;
2628:       globalMinNode = -1;
2629:     }

2631:     /* The minimum node might still be a ghost node */
2632:     MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, p->comm);
2633:     if ((allMinNode >= firstNode[rank]) && (allMinNode < firstNode[rank+1])) {
2634:       *node = allMinNode - firstNode[rank];
2635:     } else {
2636:       *node = -1;
2637:     }
2638:     if (allMinNode < 0)
2639:     PetscFunctionReturn(1);
2640:   } else {
2641:     /* Brute force check */
2642:     minNode = 0;
2643:     minDist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[0*2] - x)) +
2644:       PetscSqr(MeshPeriodicDiffY(mesh, nodes[0*2+1] - y));
2645:     for(n = 1; n < numNodes; n++) {
2646:       dist = PetscSqr(MeshPeriodicDiffX(mesh, nodes[n*2] - x)) +
2647:         PetscSqr(MeshPeriodicDiffY(mesh, nodes[n*2+1] - y));
2648:       if (dist < minDist) {
2649:         minDist = dist;
2650:         minNode = n;
2651:       }
2652:     }

2654:     /* Find the minimum distance */
2655:     MPI_Allreduce(&minDist, &allMinDist, 1, MPI_DOUBLE, MPI_MIN, p->comm);
2656:     if (minDist == allMinDist) {
2657:       PartitionLocalToGlobalNodeIndex(p, minNode, &globalMinNode);
2658:     } else {
2659:       globalMinNode = -1;
2660:     }

2662:     /* We might still have ties */
2663:     MPI_Allreduce(&globalMinNode, &allMinNode, 1, MPI_INT, MPI_MAX, p->comm);
2664:     if ((allMinNode >= firstNode[rank]) && (allMinNode < firstNode[rank+1])) {
2665:       *node = allMinNode - firstNode[rank];
2666:     } else {
2667:       *node = -1;
2668:     }
2669:     if (allMinNode < 0)
2670:       PetscFunctionReturn(1);
2671:   }

2673:   return(0);
2674: }

2676: #undef  __FUNCT__
2678: int MeshSaveMesh_Triangular_2D(Mesh mesh) {
2679:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
2680:   int              dim = mesh->dim;
2681:   int              numOverlapNodes;
2682:   int              ierr;

2685:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
2686:   if (tri->nodesOld == PETSC_NULL) {
2687:     PetscMalloc(numOverlapNodes*dim * sizeof(double), &tri->nodesOld);
2688:     PetscLogObjectMemory(mesh, numOverlapNodes*dim * sizeof(double));
2689:   }
2690:   PetscMemcpy(tri->nodesOld, tri->nodes, numOverlapNodes*dim * sizeof(double));
2691:   return(0);
2692: }

2694: #undef  __FUNCT__
2696: int MeshRestoreMesh_Triangular_2D(Mesh mesh) {
2697:   Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
2698:   int              dim = mesh->dim;
2699:   int              numOverlapNodes;
2700:   int              ierr;

2703:   ierr  = PartitionGetNumOverlapNodes(mesh->part, &numOverlapNodes);
2704:   if (tri->nodesOld != PETSC_NULL) {
2705:     PetscMemcpy(tri->nodes, tri->nodesOld, numOverlapNodes*dim * sizeof(double));
2706:   }
2707:   return(0);
2708: }

2710: #undef  __FUNCT__
2712: int MeshIsDistorted_Triangular_2D(Mesh mesh, PetscTruth *flag) {
2713:   Mesh_Triangular *tri            = (Mesh_Triangular *) mesh->data;
2714:   double          *nodes          = tri->nodes;
2715:   int             *faces          = tri->faces;
2716:   int              numCorners     = mesh->numCorners;
2717:   int              numElements    = mesh->numFaces;
2718:   double           maxAspectRatio = mesh->maxAspectRatio;
2719:   double           area;        /* Twice the area of the triangle */
2720:   double           aspectRatio; /* (Longest edge)^2/(Twice the area of the triangle) */
2721:   double           max;
2722:   PetscTruth       locFlag;
2723:   double           x21, y21, x31, y31, x32, y32;
2724:   int              elem;
2725:   int              locRetVal;
2726:   int              retVal;
2727:   int              ierr;

2730:   for(elem = 0, locFlag = PETSC_FALSE, locRetVal = 0, max = 0.0; elem < numElements; elem++) {
2731:     /* Find the area */
2732:     if (mesh->isPeriodic == PETSC_TRUE) {
2733:       x21 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2]);
2734:       y21 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1]);
2735:       x31 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2]);
2736:       y31 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1]);
2737:       x32 = MeshPeriodicDiffX(mesh, nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners+1]*2]);
2738:       y32 = MeshPeriodicDiffY(mesh, nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1]);
2739:     } else {
2740:       x21 = nodes[faces[elem*numCorners+1]*2]   - nodes[faces[elem*numCorners]*2];
2741:       y21 = nodes[faces[elem*numCorners+1]*2+1] - nodes[faces[elem*numCorners]*2+1];
2742:       x31 = nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners]*2];
2743:       y31 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners]*2+1];
2744:       x32 = nodes[faces[elem*numCorners+2]*2]   - nodes[faces[elem*numCorners+1]*2];
2745:       y32 = nodes[faces[elem*numCorners+2]*2+1] - nodes[faces[elem*numCorners+1]*2+1];
2746:     }
2747:     area = x21*y31 - x31*y21;
2748:     /* Check for inverted elements */
2749:     if (area < 0.0) {
2750:       locFlag   = PETSC_TRUE;
2751:       locRetVal = 1;
2752:       break;
2753:     } else if (area == 0.0) {
2754:       SETERRQ1(PETSC_ERR_ARG_CORRUPT, "Element %d has no area", elem);
2755:     }

2757:     /* Find the aspect ratio = (longest edge length)^2 / (2 area) */
2758:     aspectRatio  = PetscMax(x21*x21 + y21*y21, x31*x31 + y31*y31);
2759:     aspectRatio  = PetscMax(x32*x32 + y32*y32, aspectRatio);
2760:     aspectRatio /= area;
2761:     /* We cannot bail here since we must check for inverted elements */
2762:     max = PetscMax(max, aspectRatio);
2763:   }
2764:   if (max > maxAspectRatio) {
2765:     PetscLogInfo(mesh, "Aspect ratio too large in element %d: %g > %gn", elem, max, maxAspectRatio);
2766:     locFlag = PETSC_TRUE;
2767:   } else {
2768:     PetscLogInfo(mesh, "Maximum aspect ratio in element %d: %gn", elem, max);
2769:   }
2770:   PetscLogFlops(elem*19);
2771:   MPI_Allreduce(&locFlag,   flag,    1, MPI_INT, MPI_LOR, mesh->comm);
2772:   MPI_Allreduce(&locRetVal, &retVal, 1, MPI_INT, MPI_LOR, mesh->comm);
2773:   PetscFunctionReturn(retVal);
2774: }

2776: static struct _MeshOps MeshOps = {/* Generic Operations */
2777:                                   PETSC_NULL/* MeshSetup */,
2778:                                   PETSC_NULL/* MeshSetFromOptions */,
2779:                                   MeshView_Triangular_2D,
2780:                                   PETSC_NULL/* MeshCopy */,
2781:                                   PETSC_NULL/* MeshDuplicate */,
2782:                                   MeshDestroy_Triangular_2D,
2783:                                   /* Mesh-Specific Operations */
2784:                                   MeshPartition_Triangular_2D,
2785:                                   MeshCoarsen_Triangular_2D,
2786:                                   MeshRefine_Triangular_2D,
2787:                                   MeshResetNodes_Triangular_2D,
2788:                                   MeshSaveMesh_Triangular_2D,
2789:                                   MeshRestoreMesh_Triangular_2D,
2790:                                   /* Mesh Query Functions */
2791:                                   MeshUpdateBoundingBox_Triangular_2D,
2792:                                   MeshIsDistorted_Triangular_2D,
2793:                                   /* Mesh Boundary Query Functions */
2794:                                   MeshGetBoundarySize_Triangular_2D,
2795:                                   MeshGetBoundaryIndex_Triangular_2D,
2796:                                   MeshGetBoundaryStart_Triangular_2D,
2797:                                   MeshGetBoundaryNext_Triangular_2D,
2798:                                   MeshGetActiveBoundary_Triangular_2D,
2799:                                   /* Mesh Node Query Functions */
2800:                                   MeshGetNodeBoundary_Triangular_2D,
2801:                                   MeshNodeIsVertex_Triangular_2D,
2802:                                   MeshGetNodeCoords_Triangular_2D,
2803:                                   MeshSetNodeCoords_Triangular_2D,
2804:                                   MeshGetNodeCoordsSaved_Triangular_2D,
2805:                                   MeshNearestNode_Triangular_2D,
2806:                                   MeshGetNodeSupport_Triangular_2D,
2807:                                   /* Mesh Element Query Functions */
2808:                                   MeshGetElemNeighbor_Triangular_2D,
2809:                                   MeshLocatePoint_Triangular_2D,
2810:                                   /* Mesh Embedding Query Functions */
2811:                                   MeshGetNodeFromElement_Triangular_2D,
2812:                                   MeshGetNodeFromEdge_Triangular_2D,
2813:                                   /* CSR Support Functions */
2814:                                   MeshCreateLocalCSR_Triangular_2D,
2815:                                   MeshCreateFullCSR_Triangular_2D,
2816:                                   MeshCreateDualCSR_Triangular_2D};

2818: #undef  __FUNCT__
2820: /*@
2821:   MeshCreateTriangular2DCSR - Creates a basic two dimensional unstructured grid from an existing CSR representation.

2823:   Collective on MPI_Comm

2825:   Input Parameters:
2826: + comm     - The communicator that shares the grid
2827: . numNodes - The number of mesh nodes, here identical to vertices
2828: . numFaces - The number of elements in the mesh
2829: . nodes    - The coordinates for each node
2830: . offsets  - The offset of each node's adjacency list
2831: . adj      - The adjacency list for the mesh
2832: - faces    - The list of nodes for each element

2834:   Output Parameter:
2835: . mesh    - The mesh

2837:   Level: beginner

2839: .keywords unstructured mesh
2840: .seealso MeshCreateTriangular2D(), MeshCreateTriangular2DPeriodic()
2841: @*/
2842: int MeshCreateTriangular2DCSR(MPI_Comm comm, int numNodes, int numFaces, double *nodes, int *offsets, int *adj, int *faces, Mesh *mesh)
2843: {
2844:   ParameterDict   dict;
2845:   MeshBoundary2D *bdCtx;
2846:   int             ierr;

2849:   /* REWORK: We are really stretching this interface, but I need this yesterday */
2850:   PetscMalloc(sizeof(MeshBoundary2D), &bdCtx);
2851:   bdCtx->numBd       = 0;
2852:   bdCtx->numVertices = numNodes;
2853:   bdCtx->vertices    = nodes;
2854:   bdCtx->markers     = offsets;
2855:   bdCtx->segMarkers  = adj;
2856:   bdCtx->numSegments = numFaces;
2857:   bdCtx->segments    = faces;
2858:   bdCtx->numHoles    = 0;
2859:   bdCtx->holes       = PETSC_NULL;
2860:   MeshCreate(comm, mesh);
2861:   ParameterDictCreate(comm, &dict);
2862:   ParameterDictSetObject(dict, "bdCtx", bdCtx);
2863:   ParameterDictSetInteger(dict, "numLocNodes", 3);
2864:   PetscObjectSetParameterDict((PetscObject) *mesh, dict);
2865:   ParameterDictDestroy(dict);
2866:   MeshSetDimension(*mesh, 2);
2867:   /* Setup special mechanisms */
2868:   PetscObjectComposeFunction((PetscObject) *mesh, "MeshTriangular2D_Create_C", "MeshCreate_CSR",
2869:                                     (void (*)(void)) MeshCreate_CSR);
2870: 
2871:   PetscObjectComposeFunction((PetscObject) *mesh, "MeshTriangular2D_Refine_C", "MeshRefine_CSR",
2872:                                     (void (*)(void)) MeshRefine_CSR);
2873: 
2874:   MeshSetType(*mesh, MESH_TRIANGULAR_2D);

2876:   PetscObjectSetName((PetscObject) *mesh, "Constructed Mesh");
2877:   PetscFree(bdCtx);
2878:   return(0);
2879: }

2881: EXTERN_C_BEGIN
2882: #undef  __FUNCT__
2884: int MeshSerialize_Triangular_2D(MPI_Comm comm, Mesh *mesh, PetscViewer viewer, PetscTruth store)
2885: {
2886:   Mesh             m;
2887:   Mesh_Triangular *tri;
2888:   int              fd;
2889:   int              zero = 0;
2890:   int              one  = 1;
2891:   int              numNodes, numEdges, numFaces, numCorners, numBd, numBdNodes, numBdEdges, old;
2892:   int              ierr;

2895:   PetscViewerBinaryGetDescriptor(viewer, &fd);
2896:   if (store) {
2897:     m    = *mesh;
2898:     PetscBinaryWrite(fd, &m->highlightElement, 1,                   PETSC_INT,    0);
2899:     PetscBinaryWrite(fd, &m->maxDegree,        1,                   PETSC_INT,    0);
2900:     PetscBinaryWrite(fd, &m->startX,           1,                   PETSC_DOUBLE, 0);
2901:     PetscBinaryWrite(fd, &m->endX,             1,                   PETSC_DOUBLE, 0);
2902:     PetscBinaryWrite(fd, &m->sizeX,            1,                   PETSC_DOUBLE, 0);
2903:     PetscBinaryWrite(fd, &m->startY,           1,                   PETSC_DOUBLE, 0);
2904:     PetscBinaryWrite(fd, &m->endY,             1,                   PETSC_DOUBLE, 0);
2905:     PetscBinaryWrite(fd, &m->sizeY,            1,                   PETSC_DOUBLE, 0);
2906:     PetscBinaryWrite(fd, &m->startZ,           1,                   PETSC_DOUBLE, 0);
2907:     PetscBinaryWrite(fd, &m->endZ,             1,                   PETSC_DOUBLE, 0);
2908:     PetscBinaryWrite(fd, &m->sizeZ,            1,                   PETSC_DOUBLE, 0);
2909:     PetscBinaryWrite(fd, &m->locStartX,        1,                   PETSC_DOUBLE, 0);
2910:     PetscBinaryWrite(fd, &m->locEndX,          1,                   PETSC_DOUBLE, 0);
2911:     PetscBinaryWrite(fd, &m->locSizeX,         1,                   PETSC_DOUBLE, 0);
2912:     PetscBinaryWrite(fd, &m->locStartY,        1,                   PETSC_DOUBLE, 0);
2913:     PetscBinaryWrite(fd, &m->locEndY,          1,                   PETSC_DOUBLE, 0);
2914:     PetscBinaryWrite(fd, &m->locSizeY,         1,                   PETSC_DOUBLE, 0);
2915:     PetscBinaryWrite(fd, &m->locStartZ,        1,                   PETSC_DOUBLE, 0);
2916:     PetscBinaryWrite(fd, &m->locEndZ,          1,                   PETSC_DOUBLE, 0);
2917:     PetscBinaryWrite(fd, &m->locSizeZ,         1,                   PETSC_DOUBLE, 0);
2918:     PetscBinaryWrite(fd, &m->isPeriodic,       1,                   PETSC_INT,    0);
2919:     PetscBinaryWrite(fd, &m->isPeriodicDim,    3,                   PETSC_INT,    0);

2921:     PetscBinaryWrite(fd, &m->numBd,            1,                   PETSC_INT,    0);
2922:     PetscBinaryWrite(fd, &m->numVertices,      1,                   PETSC_INT,    0);
2923:     PetscBinaryWrite(fd, &m->numNodes,         1,                   PETSC_INT,    0);
2924:     PetscBinaryWrite(fd, &m->numBdNodes,       1,                   PETSC_INT,    0);
2925:     PetscBinaryWrite(fd, &m->numEdges,         1,                   PETSC_INT,    0);
2926:     PetscBinaryWrite(fd, &m->numBdEdges,       1,                   PETSC_INT,    0);
2927:     PetscBinaryWrite(fd, &m->numFaces,         1,                   PETSC_INT,    0);
2928:     PetscBinaryWrite(fd, &m->numBdFaces,       1,                   PETSC_INT,    0);
2929:     PetscBinaryWrite(fd, &m->numCells,         1,                   PETSC_INT,    0);
2930:     PetscBinaryWrite(fd, &m->numCorners,       1,                   PETSC_INT,    0);
2931:     PetscBinaryWrite(fd, &m->numCellCorners,   1,                   PETSC_INT,    0);
2932:     PetscBinaryWrite(fd, &m->numHoles,         1,                   PETSC_INT,    0);

2934:     PartitionSerialize(m, &m->part, viewer, store);

2936:     PartitionGetNumOverlapElements(m->part, &numFaces);
2937:     PartitionGetNumOverlapNodes(m->part, &numNodes);
2938:     numEdges   = m->numEdges;
2939:     numCorners = m->numCorners;
2940:     numBd      = m->numBd;
2941:     numBdNodes = m->numBdNodes;
2942:     numBdEdges = m->numBdEdges;

2944:     tri  = (Mesh_Triangular *) (*mesh)->data;
2945:     PetscBinaryWrite(fd,  tri->nodes,          numNodes*2,          PETSC_DOUBLE, 0);
2946:     if (tri->nodesOld == PETSC_NULL) {
2947:       PetscBinaryWrite(fd, &zero,              1,                   PETSC_INT,    0);
2948:     } else {
2949:       PetscBinaryWrite(fd, &one,               1,                   PETSC_INT,    0);
2950:       PetscBinaryWrite(fd,  tri->nodesOld,     numNodes*2,          PETSC_DOUBLE, 0);
2951:     }
2952:     PetscBinaryWrite(fd,  tri->markers,        numNodes,            PETSC_INT,    0);
2953:     PetscBinaryWrite(fd,  tri->degrees,        numNodes,            PETSC_INT,    0);
2954:     PetscBinaryWrite(fd,  tri->edges,          numEdges*2,          PETSC_INT,    0);
2955:     PetscBinaryWrite(fd,  tri->edgemarkers,    numEdges,            PETSC_INT,    0);
2956:     PetscBinaryWrite(fd,  tri->faces,          numFaces*numCorners, PETSC_INT,    0);
2957:     PetscBinaryWrite(fd,  tri->neighbors,      numFaces*3,          PETSC_INT,    0);
2958:     PetscBinaryWrite(fd,  tri->bdNodes,        numBdNodes,          PETSC_INT,    0);
2959:     PetscBinaryWrite(fd,  tri->bdEdges,        numBdEdges,          PETSC_INT,    0);
2960:     PetscBinaryWrite(fd,  tri->bdMarkers,      numBd,               PETSC_INT,    0);
2961:     PetscBinaryWrite(fd,  tri->bdBegin,        numBd+1,             PETSC_INT,    0);
2962:     PetscBinaryWrite(fd,  tri->bdEdgeBegin,    numBd+1,             PETSC_INT,    0);
2963:     PetscBinaryWrite(fd,  m->holes,            m->numHoles*2,       PETSC_DOUBLE, 0);

2965:     PetscBinaryWrite(fd, &m->maxAspectRatio,   1,                   PETSC_DOUBLE, 0);
2966:   } else {
2967:     /* Create the mesh context */
2968:     MeshCreate(comm, &m);
2969:     PetscNew(Mesh_Triangular, &tri);
2970:     PetscLogObjectMemory(m, sizeof(Mesh_Triangular));
2971:     PetscMemcpy(m->ops, &MeshOps, sizeof(struct _MeshOps));
2972:     PetscStrallocpy(MESH_TRIANGULAR_2D, &m->type_name);
2973:     m->data = (void *) tri;
2974:     m->dim  = 2;

2976:     PetscBinaryRead(fd, &m->highlightElement, 1, PETSC_INT);
2977:     PetscBinaryRead(fd, &m->maxDegree,        1, PETSC_INT);
2978:     PetscMalloc(m->maxDegree * sizeof(int), &m->support);

2980:     PetscBinaryRead(fd, &m->startX,           1, PETSC_DOUBLE);
2981:     PetscBinaryRead(fd, &m->endX,             1, PETSC_DOUBLE);
2982:     PetscBinaryRead(fd, &m->sizeX,            1, PETSC_DOUBLE);
2983:     PetscBinaryRead(fd, &m->startY,           1, PETSC_DOUBLE);
2984:     PetscBinaryRead(fd, &m->endY,             1, PETSC_DOUBLE);
2985:     PetscBinaryRead(fd, &m->sizeY,            1, PETSC_DOUBLE);
2986:     PetscBinaryRead(fd, &m->startZ,           1, PETSC_DOUBLE);
2987:     PetscBinaryRead(fd, &m->endZ,             1, PETSC_DOUBLE);
2988:     PetscBinaryRead(fd, &m->sizeZ,            1, PETSC_DOUBLE);
2989:     PetscBinaryRead(fd, &m->locStartX,        1, PETSC_DOUBLE);
2990:     PetscBinaryRead(fd, &m->locEndX,          1, PETSC_DOUBLE);
2991:     PetscBinaryRead(fd, &m->locSizeX,         1, PETSC_DOUBLE);
2992:     PetscBinaryRead(fd, &m->locStartY,        1, PETSC_DOUBLE);
2993:     PetscBinaryRead(fd, &m->locEndY,          1, PETSC_DOUBLE);
2994:     PetscBinaryRead(fd, &m->locSizeY,         1, PETSC_DOUBLE);
2995:     PetscBinaryRead(fd, &m->locStartZ,        1, PETSC_DOUBLE);
2996:     PetscBinaryRead(fd, &m->locEndZ,          1, PETSC_DOUBLE);
2997:     PetscBinaryRead(fd, &m->locSizeZ,         1, PETSC_DOUBLE);
2998:     PetscBinaryRead(fd, &m->isPeriodic,       1, PETSC_INT);
2999:     PetscBinaryRead(fd, &m->isPeriodicDim,    3, PETSC_INT);

3001:     m->activeBd       = -1;
3002:     m->activeBdOld    = -1;
3003:     m->activeBdNode   = -1;
3004:     tri->areas        = PETSC_NULL;
3005:     tri->aspectRatios = PETSC_NULL;

3007:     PetscBinaryRead(fd, &m->numBd,            1,                   PETSC_INT);
3008:     PetscBinaryRead(fd, &m->numVertices,      1,                   PETSC_INT);
3009:     PetscBinaryRead(fd, &m->numNodes,         1,                   PETSC_INT);
3010:     PetscBinaryRead(fd, &m->numBdNodes,       1,                   PETSC_INT);
3011:     PetscBinaryRead(fd, &m->numEdges,         1,                   PETSC_INT);
3012:     PetscBinaryRead(fd, &m->numBdEdges,       1,                   PETSC_INT);
3013:     PetscBinaryRead(fd, &m->numFaces,         1,                   PETSC_INT);
3014:     PetscBinaryRead(fd, &m->numBdFaces,       1,                   PETSC_INT);
3015:     PetscBinaryRead(fd, &m->numCells,         1,                   PETSC_INT);
3016:     PetscBinaryRead(fd, &m->numCorners,       1,                   PETSC_INT);
3017:     PetscBinaryRead(fd, &m->numCellCorners,   1,                   PETSC_INT);
3018:     PetscBinaryRead(fd, &m->numHoles,         1,                   PETSC_INT);

3020:     PartitionSerialize(m, &m->part, viewer, store);
3021:     PetscLogObjectParent(m, m->part);

3023:     PartitionGetNumOverlapElements(m->part, &numFaces);
3024:     PartitionGetNumOverlapNodes(m->part, &numNodes);
3025:     numEdges   = m->numEdges;
3026:     numCorners = m->numCorners;
3027:     numBd      = m->numBd;
3028:     numBdNodes = m->numBdNodes;
3029:     numBdEdges = m->numBdEdges;
3030:     PetscMalloc(numNodes*2          * sizeof(double), &tri->nodes);
3031:     PetscMalloc(numNodes            * sizeof(int),    &tri->markers);
3032:     PetscMalloc(numNodes            * sizeof(int),    &tri->degrees);
3033:     PetscMalloc(numEdges*2          * sizeof(int),    &tri->edges);
3034:     PetscMalloc(numEdges            * sizeof(int),    &tri->edgemarkers);
3035:     PetscMalloc(numFaces*numCorners * sizeof(int),    &tri->faces);
3036:     PetscMalloc(numFaces*3          * sizeof(int),    &tri->neighbors);
3037:     PetscMalloc(numBdNodes          * sizeof(int),    &tri->bdNodes);
3038:     PetscMalloc(numBdEdges          * sizeof(int),    &tri->bdEdges);
3039:     PetscMalloc(numBd               * sizeof(int),    &tri->bdMarkers);
3040:     PetscMalloc((numBd+1)           * sizeof(int),    &tri->bdBegin);
3041:     PetscMalloc((numBd+1)           * sizeof(int),    &tri->bdEdgeBegin);
3042:     if (m->numHoles > 0) {
3043:       PetscMalloc(m->numHoles*2     * sizeof(double), &m->holes);
3044:     } else {
3045:       m->holes = PETSC_NULL;
3046:     }
3047:     PetscLogObjectMemory(m, (numNodes*2 + m->numHoles*2) * sizeof(double) +
3048:                          (numNodes*2 + numEdges*3 + numFaces*(numCorners+1) + numBdNodes + numBdEdges + numBd*3+2)*sizeof(int));
3049:     PetscBinaryRead(fd,  tri->nodes,          numNodes*2,          PETSC_DOUBLE);
3050:     PetscBinaryRead(fd, &old,                 1,                   PETSC_INT);
3051:     if (old) {
3052:       PetscMalloc(numNodes*2 * sizeof(double), &tri->nodesOld);
3053:       PetscLogObjectMemory(m, numNodes*2 * sizeof(double));
3054:       PetscBinaryRead(fd,  tri->nodesOld,     numNodes*2,          PETSC_DOUBLE);
3055:     } else {
3056:       tri->nodesOld = PETSC_NULL;
3057:     }
3058:     PetscBinaryRead(fd,  tri->markers,        numNodes,            PETSC_INT);
3059:     PetscBinaryRead(fd,  tri->degrees,        numNodes,            PETSC_INT);
3060:     PetscBinaryRead(fd,  tri->edges,          numEdges*2,          PETSC_INT);
3061:     PetscBinaryRead(fd,  tri->edgemarkers,    numEdges,            PETSC_INT);
3062:     PetscBinaryRead(fd,  tri->faces,          numFaces*numCorners, PETSC_INT);
3063:     PetscBinaryRead(fd,  tri->neighbors,      numFaces*3,          PETSC_INT);
3064:     PetscBinaryRead(fd,  tri->bdNodes,        numBdNodes,          PETSC_INT);
3065:     PetscBinaryRead(fd,  tri->bdEdges,        numBdEdges,          PETSC_INT);
3066:     PetscBinaryRead(fd,  tri->bdMarkers,      numBd,               PETSC_INT);
3067:     PetscBinaryRead(fd,  tri->bdBegin,        numBd+1,             PETSC_INT);
3068:     PetscBinaryRead(fd,  tri->bdEdgeBegin,    numBd+1,             PETSC_INT);
3069:     PetscBinaryRead(fd,  m->holes,            m->numHoles*m->dim,  PETSC_DOUBLE);

3071:     PetscBinaryRead(fd, &m->maxAspectRatio,   1,                   PETSC_DOUBLE);

3073: #ifdef PETSC_USE_BOPT_g
3074:     /* Check mesh integrity */
3075:     MeshDebug_Triangular_2D(m, PETSC_TRUE);
3076: #endif
3077:     *mesh = m;
3078:   }

3080:   return(0);
3081: }
3082: EXTERN_C_END

3084: EXTERN_C_BEGIN
3085: #undef  __FUNCT__
3087: int MeshCreate_Triangular_2D(Mesh mesh) {
3088:   Mesh_Triangular *tri;
3089:   int            (*f)(MeshBoundary2D *, int, Mesh);
3090:   MPI_Comm         comm;
3091:   PetscTruth       opt;
3092:   int              ierr;

3095:   PetscNew(Mesh_Triangular, &tri);
3096:   PetscMemcpy(mesh->ops, &MeshOps, sizeof(struct _MeshOps));
3097:   mesh->data = (void *) tri;
3098:   PetscStrallocpy(MESH_SER_TRIANGULAR_2D_BINARY, &mesh->serialize_name);
3099:   mesh->dim  = 2;

3101:   /* Setup the periodic domain */
3102:   if (mesh->isPeriodic == PETSC_TRUE) {
3103:     MeshBoundary2DSetBoundingBox(mesh, mesh->bdCtx);
3104:   }

3106:   /* Setup default mechanisms */
3107:   PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", (PetscVoidFunction) &f);
3108:   if (f == PETSC_NULL) {
3109: #ifdef PETSC_HAVE_TRIANGLE
3110:     if (mesh->isPeriodic == PETSC_TRUE) {
3111:       PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", "MeshCreate_Periodic",
3112:                                         (void (*)(void)) MeshCreate_Periodic);
3113: 
3114:       PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", "MeshRefine_Periodic",
3115:                                         (void (*)(void)) MeshRefine_Periodic);
3116: 
3117:     } else {
3118:       PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", "MeshCreate_Triangle",
3119:                                         (void (*)(void)) MeshCreate_Triangle);
3120: 
3121:       PetscObjectComposeFunction((PetscObject) mesh, "MeshTriangular2D_Refine_C", "MeshRefine_Triangle",
3122:                                         (void (*)(void)) MeshRefine_Triangle);
3123:     }
3124: #endif
3125:   }

3127:   MeshCheckBoundary_Triangular_2D(mesh);

3129:   PetscObjectQueryFunction((PetscObject) mesh, "MeshTriangular2D_Create_C", (PetscVoidFunction) &f);
3130:   if (!f) SETERRQ(1,"Unable to find mesh generation function");
3131:   (*f)(mesh->bdCtx, mesh->numCorners, mesh);

3133:   /* Calculate maximum degree of vertices */
3134:   MeshSetupSupport_Triangular_2D(mesh);

3136:   /* Construct derived and boundary information */
3137:   MeshSetupBoundary_Triangular_2D(mesh);

3139:   /* Reorder nodes before distributing mesh */
3140:   PetscOptionsHasName(mesh->prefix, "-mesh_reorder", &opt);
3141:   if (opt == PETSC_TRUE) {
3142:     /* MUST FIX: Since we have duplicated the whole mesh, we may impersonate a serial mesh */
3143:     comm       = mesh->comm;
3144:     mesh->comm = PETSC_COMM_SELF;
3145:     MeshGetOrdering(mesh, MESH_ORDER_TRIANGULAR_2D_RCM, &mesh->nodeOrdering);
3146:     mesh->comm = comm;

3148:     /* Permute arrays implicitly numbered by node numbers */
3149:     AOApplicationToPetscPermuteReal(mesh->nodeOrdering, mesh->dim, tri->nodes);
3150:     AOApplicationToPetscPermuteInt(mesh->nodeOrdering, 1, tri->markers);
3151:     AOApplicationToPetscPermuteInt(mesh->nodeOrdering, 1, tri->degrees);
3152:     /* Renumber arrays dependent on the canonical node numbering */
3153:     AOApplicationToPetsc(mesh->nodeOrdering, mesh->numEdges*2,                tri->edges);
3154:     AOApplicationToPetsc(mesh->nodeOrdering, mesh->numFaces*mesh->numCorners, tri->faces);
3155:     AOApplicationToPetsc(mesh->nodeOrdering, mesh->numBdNodes,                tri->bdNodes);
3156:   }

3158:   /* Initialize partition */
3159:   MeshPartition(mesh);

3161:   /* Initialize save space */
3162:   tri->nodesOld = PETSC_NULL;

3164:   /* Reset the midnodes */
3165:   if (mesh->isPeriodic == PETSC_TRUE) {
3166:     MeshResetNodes(mesh, PETSC_TRUE);
3167:   }

3169:   return(0);
3170: }
3171: EXTERN_C_END