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, °ree, &support);
353: /* Check node degree */
354: PetscMalloc(numProcs * sizeof(int), °rees);
355: MPI_Allgather(°ree, 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, °ree, &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, °ree, &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, °ree, &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, °ree, &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, °ree, &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