Actual source code: complex.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/compleximpl.h> /*I "petscdmcomplex.h" I*/
5: PetscErrorCode DMComplexView_Ascii(DM dm, PetscViewer viewer)
6: {
7: DM_Complex *mesh = (DM_Complex *) dm->data;
8: PetscViewerFormat format;
9: PetscErrorCode ierr;
12: PetscViewerGetFormat(viewer, &format);
13: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
14: const char *name;
15: PetscInt maxConeSize, maxSupportSize;
16: PetscInt pStart, pEnd, p;
17: PetscMPIInt rank;
18: PetscBool hasLabel;
20: MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
21: PetscObjectGetName((PetscObject) dm, &name);
22: DMComplexGetChart(dm, &pStart, &pEnd);
23: DMComplexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
24: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
25: PetscViewerASCIIPrintf(viewer, "Mesh '%s':\n", name);
26: PetscViewerASCIISynchronizedPrintf(viewer, "Max sizes cone: %D support: %D\n", maxConeSize, maxSupportSize);
27: PetscViewerASCIIPrintf(viewer, "orientation is missing\n", name);
28: PetscViewerASCIIPrintf(viewer, "cap --> base:\n", name);
29: for(p = pStart; p < pEnd; ++p) {
30: PetscInt dof, off, s;
32: PetscSectionGetDof(mesh->supportSection, p, &dof);
33: PetscSectionGetOffset(mesh->supportSection, p, &off);
34: for(s = off; s < off+dof; ++s) {
35: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D ----> %D\n", rank, p, mesh->supports[s]);
36: }
37: }
38: PetscViewerFlush(viewer);
39: PetscViewerASCIIPrintf(viewer, "base <-- cap:\n", name);
40: for(p = pStart; p < pEnd; ++p) {
41: PetscInt dof, off, c;
43: PetscSectionGetDof(mesh->coneSection, p, &dof);
44: PetscSectionGetOffset(mesh->coneSection, p, &off);
45: for(c = off; c < off+dof; ++c) {
46: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D <---- %D (%D)\n", rank, p, mesh->cones[c], mesh->coneOrientations[c]);
47: }
48: }
49: PetscViewerFlush(viewer);
50: PetscSectionVecView(mesh->coordSection, mesh->coordinates, viewer);
51: DMComplexHasLabel(dm, "marker", &hasLabel);
52: if (hasLabel) {
53: const char *name = "marker";
54: IS ids;
55: const PetscInt *markers;
56: PetscInt num, i;
58: PetscViewerASCIIPrintf(viewer, "Label '%s':\n", name);
59: DMComplexGetLabelIdIS(dm, name, &ids);
60: ISGetSize(ids, &num);
61: ISGetIndices(ids, &markers);
62: for(i = 0; i < num; ++i) {
63: IS pIS;
64: const PetscInt *points;
65: PetscInt size, p;
67: DMComplexGetStratumIS(dm, name, markers[i], &pIS);
68: ISGetSize(pIS, &size);
69: ISGetIndices(pIS, &points);
70: for(p = 0; p < size; ++p) {
71: PetscViewerASCIISynchronizedPrintf(viewer, "[%D]: %D (%D)\n", rank, points[p], markers[i]);
72: }
73: ISRestoreIndices(pIS, &points);
74: ISDestroy(&pIS);
75: }
76: ISRestoreIndices(ids, &markers);
77: ISDestroy(&ids);
78: }
79: PetscViewerFlush(viewer);
80: } else if (format == PETSC_VIEWER_ASCII_LATEX) {
81: const char *name;
82: const char *colors[3] = {"red", "blue", "green"};
83: const int numColors = 3;
84: PetscScalar *coords;
85: PetscInt cStart, cEnd, c, vStart, vEnd, v, p;
86: PetscMPIInt rank, size;
88: MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
89: MPI_Comm_size(((PetscObject) dm)->comm, &size);
90: PetscObjectGetName((PetscObject) dm, &name);
91: PetscViewerASCIISynchronizedAllow(viewer, PETSC_TRUE);
92: PetscViewerASCIIPrintf(viewer, "\\documentclass{beamer}\n\n\
93: \\usepackage{tikz}\n\
94: \\usepackage{pgflibraryshapes}\n\
95: \\usetikzlibrary{backgrounds}\n\
96: \\usetikzlibrary{arrows}\n\
97: \\newenvironment{changemargin}[2]{%%\n\
98: \\begin{list}{}{%%\n\
99: \\setlength{\\topsep}{0pt}%%\n\
100: \\setlength{\\leftmargin}{#1}%%\n\
101: \\setlength{\\rightmargin}{#2}%%\n\
102: \\setlength{\\listparindent}{\\parindent}%%\n\
103: \\setlength{\\itemindent}{\\parindent}%%\n\
104: \\setlength{\\parsep}{\\parskip}%%\n\
105: }%%\n\
106: \\item[]}{\\end{list}}\n\n\
107: \\begin{document}\n\
108: \\begin{frame}{%s}\n\
109: \\begin{changemargin}{-1cm}{0cm}\n\
110: \\begin{center}\n\
111: \\begin{tikzpicture}[scale = 5.00,font=\\fontsize{8}{8}\\selectfont]\n", name);
112: /* Plot vertices */
113: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
114: VecGetArray(mesh->coordinates, &coords);
115: PetscViewerASCIIPrintf(viewer, "\\path\n");
116: for(v = vStart; v < vEnd; ++v) {
117: PetscInt off, dof, d;
119: PetscSectionGetDof(mesh->coordSection, v, &dof);
120: PetscSectionGetOffset(mesh->coordSection, v, &off);
121: PetscViewerASCIISynchronizedPrintf(viewer, "(");
122: for(d = 0; d < dof; ++d) {
123: if (d > 0) {PetscViewerASCIISynchronizedPrintf(viewer, ",");}
124: PetscViewerASCIISynchronizedPrintf(viewer, "%G", PetscRealPart(coords[off+d]));
125: }
126: PetscViewerASCIISynchronizedPrintf(viewer, ") node(%D_%D) [draw,shape=circle,color=%s] {%D} --\n", v, rank, colors[rank%numColors], v);
127: }
128: VecRestoreArray(mesh->coordinates, &coords);
129: PetscViewerFlush(viewer);
130: PetscViewerASCIIPrintf(viewer, "(0,0);\n");
131: /* Plot cells */
132: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
133: for(c = cStart; c < cEnd; ++c) {
134: PetscInt *closure = PETSC_NULL;
135: PetscInt closureSize, firstPoint = -1;
137: DMComplexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
138: PetscViewerASCIISynchronizedPrintf(viewer, "\\draw[color=%s] ", colors[rank%numColors]);
139: for(p = 0; p < closureSize*2; p += 2) {
140: const PetscInt point = closure[p];
142: if ((point < vStart) || (point >= vEnd)) continue;
143: if (firstPoint >= 0) {PetscViewerASCIISynchronizedPrintf(viewer, " -- ");}
144: PetscViewerASCIISynchronizedPrintf(viewer, "(%D_%D)", point, rank);
145: if (firstPoint < 0) firstPoint = point;
146: }
147: /* Why doesn't this work? PetscViewerASCIISynchronizedPrintf(viewer, " -- cycle;\n"); */
148: PetscViewerASCIISynchronizedPrintf(viewer, " -- (%D_%D);\n", firstPoint, rank);
149: }
150: PetscViewerFlush(viewer);
151: PetscViewerASCIIPrintf(viewer, "\\end{tikzpicture}\n\\end{center}\n");
152: PetscViewerASCIIPrintf(viewer, "Mesh for processes ");
153: for(p = 0; p < size; ++p) {
154: if (p == size-1) {
155: PetscViewerASCIIPrintf(viewer, ", and ", colors[p%numColors], p);
156: } else if (p > 0) {
157: PetscViewerASCIIPrintf(viewer, ", ", colors[p%numColors], p);
158: }
159: PetscViewerASCIIPrintf(viewer, "{\\textcolor{%s}%D}", colors[p%numColors], p);
160: }
161: PetscViewerASCIIPrintf(viewer, ".\n");
162: PetscViewerASCIIPrintf(viewer, "\\end{changemargin}\n\
163: \\end{frame}\n\
164: \\end{document}\n", name);
165: } else {
166: MPI_Comm comm = ((PetscObject) dm)->comm;
167: PetscInt *sizes;
168: PetscInt locDepth, depth, dim, d;
169: PetscInt pStart, pEnd, p;
170: PetscMPIInt size;
172: MPI_Comm_size(comm, &size);
173: DMComplexGetDimension(dm, &dim);
174: PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dim);
175: DMComplexGetDepth(dm, &locDepth);
176: MPI_Allreduce(&locDepth, &depth, 1, MPIU_INT, MPI_MAX, comm);
177: PetscMalloc(size * sizeof(PetscInt), &sizes);
178: if (depth == 1) {
179: DMComplexGetDepthStratum(dm, 0, &pStart, &pEnd);
180: pEnd = pEnd - pStart;
181: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
182: PetscViewerASCIIPrintf(viewer, " %D-cells:", 0);
183: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
184: PetscViewerASCIIPrintf(viewer, "\n");
185: DMComplexGetHeightStratum(dm, 0, &pStart, &pEnd);
186: pEnd = pEnd - pStart;
187: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
188: PetscViewerASCIIPrintf(viewer, " %D-cells:", dim);
189: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
190: PetscViewerASCIIPrintf(viewer, "\n");
191: } else {
192: for(d = 0; d <= dim; d++) {
193: DMComplexGetDepthStratum(dm, d, &pStart, &pEnd);
194: pEnd = pEnd - pStart;
195: MPI_Gather(&pEnd, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
196: PetscViewerASCIIPrintf(viewer, " %D-cells:", d);
197: for(p = 0; p < size; ++p) {PetscViewerASCIIPrintf(viewer, " %D", sizes[p]);}
198: PetscViewerASCIIPrintf(viewer, "\n");
199: }
200: }
201: PetscFree(sizes);
202: }
203: return(0);
204: }
208: PetscErrorCode DMView_Complex(DM dm, PetscViewer viewer)
209: {
210: PetscBool iascii, isbinary;
216: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
217: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
218: if (iascii) {
219: DMComplexView_Ascii(dm, viewer);
220: #if 0
221: } else if (isbinary) {
222: DMComplexView_Binary(dm, viewer);
223: #endif
224: } else SETERRQ1(((PetscObject)viewer)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
225: return(0);
226: }
230: PetscErrorCode DMDestroy_Complex(DM dm)
231: {
232: DM_Complex *mesh = (DM_Complex *) dm->data;
233: SieveLabel next = mesh->labels;
237: PetscSectionDestroy(&mesh->coneSection);
238: PetscFree(mesh->cones);
239: PetscFree(mesh->coneOrientations);
240: PetscSectionDestroy(&mesh->supportSection);
241: PetscFree(mesh->supports);
242: PetscSectionDestroy(&mesh->coordSection);
243: VecDestroy(&mesh->coordinates);
244: PetscFree2(mesh->meetTmpA,mesh->meetTmpB);
245: PetscFree2(mesh->joinTmpA,mesh->joinTmpB);
246: PetscFree2(mesh->closureTmpA,mesh->closureTmpB);
247: PetscFree(mesh->facesTmp);
248: while(next) {
249: SieveLabel tmp;
251: PetscFree(next->name);
252: PetscFree(next->stratumValues);
253: PetscFree(next->stratumOffsets);
254: PetscFree(next->stratumSizes);
255: PetscFree(next->points);
256: tmp = next->next;
257: PetscFree(next);
258: next = tmp;
259: }
260: return(0);
261: }
265: PetscErrorCode DMComplexGetAdjacencySingleLevel_Private(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
266: {
267: const PetscInt *support = PETSC_NULL;
268: PetscInt numAdj = 0, maxAdjSize = *adjSize, supportSize, s;
269: PetscErrorCode ierr;
272: if (useClosure) {
273: DMComplexGetConeSize(dm, p, &supportSize);
274: DMComplexGetCone(dm, p, &support);
275: for(s = 0; s < supportSize; ++s) {
276: const PetscInt *cone = PETSC_NULL;
277: PetscInt coneSize, c, q;
279: DMComplexGetSupportSize(dm, support[s], &coneSize);
280: DMComplexGetSupport(dm, support[s], &cone);
281: for(c = 0; c < coneSize; ++c) {
282: for(q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
283: if (cone[c] == adj[q]) break;
284: }
285: if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
286: }
287: }
288: } else {
289: DMComplexGetSupportSize(dm, p, &supportSize);
290: DMComplexGetSupport(dm, p, &support);
291: for(s = 0; s < supportSize; ++s) {
292: const PetscInt *cone = PETSC_NULL;
293: PetscInt coneSize, c, q;
295: DMComplexGetConeSize(dm, support[s], &coneSize);
296: DMComplexGetCone(dm, support[s], &cone);
297: for(c = 0; c < coneSize; ++c) {
298: for(q = 0; q < numAdj || (adj[numAdj++] = cone[c],0); ++q) {
299: if (cone[c] == adj[q]) break;
300: }
301: if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
302: }
303: }
304: }
305: *adjSize = numAdj;
306: return(0);
307: }
311: PetscErrorCode DMComplexGetAdjacency_Private(DM dm, PetscInt p, PetscBool useClosure, const PetscInt *tmpClosure, PetscInt *adjSize, PetscInt adj[])
312: {
313: const PetscInt *star = tmpClosure;
314: PetscInt numAdj = 0, maxAdjSize = *adjSize, starSize, s;
315: PetscErrorCode ierr;
318: DMComplexGetTransitiveClosure(dm, p, useClosure, &starSize, (PetscInt **) &star);
319: for(s = 2; s < starSize*2; s += 2) {
320: const PetscInt *closure = PETSC_NULL;
321: PetscInt closureSize, c, q;
323: DMComplexGetTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt **) &closure);
324: for(c = 0; c < closureSize*2; c += 2) {
325: for(q = 0; q < numAdj || (adj[numAdj++] = closure[c],0); ++q) {
326: if (closure[c] == adj[q]) break;
327: }
328: if (numAdj > maxAdjSize) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid mesh exceeded adjacency allocation (%D)", maxAdjSize);
329: }
330: }
331: *adjSize = numAdj;
332: return(0);
333: }
337: PetscErrorCode DMComplexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
338: {
339: DM_Complex *mesh = (DM_Complex *) dm->data;
340: MPI_Comm comm = ((PetscObject) dm)->comm;
341: PetscSF sf = dm->sf, sfDof, sfAdj;
342: PetscSection leafSectionAdj, rootSectionAdj, sectionAdj;
343: PetscInt nleaves, l, p;
344: const PetscInt *leaves;
345: const PetscSFNode *remotes;
346: PetscInt pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
347: PetscInt *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols;
348: PetscInt depth, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
349: PetscLayout rLayout;
350: PetscInt locRows, rStart, rEnd, r;
351: PetscMPIInt size;
352: PetscBool debug = PETSC_FALSE;
353: PetscErrorCode ierr;
356: PetscOptionsGetBool(PETSC_NULL, "-dm_view_preallocation", &debug, PETSC_NULL);
357: MPI_Comm_size(comm, &size);
358: /* Create dof SF based on point SF */
359: if (debug) {
360: PetscPrintf(comm, "Input Section for Preallocation:\n");
361: PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
362: PetscPrintf(comm, "Input Global Section for Preallocation:\n");
363: PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);
364: PetscPrintf(comm, "Input SF for Preallocation:\n");
365: PetscSFView(sf, PETSC_NULL);
366: }
367: PetscSFCreateSectionSF(sf, section, PETSC_NULL, section, &sfDof);
368: if (debug) {
369: PetscPrintf(comm, "Dof SF for Preallocation:\n");
370: PetscSFView(sfDof, PETSC_NULL);
371: }
372: /* Create section for dof adjacency (dof ==> # adj dof) */
373: /* Two points p and q are adjacent if q \in closure(star(p)) */
374: PetscSectionGetChart(section, &pStart, &pEnd);
375: PetscSectionGetStorageSize(section, &numDof);
376: PetscSectionCreate(comm, &leafSectionAdj);
377: PetscSectionSetChart(leafSectionAdj, 0, numDof);
378: PetscSectionCreate(comm, &rootSectionAdj);
379: PetscSectionSetChart(rootSectionAdj, 0, numDof);
380: /* Fill in the ghost dofs on the interface */
381: PetscSFGetGraph(sf, PETSC_NULL, &nleaves, &leaves, &remotes);
382: DMComplexGetDepth(dm, &depth);
383: DMComplexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
384: maxClosureSize = (PetscInt) (2*PetscMax(pow((PetscReal) mesh->maxConeSize, depth)+1, pow((PetscReal) mesh->maxSupportSize, depth)+1));
385: maxAdjSize = (PetscInt) (pow((PetscReal) mesh->maxConeSize, depth)*pow((PetscReal) mesh->maxSupportSize, depth)+1);
386: PetscMalloc2(maxClosureSize,PetscInt,&tmpClosure,maxAdjSize,PetscInt,&tmpAdj);
388: /*
389: ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
390: 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
391: Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
392: 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
393: Create sfAdj connecting rootSectionAdj and leafSectionAdj
394: 3. Visit unowned points on interface, write adjacencies to adj
395: Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
396: 4. Visit owned points on interface, write adjacencies to rootAdj
397: Remove redundancy in rootAdj
398: ** The last two traversals use transitive closure
399: 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
400: Allocate memory addressed by sectionAdj (cols)
401: 6. Visit all owned points in the subdomain, insert dof adjacencies into cols
402: ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
403: */
405: for(l = 0; l < nleaves; ++l) {
406: PetscInt dof, off, d, q;
407: PetscInt p = leaves[l], numAdj = maxAdjSize;
409: PetscSectionGetDof(section, p, &dof);
410: PetscSectionGetOffset(section, p, &off);
411: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
412: for(q = 0; q < numAdj; ++q) {
413: PetscInt ndof, ncdof;
415: PetscSectionGetDof(section, tmpAdj[q], &ndof);
416: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
417: for(d = off; d < off+dof; ++d) {
418: PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);
419: }
420: }
421: }
422: PetscSectionSetUp(leafSectionAdj);
423: if (debug) {
424: PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
425: PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
426: }
427: /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
428: if (size > 1) {
429: PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
430: PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
431: }
432: if (debug) {
433: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
434: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
435: }
436: /* Add in local adjacency sizes for owned dofs on interface (roots) */
437: for(p = pStart; p < pEnd; ++p) {
438: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
440: PetscSectionGetDof(section, p, &dof);
441: PetscSectionGetOffset(section, p, &off);
442: if (!dof) continue;
443: PetscSectionGetDof(rootSectionAdj, off, &adof);
444: if (adof <= 0) continue;
445: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
446: for(q = 0; q < numAdj; ++q) {
447: PetscInt ndof, ncdof;
449: PetscSectionGetDof(section, tmpAdj[q], &ndof);
450: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
451: for(d = off; d < off+dof; ++d) {
452: PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);
453: }
454: }
455: }
456: PetscSectionSetUp(rootSectionAdj);
457: if (debug) {
458: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
459: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
460: }
461: /* Create adj SF based on dof SF */
462: PetscSFCreateSectionSF(sfDof, rootSectionAdj, PETSC_NULL, leafSectionAdj, &sfAdj);
463: if (debug) {
464: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
465: PetscSFView(sfAdj, PETSC_NULL);
466: }
467: PetscSFDestroy(&sfDof);
468: /* Create leaf adjacency */
469: PetscSectionSetUp(leafSectionAdj);
470: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
471: PetscMalloc(adjSize * sizeof(PetscInt), &adj);
472: PetscMemzero(adj, adjSize * sizeof(PetscInt));
473: for(l = 0; l < nleaves; ++l) {
474: PetscInt dof, off, d, q;
475: PetscInt p = leaves[l], numAdj = maxAdjSize;
477: PetscSectionGetDof(section, p, &dof);
478: PetscSectionGetOffset(section, p, &off);
479: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
480: for(d = off; d < off+dof; ++d) {
481: PetscInt aoff, i = 0;
483: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
484: for(q = 0; q < numAdj; ++q) {
485: PetscInt ndof, ncdof, ngoff, nd;
487: PetscSectionGetDof(section, tmpAdj[q], &ndof);
488: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
489: PetscSectionGetOffset(sectionGlobal, tmpAdj[q], &ngoff);
490: for(nd = 0; nd < ndof-ncdof; ++nd) {
491: adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
492: ++i;
493: }
494: }
495: }
496: }
497: /* Debugging */
498: if (debug) {
499: IS tmp;
500: PetscPrintf(comm, "Leaf adjacency indices\n");
501: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
502: ISView(tmp, PETSC_NULL);
503: }
504: /* Gather adjacenct indices to root */
505: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
506: PetscMalloc(adjSize * sizeof(PetscInt), &rootAdj);
507: for(r = 0; r < adjSize; ++r) {
508: rootAdj[r] = -1;
509: }
510: if (size > 1) {
511: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);
512: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);
513: }
514: PetscSFDestroy(&sfAdj);
515: PetscFree(adj);
516: /* Debugging */
517: if (debug) {
518: IS tmp;
519: PetscPrintf(comm, "Root adjacency indices after gather\n");
520: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
521: ISView(tmp, PETSC_NULL);
522: }
523: /* Add in local adjacency indices for owned dofs on interface (roots) */
524: for(p = pStart; p < pEnd; ++p) {
525: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
527: PetscSectionGetDof(section, p, &dof);
528: PetscSectionGetOffset(section, p, &off);
529: if (!dof) continue;
530: PetscSectionGetDof(rootSectionAdj, off, &adof);
531: if (adof <= 0) continue;
532: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
533: for(d = off; d < off+dof; ++d) {
534: PetscInt adof, aoff, i;
536: PetscSectionGetDof(rootSectionAdj, d, &adof);
537: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
538: i = adof-1;
539: for(q = 0; q < numAdj; ++q) {
540: PetscInt ndof, ncdof, ngoff, nd;
542: PetscSectionGetDof(section, tmpAdj[q], &ndof);
543: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
544: PetscSectionGetOffset(sectionGlobal, tmpAdj[q], &ngoff);
545: for(nd = 0; nd < ndof-ncdof; ++nd) {
546: rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd: ngoff+nd;
547: --i;
548: }
549: }
550: }
551: }
552: /* Debugging */
553: if (debug) {
554: IS tmp;
555: PetscPrintf(comm, "Root adjacency indices\n");
556: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
557: ISView(tmp, PETSC_NULL);
558: }
559: /* Compress indices */
560: PetscSectionSetUp(rootSectionAdj);
561: for(p = pStart; p < pEnd; ++p) {
562: PetscInt dof, cdof, off, d;
563: PetscInt adof, aoff;
565: PetscSectionGetDof(section, p, &dof);
566: PetscSectionGetConstraintDof(section, p, &cdof);
567: PetscSectionGetOffset(section, p, &off);
568: if (!dof) continue;
569: PetscSectionGetDof(rootSectionAdj, off, &adof);
570: if (adof <= 0) continue;
571: for(d = off; d < off+dof-cdof; ++d) {
572: PetscSectionGetDof(rootSectionAdj, d, &adof);
573: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
574: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
575: PetscSectionSetDof(rootSectionAdj, d, adof);
576: }
577: }
578: /* Debugging */
579: if (debug) {
580: IS tmp;
581: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
582: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
583: PetscPrintf(comm, "Root adjacency indices after compression\n");
584: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
585: ISView(tmp, PETSC_NULL);
586: }
587: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
588: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
589: PetscSectionCreate(comm, §ionAdj);
590: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
591: for(p = pStart; p < pEnd; ++p) {
592: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
593: PetscBool found = PETSC_TRUE;
595: PetscSectionGetDof(section, p, &dof);
596: PetscSectionGetConstraintDof(section, p, &cdof);
597: PetscSectionGetOffset(section, p, &off);
598: PetscSectionGetOffset(sectionGlobal, p, &goff);
599: for(d = 0; d < dof-cdof; ++d) {
600: PetscInt ldof, rdof;
602: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
603: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
604: if (ldof > 0) {
605: /* We do not own this point */
606: } else if (rdof > 0) {
607: PetscSectionSetDof(sectionAdj, goff+d, rdof);
608: } else {
609: found = PETSC_FALSE;
610: }
611: }
612: if (found) continue;
613: PetscSectionGetDof(section, p, &dof);
614: PetscSectionGetOffset(sectionGlobal, p, &goff);
615: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
616: for(q = 0; q < numAdj; ++q) {
617: PetscInt ndof, ncdof, noff;
619: PetscSectionGetDof(section, tmpAdj[q], &ndof);
620: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
621: PetscSectionGetOffset(section, tmpAdj[q], &noff);
622: for(d = goff; d < goff+dof-cdof; ++d) {
623: PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
624: }
625: }
626: }
627: PetscSectionSetUp(sectionAdj);
628: if (debug) {
629: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
630: PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
631: }
632: /* Get adjacent indices */
633: PetscSectionGetStorageSize(sectionAdj, &numCols);
634: PetscMalloc(numCols * sizeof(PetscInt), &cols);
635: for(p = pStart; p < pEnd; ++p) {
636: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
637: PetscBool found = PETSC_TRUE;
639: PetscSectionGetDof(section, p, &dof);
640: PetscSectionGetConstraintDof(section, p, &cdof);
641: PetscSectionGetOffset(section, p, &off);
642: PetscSectionGetOffset(sectionGlobal, p, &goff);
643: for(d = 0; d < dof-cdof; ++d) {
644: PetscInt ldof, rdof;
646: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
647: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
648: if (ldof > 0) {
649: /* We do not own this point */
650: } else if (rdof > 0) {
651: PetscInt aoff, roff;
653: PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
654: PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
655: PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
656: } else {
657: found = PETSC_FALSE;
658: }
659: }
660: if (found) continue;
661: DMComplexGetAdjacency_Private(dm, p, PETSC_FALSE, tmpClosure, &numAdj, tmpAdj);
662: for(d = goff; d < goff+dof-cdof; ++d) {
663: PetscInt adof, aoff, i = 0;
665: PetscSectionGetDof(sectionAdj, d, &adof);
666: PetscSectionGetOffset(sectionAdj, d, &aoff);
667: for(q = 0; q < numAdj; ++q) {
668: PetscInt ndof, ncdof, ngoff, nd;
669: PetscInt *ncind;
671: PetscSectionGetDof(section, tmpAdj[q], &ndof);
672: PetscSectionGetConstraintDof(section, tmpAdj[q], &ncdof);
673: PetscSectionGetConstraintIndices(section, tmpAdj[q], &ncind);
674: PetscSectionGetOffset(sectionGlobal, tmpAdj[q], &ngoff);
675: for(nd = 0; nd < ndof-ncdof; ++nd, ++i) {
676: cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd: ngoff+nd;
677: }
678: }
679: if (i != adof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %D != %D for dof %D (point %D)", i, adof, d, p);
680: }
681: }
682: PetscSectionDestroy(&leafSectionAdj);
683: PetscSectionDestroy(&rootSectionAdj);
684: PetscFree(rootAdj);
685: PetscFree2(tmpClosure, tmpAdj);
686: /* Debugging */
687: if (debug) {
688: IS tmp;
689: PetscPrintf(comm, "Column indices\n");
690: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
691: ISView(tmp, PETSC_NULL);
692: }
693: /* Create allocation vectors from adjacency graph */
694: MatGetLocalSize(A, &locRows, PETSC_NULL);
695: PetscLayoutCreate(((PetscObject) A)->comm, &rLayout);
696: PetscLayoutSetLocalSize(rLayout, locRows);
697: PetscLayoutSetBlockSize(rLayout, 1);
698: PetscLayoutSetUp(rLayout);
699: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
700: PetscLayoutDestroy(&rLayout);
701: for(r = rStart; r < rEnd; ++r) {
702: PetscInt numCols, cStart, c;
704: PetscSectionGetDof(sectionAdj, r, &numCols);
705: PetscSectionGetOffset(sectionAdj, r, &cStart);
706: for(c = cStart; c < cStart+numCols; ++c) {
707: if ((cols[c] >= rStart) && (cols[c] < rEnd)) {
708: ++dnz[r-rStart];
709: if (cols[c] >= r) {++dnzu[r-rStart];}
710: } else {
711: ++onz[r-rStart];
712: if (cols[c] >= r) {++onzu[r-rStart];}
713: }
714: }
715: }
716: /* Set matrix pattern */
717: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
718: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
719: /* Fill matrix with zeros */
720: if (fillMatrix) {
721: PetscScalar *values;
722: PetscInt maxRowLen = 0;
724: for(r = rStart; r < rEnd; ++r) {
725: PetscInt len;
727: PetscSectionGetDof(sectionAdj, r, &len);
728: maxRowLen = PetscMax(maxRowLen, len);
729: }
730: PetscMalloc(maxRowLen * sizeof(PetscScalar), &values);
731: PetscMemzero(values, maxRowLen * sizeof(PetscScalar));
732: for(r = rStart; r < rEnd; ++r) {
733: PetscInt numCols, cStart;
735: PetscSectionGetDof(sectionAdj, r, &numCols);
736: PetscSectionGetOffset(sectionAdj, r, &cStart);
737: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
738: }
739: PetscFree(values);
740: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
741: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
742: }
743: PetscSectionDestroy(§ionAdj);
744: PetscFree(cols);
745: return(0);
746: }
748: #if 0
751: PetscErrorCode DMComplexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
752: {
754: PetscInt c,cStart,cEnd,pStart,pEnd;
755: PetscInt *tmpClosure,*tmpAdj,*visits;
758: DMComplexGetDimension(dm, &dim);
759: DMComplexGetDepth(dm, &depth);
760: DMComplexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
761: maxClosureSize = 2*PetscMax(pow(mesh->maxConeSize, depth)+1, pow(mesh->maxSupportSize, depth)+1);
762: PetscSectionGetChart(section, &pStart, &pEnd);
763: npoints = pEnd - pStart;
764: PetscMalloc3(maxClosureSize,PetscInt,&tmpClosure,npoints,PetscInt,&lvisits,npoints,PetscInt,&visits);
765: PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));
766: PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));
767: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
768: for (c=cStart; c<cEnd; c++) {
769: PetscInt *support = tmpClosure;
770: DMComplexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
771: for (p=0; p<supportSize; p++) {
772: lvisits[support[p]]++;
773: }
774: }
775: PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
776: PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM);
777: PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);
778: PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits);
780: PetscSFGetRanks();
783: PetscMalloc2(maxClosureSize*maxClosureSize,PetscInt,&cellmat,npoints,PetscInt,&owner);
784: for (c=cStart; c<cEnd; c++) {
785: PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));
786: /*
787: Depth-first walk of transitive closure.
788: At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
789: This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
790: */
791: }
793: PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
794: PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM);
795: return(0);
796: }
797: #endif
801: PetscErrorCode DMCreateMatrix_Complex(DM dm, const MatType mtype, Mat *J)
802: {
803: PetscSection section, sectionGlobal;
804: PetscInt bs = -1;
805: PetscInt localSize;
806: PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
810: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
811: MatInitializePackage(PETSC_NULL);
812: #endif
813: if (!mtype) mtype = MATAIJ;
814: DMGetDefaultSection(dm, §ion);
815: DMGetDefaultGlobalSection(dm, §ionGlobal);
816: /* PetscSectionGetStorageSize(sectionGlobal, &localSize); */
817: PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);
818: MatCreate(((PetscObject) dm)->comm, J);
819: MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);
820: MatSetType(*J, mtype);
821: MatSetFromOptions(*J);
822: PetscStrcmp(mtype, MATSHELL, &isShell);
823: PetscStrcmp(mtype, MATBAIJ, &isBlock);
824: PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);
825: PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);
826: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
827: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
828: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
829: /* Check for symmetric storage */
830: isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
831: if (isSymmetric) {
832: MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
833: }
834: if (!isShell) {
835: PetscBool fillMatrix = (PetscBool) !dm->prealloc_only;
836: PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal;
838: if (bs < 0) {
839: if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
840: PetscInt pStart, pEnd, p, dof;
842: DMComplexGetChart(dm, &pStart, &pEnd);
843: for(p = pStart; p < pEnd; ++p) {
844: PetscSectionGetDof(sectionGlobal, p, &dof);
845: if (dof) {
846: bs = dof;
847: break;
848: }
849: }
850: } else {
851: bs = 1;
852: }
853: /* Must have same blocksize on all procs (some might have no points) */
854: bsLocal = bs;
855: MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, ((PetscObject) dm)->comm);
856: }
857: PetscMalloc4(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz, localSize/bs, PetscInt, &dnzu, localSize/bs, PetscInt, &onzu);
858: PetscMemzero(dnz, localSize/bs * sizeof(PetscInt));
859: PetscMemzero(onz, localSize/bs * sizeof(PetscInt));
860: PetscMemzero(dnzu, localSize/bs * sizeof(PetscInt));
861: PetscMemzero(onzu, localSize/bs * sizeof(PetscInt));
862: DMComplexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);
863: PetscFree4(dnz, onz, dnzu, onzu);
864: }
865: return(0);
866: }
870: /*@
871: DMComplexGetDimension - Return the topological mesh dimension
873: Not collective
875: Input Parameter:
876: . mesh - The DMComplex
878: Output Parameter:
879: . dim - The topological mesh dimension
881: Level: beginner
883: .seealso: DMComplexCreate()
884: @*/
885: PetscErrorCode DMComplexGetDimension(DM dm, PetscInt *dim)
886: {
887: DM_Complex *mesh = (DM_Complex *) dm->data;
892: *dim = mesh->dim;
893: return(0);
894: }
898: /*@
899: DMComplexSetDimension - Set the topological mesh dimension
901: Collective on mesh
903: Input Parameters:
904: + mesh - The DMComplex
905: - dim - The topological mesh dimension
907: Level: beginner
909: .seealso: DMComplexCreate()
910: @*/
911: PetscErrorCode DMComplexSetDimension(DM dm, PetscInt dim)
912: {
913: DM_Complex *mesh = (DM_Complex *) dm->data;
918: mesh->dim = dim;
919: return(0);
920: }
924: /*@
925: DMComplexGetChart - Return the interval for all mesh points [pStart, pEnd)
927: Not collective
929: Input Parameter:
930: . mesh - The DMComplex
932: Output Parameters:
933: + pStart - The first mesh point
934: - pEnd - The upper bound for mesh points
936: Level: beginner
938: .seealso: DMComplexCreate(), DMComplexSetChart()
939: @*/
940: PetscErrorCode DMComplexGetChart(DM dm, PetscInt *pStart, PetscInt *pEnd)
941: {
942: DM_Complex *mesh = (DM_Complex *) dm->data;
947: PetscSectionGetChart(mesh->coneSection, pStart, pEnd);
948: return(0);
949: }
953: /*@
954: DMComplexSetChart - Set the interval for all mesh points [pStart, pEnd)
956: Not collective
958: Input Parameters:
959: + mesh - The DMComplex
960: . pStart - The first mesh point
961: - pEnd - The upper bound for mesh points
963: Output Parameters:
965: Level: beginner
967: .seealso: DMComplexCreate(), DMComplexGetChart()
968: @*/
969: PetscErrorCode DMComplexSetChart(DM dm, PetscInt pStart, PetscInt pEnd)
970: {
971: DM_Complex *mesh = (DM_Complex *) dm->data;
976: PetscSectionSetChart(mesh->coneSection, pStart, pEnd);
977: return(0);
978: }
982: /*@
983: DMComplexGetConeSize - Return the number of in-edges for this point in the Sieve DAG
985: Not collective
987: Input Parameters:
988: + mesh - The DMComplex
989: - p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
991: Output Parameter:
992: . size - The cone size for point p
994: Level: beginner
996: .seealso: DMComplexCreate(), DMComplexSetConeSize(), DMComplexSetChart()
997: @*/
998: PetscErrorCode DMComplexGetConeSize(DM dm, PetscInt p, PetscInt *size)
999: {
1000: DM_Complex *mesh = (DM_Complex *) dm->data;
1006: PetscSectionGetDof(mesh->coneSection, p, size);
1007: return(0);
1008: }
1012: /*@
1013: DMComplexSetConeSize - Set the number of in-edges for this point in the Sieve DAG
1015: Not collective
1017: Input Parameters:
1018: + mesh - The DMComplex
1019: . p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1020: - size - The cone size for point p
1022: Output Parameter:
1024: Note:
1025: This should be called after DMComplexSetChart().
1027: Level: beginner
1029: .seealso: DMComplexCreate(), DMComplexGetConeSize(), DMComplexSetChart()
1030: @*/
1031: PetscErrorCode DMComplexSetConeSize(DM dm, PetscInt p, PetscInt size)
1032: {
1033: DM_Complex *mesh = (DM_Complex *) dm->data;
1038: PetscSectionSetDof(mesh->coneSection, p, size);
1039: mesh->maxConeSize = PetscMax(mesh->maxConeSize, size);
1040: return(0);
1041: }
1045: /*@C
1046: DMComplexGetCone - Return the points on the in-edges for this point in the Sieve DAG
1048: Not collective
1050: Input Parameters:
1051: + mesh - The DMComplex
1052: - p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1054: Output Parameter:
1055: . cone - An array of points which are on the in-edges for point p
1057: Level: beginner
1059: Note:
1060: This routine is not available in Fortran.
1062: .seealso: DMComplexCreate(), DMComplexSetCone(), DMComplexSetChart()
1063: @*/
1064: PetscErrorCode DMComplexGetCone(DM dm, PetscInt p, const PetscInt *cone[])
1065: {
1066: DM_Complex *mesh = (DM_Complex *) dm->data;
1067: PetscInt off;
1073: PetscSectionGetOffset(mesh->coneSection, p, &off);
1074: *cone = &mesh->cones[off];
1075: return(0);
1076: }
1080: /*@
1081: DMComplexSetCone - Set the points on the in-edges for this point in the Sieve DAG
1083: Not collective
1085: Input Parameters:
1086: + mesh - The DMComplex
1087: . p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1088: - cone - An array of points which are on the in-edges for point p
1090: Output Parameter:
1092: Note:
1093: This should be called after all calls to DMComplexSetConeSize() and DMSetUp().
1095: Level: beginner
1097: .seealso: DMComplexCreate(), DMComplexGetCone(), DMComplexSetChart(), DMComplexSetConeSize(), DMSetUp()
1098: @*/
1099: PetscErrorCode DMComplexSetCone(DM dm, PetscInt p, const PetscInt cone[])
1100: {
1101: DM_Complex *mesh = (DM_Complex *) dm->data;
1102: PetscInt pStart, pEnd;
1103: PetscInt dof, off, c;
1109: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1110: PetscSectionGetDof(mesh->coneSection, p, &dof);
1111: PetscSectionGetOffset(mesh->coneSection, p, &off);
1112: if ((p < pStart) || (p >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1113: for(c = 0; c < dof; ++c) {
1114: if ((cone[c] < pStart) || (cone[c] >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D. %D)", cone[c], pStart, pEnd);
1115: mesh->cones[off+c] = cone[c];
1116: }
1117: return(0);
1118: }
1122: /*@C
1123: DMComplexGetConeOrientation - Return the orientations on the in-edges for this point in the Sieve DAG
1125: Not collective
1127: Input Parameters:
1128: + mesh - The DMComplex
1129: - p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1131: Output Parameter:
1132: . coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1133: integer giving the prescription for cone traversal. If it is negative, the cone is
1134: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1135: the index of the cone point on which to start.
1137: Level: beginner
1139: Note:
1140: This routine is not available in Fortran.
1142: .seealso: DMComplexCreate(), DMComplexGetCone(), DMComplexSetCone(), DMComplexSetChart()
1143: @*/
1144: PetscErrorCode DMComplexGetConeOrientation(DM dm, PetscInt p, const PetscInt *coneOrientation[])
1145: {
1146: DM_Complex *mesh = (DM_Complex *) dm->data;
1147: PetscInt off;
1153: PetscSectionGetOffset(mesh->coneSection, p, &off);
1154: *coneOrientation = &mesh->coneOrientations[off];
1155: return(0);
1156: }
1160: /*@
1161: DMComplexSetConeOrientation - Set the orientations on the in-edges for this point in the Sieve DAG
1163: Not collective
1165: Input Parameters:
1166: + mesh - The DMComplex
1167: . p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1168: - coneOrientation - An array of orientations which are on the in-edges for point p. An orientation is an
1169: integer giving the prescription for cone traversal. If it is negative, the cone is
1170: traversed in the opposite direction. Its value 'o', or if negative '-(o+1)', gives
1171: the index of the cone point on which to start.
1173: Output Parameter:
1175: Note:
1176: This should be called after all calls to DMComplexSetConeSize() and DMSetUp().
1178: Level: beginner
1180: .seealso: DMComplexCreate(), DMComplexGetConeOrientation(), DMComplexSetCone(), DMComplexSetChart(), DMComplexSetConeSize(), DMSetUp()
1181: @*/
1182: PetscErrorCode DMComplexSetConeOrientation(DM dm, PetscInt p, const PetscInt coneOrientation[])
1183: {
1184: DM_Complex *mesh = (DM_Complex *) dm->data;
1185: PetscInt pStart, pEnd;
1186: PetscInt dof, off, c;
1192: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1193: PetscSectionGetDof(mesh->coneSection, p, &dof);
1194: PetscSectionGetOffset(mesh->coneSection, p, &off);
1195: if ((p < pStart) || (p >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1196: for(c = 0; c < dof; ++c) {
1197: PetscInt cdof, o = coneOrientation[c];
1199: PetscSectionGetDof(mesh->coneSection, mesh->cones[off+c], &cdof);
1200: if ((o < -(cdof+1)) || (o >= cdof)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone orientation %D is not in the valid range [%D. %D)", o, -(cdof+1), cdof);
1201: mesh->coneOrientations[off+c] = o;
1202: }
1203: return(0);
1204: }
1208: PetscErrorCode DMComplexInsertCone(DM dm, PetscInt p, PetscInt conePos, PetscInt conePoint)
1209: {
1210: DM_Complex *mesh = (DM_Complex *) dm->data;
1211: PetscInt pStart, pEnd;
1212: PetscInt dof, off;
1217: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1218: PetscSectionGetDof(mesh->coneSection, p, &dof);
1219: PetscSectionGetOffset(mesh->coneSection, p, &off);
1220: if ((p < pStart) || (p >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Mesh point %D is not in the valid range [%D, %D)", p, pStart, pEnd);
1221: if ((conePoint < pStart) || (conePoint >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone point %D is not in the valid range [%D, %D)", conePoint, pStart, pEnd);
1222: if (conePos >= dof) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone position %D is not in the valid range [0, %D)", conePos, dof);
1223: mesh->cones[off+conePos] = conePoint;
1224: return(0);
1225: }
1229: /*@
1230: DMComplexGetSupportSize - Return the number of out-edges for this point in the Sieve DAG
1232: Not collective
1234: Input Parameters:
1235: + mesh - The DMComplex
1236: - p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1238: Output Parameter:
1239: . size - The support size for point p
1241: Level: beginner
1243: .seealso: DMComplexCreate(), DMComplexSetConeSize(), DMComplexSetChart(), DMComplexGetConeSize()
1244: @*/
1245: PetscErrorCode DMComplexGetSupportSize(DM dm, PetscInt p, PetscInt *size)
1246: {
1247: DM_Complex *mesh = (DM_Complex *) dm->data;
1253: PetscSectionGetDof(mesh->supportSection, p, size);
1254: return(0);
1255: }
1259: /*@C
1260: DMComplexGetSupport - Return the points on the out-edges for this point in the Sieve DAG
1262: Not collective
1264: Input Parameters:
1265: + mesh - The DMComplex
1266: - p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1268: Output Parameter:
1269: . support - An array of points which are on the out-edges for point p
1271: Level: beginner
1273: .seealso: DMComplexCreate(), DMComplexSetCone(), DMComplexSetChart(), DMComplexGetCone()
1274: @*/
1275: PetscErrorCode DMComplexGetSupport(DM dm, PetscInt p, const PetscInt *support[])
1276: {
1277: DM_Complex *mesh = (DM_Complex *) dm->data;
1278: PetscInt off;
1284: PetscSectionGetOffset(mesh->supportSection, p, &off);
1285: *support = &mesh->supports[off];
1286: return(0);
1287: }
1291: /*@C
1292: DMComplexGetTransitiveClosure - Return the points on the transitive closure of the in-edges or out-edges for this point in the Sieve DAG
1294: Not collective
1296: Input Parameters:
1297: + mesh - The DMComplex
1298: . p - The Sieve point, which must lie in the chart set with DMComplexSetChart()
1299: . useCone - PETSC_TRUE for in-edges, otherwise use out-edges
1300: - points - If points is PETSC_NULL on input, internal storage will be returned, otherwise the provided array is used
1302: Output Parameters:
1303: + numPoints - The number of points in the closure, so points[] is of size 2*numPoints
1304: - points - The points and point orientations, interleaved as pairs [p0, o0, p1, o1, ...]
1306: Note:
1307: If using internal storage (points is PETSC_NULL on input), each call overwrites the last output.
1309: Level: beginner
1311: .seealso: DMComplexCreate(), DMComplexSetCone(), DMComplexSetChart(), DMComplexGetCone()
1312: @*/
1313: PetscErrorCode DMComplexGetTransitiveClosure(DM dm, PetscInt p, PetscBool useCone, PetscInt *numPoints, PetscInt *points[])
1314: {
1315: DM_Complex *mesh = (DM_Complex *) dm->data;
1316: PetscInt *closure, *fifo;
1317: const PetscInt *tmp, *tmpO = PETSC_NULL;
1318: PetscInt tmpSize, t;
1319: PetscInt closureSize = 2, fifoSize = 0, fifoStart = 0;
1320: PetscErrorCode ierr;
1324: if (!mesh->closureTmpA) {
1325: PetscInt depth, maxSize;
1327: DMComplexGetDepth(dm, &depth);
1328: maxSize = (PetscInt) (2*PetscMax(pow((PetscReal) mesh->maxConeSize, depth)+1, pow((PetscReal) mesh->maxSupportSize, depth)+1));
1329: PetscMalloc2(maxSize,PetscInt,&mesh->closureTmpA,maxSize,PetscInt,&mesh->closureTmpB);
1330: }
1331: closure = *points ? *points : mesh->closureTmpA;
1332: fifo = mesh->closureTmpB;
1333: closure[0] = p; closure[1] = 0;
1334: /* This is only 1-level */
1335: if (useCone) {
1336: DMComplexGetConeSize(dm, p, &tmpSize);
1337: DMComplexGetCone(dm, p, &tmp);
1338: DMComplexGetConeOrientation(dm, p, &tmpO);
1339: } else {
1340: DMComplexGetSupportSize(dm, p, &tmpSize);
1341: DMComplexGetSupport(dm, p, &tmp);
1342: }
1343: for(t = 0; t < tmpSize; ++t, closureSize += 2, fifoSize += 2) {
1344: const PetscInt cp = tmp[t];
1345: const PetscInt co = tmpO ? tmpO[t] : 0;
1347: closure[closureSize] = cp;
1348: closure[closureSize+1] = co;
1349: fifo[fifoSize] = cp;
1350: fifo[fifoSize+1] = co;
1351: }
1352: while(fifoSize - fifoStart) {
1353: const PetscInt q = fifo[fifoStart];
1354: const PetscInt o = fifo[fifoStart+1];
1355: const PetscInt rev = o >= 0 ? 0 : 1;
1356: const PetscInt off = rev ? -(o+1) : o;
1358: if (useCone) {
1359: DMComplexGetConeSize(dm, q, &tmpSize);
1360: DMComplexGetCone(dm, q, &tmp);
1361: DMComplexGetConeOrientation(dm, q, &tmpO);
1362: } else {
1363: DMComplexGetSupportSize(dm, q, &tmpSize);
1364: DMComplexGetSupport(dm, q, &tmp);
1365: tmpO = PETSC_NULL;
1366: }
1367: for(t = 0; t < tmpSize; ++t) {
1368: const PetscInt i = ((rev ? tmpSize-t : t) + off)%tmpSize;
1369: const PetscInt cp = tmp[i];
1370: const PetscInt co = tmpO ? tmpO[i] : 0;
1371: PetscInt c;
1373: /* Check for duplicate */
1374: for(c = 0; c < closureSize; c += 2) {
1375: if (closure[c] == cp) break;
1376: }
1377: if (c == closureSize) {
1378: closure[closureSize] = cp;
1379: closure[closureSize+1] = co;
1380: fifo[fifoSize] = cp;
1381: fifo[fifoSize+1] = co;
1382: closureSize += 2;
1383: fifoSize += 2;
1384: }
1385: }
1386: fifoStart += 2;
1387: }
1388: if (numPoints) *numPoints = closureSize/2;
1389: if (points) *points = closure;
1390: return(0);
1391: }
1395: /*
1396: DMComplexGetFaces -
1398: Note: This will only work for cell-vertex meshes.
1399: */
1400: PetscErrorCode DMComplexGetFaces(DM dm, PetscInt p, PetscInt *numFaces, PetscInt *faceSize, const PetscInt *faces[])
1401: {
1402: DM_Complex *mesh = (DM_Complex *) dm->data;
1403: const PetscInt *cone;
1404: PetscInt depth, dim, coneSize;
1405: PetscErrorCode ierr;
1409: DMComplexGetDimension(dm, &dim);
1410: DMComplexGetDepth(dm, &depth);
1411: if (depth > 1) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Faces can only be returned for cell-vertex meshes.");
1412: if (!mesh->facesTmp) {PetscMalloc(mesh->maxConeSize*mesh->maxSupportSize * sizeof(PetscInt), &mesh->facesTmp);}
1413: DMComplexGetConeSize(dm, p, &coneSize);
1414: DMComplexGetCone(dm, p, &cone);
1415: switch(dim) {
1416: case 2:
1417: switch(coneSize) {
1418: case 3:
1419: mesh->facesTmp[0] = cone[0]; mesh->facesTmp[1] = cone[1];
1420: mesh->facesTmp[2] = cone[1]; mesh->facesTmp[3] = cone[2];
1421: mesh->facesTmp[4] = cone[2]; mesh->facesTmp[5] = cone[0];
1422: *numFaces = 3;
1423: *faceSize = 2;
1424: *faces = mesh->facesTmp;
1425: break;
1426: default:
1427: SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Cone size %D not supported for dimension %", coneSize, dim);
1428: }
1429: break;
1430: default:
1431: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension % not supported", dim);
1432: }
1433: return(0);
1434: }
1438: /*@
1439: DMComplexGetMaxSizes - Return the maximum number of in-edges (cone) and out-edges (support) for any point in the Sieve DAG
1441: Not collective
1443: Input Parameter:
1444: . mesh - The DMComplex
1446: Output Parameters:
1447: + maxConeSize - The maximum number of in-edges
1448: - maxSupportSize - The maximum number of out-edges
1450: Level: beginner
1452: .seealso: DMComplexCreate(), DMComplexSetConeSize(), DMComplexSetChart()
1453: @*/
1454: PetscErrorCode DMComplexGetMaxSizes(DM dm, PetscInt *maxConeSize, PetscInt *maxSupportSize)
1455: {
1456: DM_Complex *mesh = (DM_Complex *) dm->data;
1460: if (maxConeSize) *maxConeSize = mesh->maxConeSize;
1461: if (maxSupportSize) *maxSupportSize = mesh->maxSupportSize;
1462: return(0);
1463: }
1467: PetscErrorCode DMSetUp_Complex(DM dm)
1468: {
1469: DM_Complex *mesh = (DM_Complex *) dm->data;
1470: PetscInt size;
1475: PetscSectionSetUp(mesh->coneSection);
1476: PetscSectionGetStorageSize(mesh->coneSection, &size);
1477: PetscMalloc(size * sizeof(PetscInt), &mesh->cones);
1478: PetscMalloc(size * sizeof(PetscInt), &mesh->coneOrientations);
1479: PetscMemzero(mesh->coneOrientations, size * sizeof(PetscInt));
1480: return(0);
1481: }
1485: /*@
1486: DMComplexSymmetrize - Creates support (out-edge) information from cone (in-edge) inoformation
1488: Not collective
1490: Input Parameter:
1491: . mesh - The DMComplex
1493: Output Parameter:
1495: Note:
1496: This should be called after all calls to DMComplexSetCone()
1498: Level: beginner
1500: .seealso: DMComplexCreate(), DMComplexSetChart(), DMComplexSetConeSize(), DMComplexSetCone()
1501: @*/
1502: PetscErrorCode DMComplexSymmetrize(DM dm)
1503: {
1504: DM_Complex *mesh = (DM_Complex *) dm->data;
1505: PetscInt *offsets;
1506: PetscInt supportSize;
1507: PetscInt pStart, pEnd, p;
1512: /* Calculate support sizes */
1513: DMComplexGetChart(dm, &pStart, &pEnd);
1514: PetscSectionSetChart(mesh->supportSection, pStart, pEnd);
1515: for(p = pStart; p < pEnd; ++p) {
1516: PetscInt dof, off, c;
1518: PetscSectionGetDof(mesh->coneSection, p, &dof);
1519: PetscSectionGetOffset(mesh->coneSection, p, &off);
1520: for(c = off; c < off+dof; ++c) {
1521: PetscSectionAddDof(mesh->supportSection, mesh->cones[c], 1);
1522: }
1523: }
1524: for(p = pStart; p < pEnd; ++p) {
1525: PetscInt dof;
1527: PetscSectionGetDof(mesh->supportSection, p, &dof);
1528: mesh->maxSupportSize = PetscMax(mesh->maxSupportSize, dof);
1529: }
1530: PetscSectionSetUp(mesh->supportSection);
1531: /* Calculate supports */
1532: PetscSectionGetStorageSize(mesh->supportSection, &supportSize);
1533: PetscMalloc(supportSize * sizeof(PetscInt), &mesh->supports);
1534: PetscMalloc((pEnd - pStart) * sizeof(PetscInt), &offsets);
1535: PetscMemzero(offsets, (pEnd - pStart) * sizeof(PetscInt));
1536: for(p = pStart; p < pEnd; ++p) {
1537: PetscInt dof, off, c;
1539: PetscSectionGetDof(mesh->coneSection, p, &dof);
1540: PetscSectionGetOffset(mesh->coneSection, p, &off);
1541: for(c = off; c < off+dof; ++c) {
1542: const PetscInt q = mesh->cones[c];
1543: PetscInt offS;
1545: PetscSectionGetOffset(mesh->supportSection, q, &offS);
1546: mesh->supports[offS+offsets[q]] = p;
1547: ++offsets[q];
1548: }
1549: }
1550: PetscFree(offsets);
1551: return(0);
1552: }
1556: PetscErrorCode DMComplexSetDepth_Private(DM dm, PetscInt p, PetscInt *depth)
1557: {
1558: PetscInt d;
1562: DMComplexGetLabelValue(dm, "depth", p, &d);
1563: if (d < 0) {
1564: /* We are guaranteed that the point has a cone since the depth was not yet set */
1565: const PetscInt *cone;
1566: PetscInt dCone;
1568: DMComplexGetCone(dm, p, &cone);
1569: DMComplexSetDepth_Private(dm, cone[0], &dCone);
1570: d = dCone+1;
1571: DMComplexSetLabelValue(dm, "depth", p, d);
1572: }
1573: *depth = d;
1574: return(0);
1575: }
1579: /*@
1580: DMComplexStratify - The Sieve DAG for most topologies is a graded poset (http://en.wikipedia.org/wiki/Graded_poset), and
1581: can be illustrated by Hasse Diagram (a http://en.wikipedia.org/wiki/Hasse_diagram). The strata group all points of the
1582: same grade, and this function calculates the strata. This grade can be seen as the height (or depth) of the point in
1583: the DAG.
1585: Not collective
1587: Input Parameter:
1588: . mesh - The DMComplex
1590: Output Parameter:
1592: Notes:
1593: The normal association for the point grade is element dimension (or co-dimension). For instance, all vertices would
1594: have depth 0, and all edges depth 1. Likewise, all cells heights would have height 0, and all faces height 1.
1596: This should be called after all calls to DMComplexSymmetrize()
1598: Level: beginner
1600: .seealso: DMComplexCreate(), DMComplexSymmetrize()
1601: @*/
1602: PetscErrorCode DMComplexStratify(DM dm)
1603: {
1604: DM_Complex *mesh = (DM_Complex *) dm->data;
1605: PetscInt pStart, pEnd, p;
1606: PetscInt numRoots = 0, numLeaves = 0;
1611: /* Calculate depth */
1612: PetscSectionGetChart(mesh->coneSection, &pStart, &pEnd);
1613: /* Initialize roots and count leaves */
1614: for(p = pStart; p < pEnd; ++p) {
1615: PetscInt coneSize, supportSize;
1617: DMComplexGetConeSize(dm, p, &coneSize);
1618: DMComplexGetSupportSize(dm, p, &supportSize);
1619: if (!coneSize && supportSize) {
1620: ++numRoots;
1621: DMComplexSetLabelValue(dm, "depth", p, 0);
1622: } else if (!supportSize && coneSize) {
1623: ++numLeaves;
1624: }
1625: }
1626: if (numRoots + numLeaves == (pEnd - pStart)) {
1627: for(p = pStart; p < pEnd; ++p) {
1628: PetscInt coneSize, supportSize;
1630: DMComplexGetConeSize(dm, p, &coneSize);
1631: DMComplexGetSupportSize(dm, p, &supportSize);
1632: if (!supportSize && coneSize) {
1633: DMComplexSetLabelValue(dm, "depth", p, 1);
1634: }
1635: }
1636: } else {
1637: /* This might be slow since lookup is not fast */
1638: for(p = pStart; p < pEnd; ++p) {
1639: PetscInt depth;
1641: DMComplexSetDepth_Private(dm, p, &depth);
1642: }
1643: }
1644: return(0);
1645: }
1649: /*@C
1650: DMComplexHasLabel - Determine whether the mesh has a label of a given name
1652: Not Collective
1654: Input Parameters:
1655: + dm - The DMComplex object
1656: - name - The label name
1658: Output Parameter:
1659: . hasLabel - PETSC_TRUE if the label is present
1661: Level: intermediate
1663: .keywords: mesh
1664: .seealso: DMComplexGetLabelValue(), DMComplexSetLabelValue(), DMComplexGetLabelStratum()
1665: @*/
1666: PetscErrorCode DMComplexHasLabel(DM dm, const char name[], PetscBool *hasLabel)
1667: {
1668: DM_Complex *mesh = (DM_Complex *) dm->data;
1669: SieveLabel next = mesh->labels;
1676: *hasLabel = PETSC_FALSE;
1677: while(next) {
1678: PetscStrcmp(name, next->name, hasLabel);
1679: if (*hasLabel) break;
1680: next = next->next;
1681: }
1682: return(0);
1683: }
1687: /*@C
1688: DMComplexGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default
1690: Not Collective
1692: Input Parameters:
1693: + dm - The DMComplex object
1694: . name - The label name
1695: - point - The mesh point
1697: Output Parameter:
1698: . value - The label value for this point, or -1 if the point is not in the label
1700: Level: beginner
1702: .keywords: mesh
1703: .seealso: DMComplexSetLabelValue(), DMComplexGetLabelStratum()
1704: @*/
1705: PetscErrorCode DMComplexGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
1706: {
1707: DM_Complex *mesh = (DM_Complex *) dm->data;
1708: SieveLabel next = mesh->labels;
1709: PetscBool flg;
1710: PetscInt v, p;
1716: *value = -1;
1717: DMComplexHasLabel(dm, name, &flg);
1718: if (!flg) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
1719: /* We should have a generic GetLabel() and a Label class */
1720: while(next) {
1721: PetscStrcmp(name, next->name, &flg);
1722: if (flg) break;
1723: next = next->next;
1724: }
1725: /* Find, or add, label value */
1726: for(v = 0; v < next->numStrata; ++v) {
1727: for(p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1728: if (next->points[p] == point) {
1729: *value = next->stratumValues[v];
1730: break;
1731: }
1732: }
1733: }
1734: return(0);
1735: }
1739: /*@C
1740: DMComplexSetLabelValue - Add a point to a Sieve Label with given value
1742: Not Collective
1744: Input Parameters:
1745: + dm - The DMComplex object
1746: . name - The label name
1747: . point - The mesh point
1748: - value - The label value for this point
1750: Output Parameter:
1752: Level: beginner
1754: .keywords: mesh
1755: .seealso: DMComplexGetLabelStratum()
1756: @*/
1757: PetscErrorCode DMComplexSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
1758: {
1759: DM_Complex *mesh = (DM_Complex *) dm->data;
1760: SieveLabel next = mesh->labels;
1761: PetscBool flg = PETSC_FALSE;
1762: PetscInt v, p;
1768: /* Find, or create, label */
1769: while(next) {
1770: PetscStrcmp(name, next->name, &flg);
1771: if (flg) break;
1772: next = next->next;
1773: }
1774: if (!flg) {
1775: SieveLabel tmpLabel = mesh->labels;
1776: PetscNew(struct Sieve_Label, &mesh->labels);
1777: mesh->labels->next = tmpLabel;
1778: next = mesh->labels;
1779: PetscStrallocpy(name, &next->name);
1780: }
1781: /* Find, or add, label value */
1782: for(v = 0; v < next->numStrata; ++v) {
1783: if (next->stratumValues[v] == value) break;
1784: }
1785: if (v >= next->numStrata) {
1786: PetscInt *tmpV, *tmpO, *tmpS;
1787: PetscMalloc3(next->numStrata+1,PetscInt,&tmpV,next->numStrata+2,PetscInt,&tmpO,next->numStrata+1,PetscInt,&tmpS);
1788: for(v = 0; v < next->numStrata; ++v) {
1789: tmpV[v] = next->stratumValues[v];
1790: tmpO[v] = next->stratumOffsets[v];
1791: tmpS[v] = next->stratumSizes[v];
1792: }
1793: tmpV[v] = value;
1794: tmpO[v] = v == 0 ? 0 : next->stratumOffsets[v];
1795: tmpS[v] = 0;
1796: tmpO[v+1] = tmpO[v];
1797: ++next->numStrata;
1798: PetscFree3(next->stratumValues,next->stratumOffsets,next->stratumSizes);
1799: next->stratumValues = tmpV;
1800: next->stratumOffsets = tmpO;
1801: next->stratumSizes = tmpS;
1802: }
1803: /* Check whether point exists */
1804: for(p = next->stratumOffsets[v]; p < next->stratumOffsets[v]+next->stratumSizes[v]; ++p) {
1805: if (next->points[p] == point) {
1806: break;
1807: }
1808: }
1809: /* Add point: NEED TO OPTIMIZE */
1810: if (p >= next->stratumOffsets[v]+next->stratumSizes[v]) {
1811: /* Check for reallocation */
1812: if (next->stratumSizes[v] >= next->stratumOffsets[v+1]-next->stratumOffsets[v]) {
1813: PetscInt oldSize = next->stratumOffsets[v+1]-next->stratumOffsets[v];
1814: PetscInt newSize = PetscMax(10, 2*oldSize); /* Double the size, since 2 is the optimal base for this online algorithm */
1815: PetscInt shift = newSize - oldSize;
1816: PetscInt allocSize = next->stratumOffsets[next->numStrata] + shift;
1817: PetscInt *newPoints;
1818: PetscInt w, q;
1820: PetscMalloc(allocSize * sizeof(PetscInt), &newPoints);
1821: for(q = 0; q < next->stratumOffsets[v]+next->stratumSizes[v]; ++q) {
1822: newPoints[q] = next->points[q];
1823: }
1824: for(w = v+1; w < next->numStrata; ++w) {
1825: for(q = next->stratumOffsets[w]; q < next->stratumOffsets[w]+next->stratumSizes[w]; ++q) {
1826: newPoints[q+shift] = next->points[q];
1827: }
1828: next->stratumOffsets[w] += shift;
1829: }
1830: next->stratumOffsets[next->numStrata] += shift;
1831: PetscFree(next->points);
1832: next->points = newPoints;
1833: }
1834: /* Insert point and resort */
1835: next->points[next->stratumOffsets[v]+next->stratumSizes[v]] = point;
1836: ++next->stratumSizes[v];
1837: PetscSortInt(next->stratumSizes[v], &next->points[next->stratumOffsets[v]]);
1838: }
1839: return(0);
1840: }
1844: /*@C
1845: DMComplexGetLabelSize - Get the number of different integer ids in a Label
1847: Not Collective
1849: Input Parameters:
1850: + dm - The DMComplex object
1851: - name - The label name
1853: Output Parameter:
1854: . size - The label size (number of different integer ids)
1856: Level: beginner
1858: .keywords: mesh
1859: .seealso: DMComplexSetLabelValue()
1860: @*/
1861: PetscErrorCode DMComplexGetLabelSize(DM dm, const char name[], PetscInt *size)
1862: {
1863: DM_Complex *mesh = (DM_Complex *) dm->data;
1864: SieveLabel next = mesh->labels;
1865: PetscBool flg;
1872: *size = 0;
1873: while(next) {
1874: PetscStrcmp(name, next->name, &flg);
1875: if (flg) {
1876: *size = next->numStrata;
1877: break;
1878: }
1879: next = next->next;
1880: }
1881: return(0);
1882: }
1886: /*@C
1887: DMComplexGetLabelIdIS - Get the integer ids in a label
1889: Not Collective
1891: Input Parameters:
1892: + mesh - The DMComplex object
1893: - name - The label name
1895: Output Parameter:
1896: . ids - The integer ids
1898: Level: beginner
1900: .keywords: mesh
1901: .seealso: DMComplexGetLabelSize()
1902: @*/
1903: PetscErrorCode DMComplexGetLabelIdIS(DM dm, const char name[], IS *ids)
1904: {
1905: DM_Complex *mesh = (DM_Complex *) dm->data;
1906: SieveLabel next = mesh->labels;
1907: PetscInt *values;
1908: PetscInt size=-1, i = 0;
1909: PetscBool flg;
1916: while(next) {
1917: PetscStrcmp(name, next->name, &flg);
1918: if (flg) {
1919: size = next->numStrata;
1920: PetscMalloc(size * sizeof(PetscInt), &values);
1921: for(i = 0; i < next->numStrata; ++i) {
1922: values[i] = next->stratumValues[i];
1923: }
1924: break;
1925: }
1926: next = next->next;
1927: }
1928: if (!next) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label with name %s exists in this mesh", name);
1929: ISCreateGeneral(((PetscObject) dm)->comm, size, values, PETSC_OWN_POINTER, ids);
1930: return(0);
1931: }
1935: /*@C
1936: DMComplexGetStratumSize - Get the number of points in a label stratum
1938: Not Collective
1940: Input Parameters:
1941: + dm - The DMComplex object
1942: . name - The label name
1943: - value - The stratum value
1945: Output Parameter:
1946: . size - The stratum size
1948: Level: beginner
1950: .keywords: mesh
1951: .seealso: DMComplexGetLabelSize(), DMComplexGetLabelIds()
1952: @*/
1953: PetscErrorCode DMComplexGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
1954: {
1955: DM_Complex *mesh = (DM_Complex *) dm->data;
1956: SieveLabel next = mesh->labels;
1957: PetscBool flg;
1964: *size = 0;
1965: while(next) {
1966: PetscStrcmp(name, next->name, &flg);
1967: if (flg) {
1968: PetscInt v;
1970: for(v = 0; v < next->numStrata; ++v) {
1971: if (next->stratumValues[v] == value) {
1972: *size = next->stratumSizes[v];
1973: break;
1974: }
1975: }
1976: break;
1977: }
1978: next = next->next;
1979: }
1980: return(0);
1981: }
1985: /*@C
1986: DMComplexGetStratumIS - Get the points in a label stratum
1988: Not Collective
1990: Input Parameters:
1991: + dm - The DMComplex object
1992: . name - The label name
1993: - value - The stratum value
1995: Output Parameter:
1996: . is - The stratum points
1998: Level: beginner
2000: .keywords: mesh
2001: .seealso: DMComplexGetStratumSize()
2002: @*/
2003: PetscErrorCode DMComplexGetStratumIS(DM dm, const char name[], PetscInt value, IS *is) {
2004: DM_Complex *mesh = (DM_Complex *) dm->data;
2005: SieveLabel next = mesh->labels;
2006: PetscBool flg;
2013: *is = PETSC_NULL;
2014: while(next) {
2015: PetscStrcmp(name, next->name, &flg);
2016: if (flg) {
2017: PetscInt v;
2019: for(v = 0; v < next->numStrata; ++v) {
2020: if (next->stratumValues[v] == value) {
2021: ISCreateGeneral(PETSC_COMM_SELF, next->stratumSizes[v], &next->points[next->stratumOffsets[v]], PETSC_COPY_VALUES, is);
2022: break;
2023: }
2024: }
2025: break;
2026: }
2027: next = next->next;
2028: }
2029: return(0);
2030: }
2034: /* This is a 1-level join */
2035: PetscErrorCode DMComplexJoinPoints(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveredPoints, const PetscInt **coveredPoints)
2036: {
2037: DM_Complex *mesh = (DM_Complex *) dm->data;
2038: PetscInt *join[2];
2039: PetscInt joinSize, i = 0;
2040: PetscInt dof, off, p, c, m;
2048: if (!mesh->joinTmpA) {PetscMalloc2(mesh->maxSupportSize,PetscInt,&mesh->joinTmpA,mesh->maxSupportSize,PetscInt,&mesh->joinTmpB);}
2049: join[0] = mesh->joinTmpA; join[1] = mesh->joinTmpB;
2050: /* Copy in support of first point */
2051: PetscSectionGetDof(mesh->supportSection, points[0], &dof);
2052: PetscSectionGetOffset(mesh->supportSection, points[0], &off);
2053: for(joinSize = 0; joinSize < dof; ++joinSize) {
2054: join[i][joinSize] = mesh->supports[off+joinSize];
2055: }
2056: /* Check each successive cone */
2057: for(p = 1; p < numPoints; ++p) {
2058: PetscInt newJoinSize = 0;
2060: PetscSectionGetDof(mesh->supportSection, points[p], &dof);
2061: PetscSectionGetOffset(mesh->supportSection, points[p], &off);
2062: for(c = 0; c < dof; ++c) {
2063: const PetscInt point = mesh->supports[off+c];
2065: for(m = 0; m < joinSize; ++m) {
2066: if (point == join[i][m]) {
2067: join[1-i][newJoinSize++] = point;
2068: break;
2069: }
2070: }
2071: }
2072: joinSize = newJoinSize;
2073: i = 1-i;
2074: }
2075: *numCoveredPoints = joinSize;
2076: *coveredPoints = join[i];
2077: return(0);
2078: }
2082: /* This is a 1-level meet */
2083: PetscErrorCode DMComplexMeetPoints(DM dm, PetscInt numPoints, const PetscInt points[], PetscInt *numCoveringPoints, const PetscInt **coveringPoints)
2084: {
2085: DM_Complex *mesh = (DM_Complex *) dm->data;
2086: PetscInt *meet[2];
2087: PetscInt meetSize, i = 0;
2088: PetscInt dof, off, p, c, m;
2096: if (!mesh->meetTmpA) {PetscMalloc2(mesh->maxConeSize,PetscInt,&mesh->meetTmpA,mesh->maxConeSize,PetscInt,&mesh->meetTmpB);}
2097: meet[0] = mesh->meetTmpA; meet[1] = mesh->meetTmpB;
2098: /* Copy in cone of first point */
2099: PetscSectionGetDof(mesh->coneSection, points[0], &dof);
2100: PetscSectionGetOffset(mesh->coneSection, points[0], &off);
2101: for(meetSize = 0; meetSize < dof; ++meetSize) {
2102: meet[i][meetSize] = mesh->cones[off+meetSize];
2103: }
2104: /* Check each successive cone */
2105: for(p = 1; p < numPoints; ++p) {
2106: PetscInt newMeetSize = 0;
2108: PetscSectionGetDof(mesh->coneSection, points[p], &dof);
2109: PetscSectionGetOffset(mesh->coneSection, points[p], &off);
2110: for(c = 0; c < dof; ++c) {
2111: const PetscInt point = mesh->cones[off+c];
2113: for(m = 0; m < meetSize; ++m) {
2114: if (point == meet[i][m]) {
2115: meet[1-i][newMeetSize++] = point;
2116: break;
2117: }
2118: }
2119: }
2120: meetSize = newMeetSize;
2121: i = 1-i;
2122: }
2123: *numCoveringPoints = meetSize;
2124: *coveringPoints = meet[i];
2125: return(0);
2126: }
2130: PetscErrorCode DMComplexCreateNeighborCSR(DM dm, PetscInt *numVertices, PetscInt **offsets, PetscInt **adjacency) {
2131: const PetscInt maxFaceCases = 30;
2132: PetscInt numFaceCases = 0;
2133: PetscInt numFaceVertices[30]; /* maxFaceCases, C89 sucks sucks sucks */
2134: PetscInt *off, *adj;
2135: PetscInt *neighborCells, *tmpClosure;
2136: PetscInt maxConeSize, maxSupportSize, maxClosure, maxNeighbors;
2137: PetscInt dim, depth, cStart, cEnd, c, numCells, cell;
2141: /* For parallel partitioning, I think you have to communicate supports */
2142: DMComplexGetDimension(dm, &dim);
2143: DMComplexGetDepth(dm, &depth);
2144: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
2145: DMComplexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
2146: if (cEnd - cStart == 0) {
2147: if (numVertices) *numVertices = 0;
2148: if (offsets) *offsets = PETSC_NULL;
2149: if (adjacency) *adjacency = PETSC_NULL;
2150: return(0);
2151: }
2152: numCells = cEnd - cStart;
2153: /* Setup face recognition */
2154: {
2155: PetscInt cornersSeen[30] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; /* Could use PetscBT */
2157: for(c = cStart; c < cEnd; ++c) {
2158: PetscInt corners;
2160: DMComplexGetConeSize(dm, c, &corners);
2161: if (!cornersSeen[corners]) {
2162: if (numFaceCases >= maxFaceCases) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Exceeded maximum number of face recognition cases");
2163: cornersSeen[corners] = 1;
2164: if (corners == dim+1) {
2165: numFaceVertices[numFaceCases] = dim;
2166: PetscInfo(dm, "Recognizing simplices\n");
2167: } else if ((dim == 1) && (corners == 3)) {
2168: numFaceVertices[numFaceCases] = 3;
2169: PetscInfo(dm, "Recognizing quadratic edges\n");
2170: } else if ((dim == 2) && (corners == 4)) {
2171: numFaceVertices[numFaceCases] = 2;
2172: PetscInfo(dm, "Recognizing quads\n");
2173: } else if ((dim == 2) && (corners == 6)) {
2174: numFaceVertices[numFaceCases] = 3;
2175: PetscInfo(dm, "Recognizing tri and quad cohesive Lagrange cells\n");
2176: } else if ((dim == 2) && (corners == 9)) {
2177: numFaceVertices[numFaceCases] = 3;
2178: PetscInfo(dm, "Recognizing quadratic quads and quadratic quad cohesive Lagrange cells\n");
2179: } else if ((dim == 3) && (corners == 6)) {
2180: numFaceVertices[numFaceCases] = 4;
2181: PetscInfo(dm, "Recognizing tet cohesive cells\n");
2182: } else if ((dim == 3) && (corners == 8)) {
2183: numFaceVertices[numFaceCases] = 4;
2184: PetscInfo(dm, "Recognizing hexes\n");
2185: } else if ((dim == 3) && (corners == 9)) {
2186: numFaceVertices[numFaceCases] = 6;
2187: PetscInfo(dm, "Recognizing tet cohesive Lagrange cells\n");
2188: } else if ((dim == 3) && (corners == 10)) {
2189: numFaceVertices[numFaceCases] = 6;
2190: PetscInfo(dm, "Recognizing quadratic tets\n");
2191: } else if ((dim == 3) && (corners == 12)) {
2192: numFaceVertices[numFaceCases] = 6;
2193: PetscInfo(dm, "Recognizing hex cohesive Lagrange cells\n");
2194: } else if ((dim == 3) && (corners == 18)) {
2195: numFaceVertices[numFaceCases] = 6;
2196: PetscInfo(dm, "Recognizing quadratic tet cohesive Lagrange cells\n");
2197: } else if ((dim == 3) && (corners == 27)) {
2198: numFaceVertices[numFaceCases] = 9;
2199: PetscInfo(dm, "Recognizing quadratic hexes and quadratic hex cohesive Lagrange cells\n");
2200: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Could not recognize number of face vertices for %D corners", corners);
2201: ++numFaceCases;
2202: }
2203: }
2204: }
2205: maxClosure = (PetscInt) (2*PetscMax(pow((PetscReal) maxConeSize, depth)+1, pow((PetscReal) maxSupportSize, depth)+1));
2206: maxNeighbors = (PetscInt) (pow((PetscReal) maxConeSize, depth)*pow((PetscReal) maxSupportSize, depth)+1);
2207: PetscMalloc2(maxNeighbors,PetscInt,&neighborCells,maxClosure,PetscInt,&tmpClosure);
2208: PetscMalloc((numCells+1) * sizeof(PetscInt), &off);
2209: PetscMemzero(off, (numCells+1) * sizeof(PetscInt));
2210: /* Count neighboring cells */
2211: for(cell = cStart; cell < cEnd; ++cell) {
2212: PetscInt numNeighbors = maxNeighbors, n;
2214: DMComplexGetAdjacencySingleLevel_Private(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);
2215: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2216: for(n = 0; n < numNeighbors; ++n) {
2217: PetscInt cellPair[2] = {cell, neighborCells[n]};
2218: PetscBool found = depth > 1 ? PETSC_TRUE : PETSC_FALSE;
2219: PetscInt meetSize;
2220: const PetscInt *meet;
2222: if (cellPair[0] == cellPair[1]) continue;
2223: if (!found) {
2224: DMComplexMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2225: if (meetSize) {
2226: PetscInt f;
2228: for(f = 0; f < numFaceCases; ++f) {
2229: if (numFaceVertices[f] == meetSize) {
2230: found = PETSC_TRUE;
2231: break;
2232: }
2233: }
2234: }
2235: }
2236: if (found) {
2237: ++off[cell-cStart+1];
2238: }
2239: }
2240: }
2241: /* Prefix sum */
2242: for(cell = 1; cell <= numCells; ++cell) {
2243: off[cell] += off[cell-1];
2244: }
2245: if (adjacency) {
2246: PetscMalloc(off[numCells] * sizeof(PetscInt), &adj);
2247: /* Get neighboring cells */
2248: for(cell = cStart; cell < cEnd; ++cell) {
2249: PetscInt numNeighbors = maxNeighbors, n;
2250: PetscInt cellOffset = 0;
2252: DMComplexGetAdjacencySingleLevel_Private(dm, cell, PETSC_TRUE, tmpClosure, &numNeighbors, neighborCells);
2253: /* Get meet with each cell, and check with recognizer (could optimize to check each pair only once) */
2254: for(n = 0; n < numNeighbors; ++n) {
2255: PetscInt cellPair[2] = {cell, neighborCells[n]};
2256: PetscBool found = depth > 1 ? PETSC_TRUE : PETSC_FALSE;
2257: PetscInt meetSize;
2258: const PetscInt *meet;
2260: if (cellPair[0] == cellPair[1]) continue;
2261: if (!found) {
2262: DMComplexMeetPoints(dm, 2, cellPair, &meetSize, &meet);
2263: if (meetSize) {
2264: PetscInt f;
2266: for(f = 0; f < numFaceCases; ++f) {
2267: if (numFaceVertices[f] == meetSize) {
2268: found = PETSC_TRUE;
2269: break;
2270: }
2271: }
2272: }
2273: }
2274: if (found) {
2275: adj[off[cell-cStart]+cellOffset] = neighborCells[n];
2276: ++cellOffset;
2277: }
2278: }
2279: }
2280: }
2281: PetscFree2(neighborCells,tmpClosure);
2282: if (numVertices) *numVertices = numCells;
2283: if (offsets) *offsets = off;
2284: if (adjacency) *adjacency = adj;
2285: return(0);
2286: }
2288: #ifdef PETSC_HAVE_CHACO
2289: #ifdef PETSC_HAVE_UNISTD_H
2290: #include <unistd.h>
2291: #endif
2292: /* Chaco does not have an include file */
2294: float *ewgts, float *x, float *y, float *z, char *outassignname,
2295: char *outfilename, short *assignment, int architecture, int ndims_tot,
2296: int mesh_dims[3], double *goal, int global_method, int local_method,
2297: int rqi_flag, int vmax, int ndims, double eigtol, long seed);
2299: extern int FREE_GRAPH;
2303: PetscErrorCode DMComplexPartition_Chaco(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2304: {
2305: enum {DEFAULT_METHOD = 1, INERTIAL_METHOD = 3};
2306: MPI_Comm comm = ((PetscObject) dm)->comm;
2307: int nvtxs = numVertices; /* number of vertices in full graph */
2308: int *vwgts = NULL; /* weights for all vertices */
2309: float *ewgts = NULL; /* weights for all edges */
2310: float *x = NULL, *y = NULL, *z = NULL; /* coordinates for inertial method */
2311: char *outassignname = NULL; /* name of assignment output file */
2312: char *outfilename = NULL; /* output file name */
2313: int architecture = 1; /* 0 => hypercube, d => d-dimensional mesh */
2314: int ndims_tot = 0; /* total number of cube dimensions to divide */
2315: int mesh_dims[3]; /* dimensions of mesh of processors */
2316: double *goal = NULL; /* desired set sizes for each set */
2317: int global_method = 1; /* global partitioning algorithm */
2318: int local_method = 1; /* local partitioning algorithm */
2319: int rqi_flag = 0; /* should I use RQI/Symmlq eigensolver? */
2320: int vmax = 200; /* how many vertices to coarsen down to? */
2321: int ndims = 1; /* number of eigenvectors (2^d sets) */
2322: double eigtol = 0.001; /* tolerance on eigenvectors */
2323: long seed = 123636512; /* for random graph mutations */
2324: short int *assignment; /* Output partition */
2325: int fd_stdout, fd_pipe[2];
2326: PetscInt *points;
2327: PetscMPIInt commSize;
2328: int i, v, p;
2332: MPI_Comm_size(comm, &commSize);
2333: if (!numVertices) {
2334: PetscSectionCreate(comm, partSection);
2335: PetscSectionSetChart(*partSection, 0, commSize);
2336: PetscSectionSetUp(*partSection);
2337: ISCreateGeneral(comm, 0, PETSC_NULL, PETSC_OWN_POINTER, partition);
2338: return(0);
2339: }
2340: FREE_GRAPH = 0; /* Do not let Chaco free my memory */
2341: for(i = 0; i < start[numVertices]; ++i) {
2342: ++adjacency[i];
2343: }
2344: if (global_method == INERTIAL_METHOD) {
2345: /* manager.createCellCoordinates(nvtxs, &x, &y, &z); */
2346: SETERRQ(comm, PETSC_ERR_SUP, "Inertial partitioning not yet supported");
2347: }
2348: mesh_dims[0] = commSize;
2349: mesh_dims[1] = 1;
2350: mesh_dims[2] = 1;
2351: PetscMalloc(nvtxs * sizeof(short int), &assignment);
2352: /* Chaco outputs to stdout. We redirect this to a buffer. */
2353: /* TODO: check error codes for UNIX calls */
2354: #ifdef PETSC_HAVE_UNISTD_H
2355: {
2356: fd_stdout = dup(1);
2357: pipe(fd_pipe);
2358: close(1);
2359: dup2(fd_pipe[1], 1);
2360: }
2361: #endif
2362: interface(nvtxs, (int *) start, (int *) adjacency, vwgts, ewgts, x, y, z, outassignname, outfilename,
2363: assignment, architecture, ndims_tot, mesh_dims, goal, global_method, local_method, rqi_flag,
2364: vmax, ndims, eigtol, seed);
2365: #ifdef PETSC_HAVE_UNISTD_H
2366: {
2367: char msgLog[10000];
2368: int count;
2370: fflush(stdout);
2371: count = read(fd_pipe[0], msgLog, (10000-1)*sizeof(char));
2372: if (count < 0) count = 0;
2373: msgLog[count] = 0;
2374: close(1);
2375: dup2(fd_stdout, 1);
2376: close(fd_stdout);
2377: close(fd_pipe[0]);
2378: close(fd_pipe[1]);
2379: if (ierr) SETERRQ1(comm, PETSC_ERR_LIB, "Error in Chaco library: %s", msgLog);
2380: }
2381: #endif
2382: /* Convert to PetscSection+IS */
2383: PetscSectionCreate(comm, partSection);
2384: PetscSectionSetChart(*partSection, 0, commSize);
2385: for(v = 0; v < nvtxs; ++v) {
2386: PetscSectionAddDof(*partSection, assignment[v], 1);
2387: }
2388: PetscSectionSetUp(*partSection);
2389: PetscMalloc(nvtxs * sizeof(PetscInt), &points);
2390: for(p = 0, i = 0; p < commSize; ++p) {
2391: for(v = 0; v < nvtxs; ++v) {
2392: if (assignment[v] == p) points[i++] = v;
2393: }
2394: }
2395: if (i != nvtxs) SETERRQ2(comm, PETSC_ERR_PLIB, "Number of points %D should be %D", i, nvtxs);
2396: ISCreateGeneral(comm, nvtxs, points, PETSC_OWN_POINTER, partition);
2397: if (global_method == INERTIAL_METHOD) {
2398: /* manager.destroyCellCoordinates(nvtxs, &x, &y, &z); */
2399: }
2400: PetscFree(assignment);
2401: for(i = 0; i < start[numVertices]; ++i) {
2402: --adjacency[i];
2403: }
2404: return(0);
2405: }
2406: #endif
2408: #ifdef PETSC_HAVE_PARMETIS
2411: PetscErrorCode DMComplexPartition_ParMetis(DM dm, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection *partSection, IS *partition)
2412: {
2414: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "ParMetis not yet supported");
2415: return(0);
2416: }
2417: #endif
2421: PetscErrorCode DMComplexCreatePartition(DM dm, PetscSection *partSection, IS *partition, PetscInt height) {
2422: PetscMPIInt size;
2426: MPI_Comm_size(((PetscObject) dm)->comm, &size);
2427: if (size == 1) {
2428: PetscInt *points;
2429: PetscInt cStart, cEnd, c;
2431: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
2432: PetscSectionCreate(((PetscObject) dm)->comm, partSection);
2433: PetscSectionSetChart(*partSection, 0, size);
2434: PetscSectionSetDof(*partSection, 0, cEnd-cStart);
2435: PetscSectionSetUp(*partSection);
2436: PetscMalloc((cEnd - cStart) * sizeof(PetscInt), &points);
2437: for(c = cStart; c < cEnd; ++c) {
2438: points[c] = c;
2439: }
2440: ISCreateGeneral(((PetscObject) dm)->comm, cEnd-cStart, points, PETSC_OWN_POINTER, partition);
2441: return(0);
2442: }
2443: if (height == 0) {
2444: PetscInt numVertices;
2445: PetscInt *start = PETSC_NULL;
2446: PetscInt *adjacency = PETSC_NULL;
2448: if (1) {
2449: DMComplexCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2450: #ifdef PETSC_HAVE_CHACO
2451: DMComplexPartition_Chaco(dm, numVertices, start, adjacency, partSection, partition);
2452: #endif
2453: } else {
2454: DMComplexCreateNeighborCSR(dm, &numVertices, &start, &adjacency);
2455: #ifdef PETSC_HAVE_PARMETIS
2456: DMComplexPartition_ParMetis(dm, numVertices, start, adjacency, partSection, partition);
2457: #endif
2458: }
2459: PetscFree(start);
2460: PetscFree(adjacency);
2461: # if 0
2462: } else if (height == 1) {
2463: /* Build the dual graph for faces and partition the hypergraph */
2464: PetscInt numEdges;
2466: buildFaceCSRV(mesh, mesh->getFactory()->getNumbering(mesh, mesh->depth()-1), &numEdges, &start, &adjacency, GraphPartitioner::zeroBase());
2467: GraphPartitioner().partition(numEdges, start, adjacency, partition, manager);
2468: destroyCSR(numEdges, start, adjacency);
2469: #endif
2470: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid partition height %D", height);
2471: return(0);
2472: }
2476: PetscErrorCode DMComplexCreatePartitionClosure(DM dm, PetscSection pointSection, IS pointPartition, PetscSection *section, IS *partition) {
2477: /* const PetscInt height = 0; */
2478: const PetscInt *partArray;
2479: PetscInt *allPoints, *partPoints = PETSC_NULL;
2480: PetscInt rStart, rEnd, rank, maxPartSize = 0, newSize;
2481: PetscErrorCode ierr;
2484: PetscSectionGetChart(pointSection, &rStart, &rEnd);
2485: ISGetIndices(pointPartition, &partArray);
2486: PetscSectionCreate(((PetscObject) dm)->comm, section);
2487: PetscSectionSetChart(*section, rStart, rEnd);
2488: for(rank = rStart; rank < rEnd; ++rank) {
2489: PetscInt partSize = 0;
2490: PetscInt numPoints, offset, p;
2492: PetscSectionGetDof(pointSection, rank, &numPoints);
2493: PetscSectionGetOffset(pointSection, rank, &offset);
2494: for(p = 0; p < numPoints; ++p) {
2495: PetscInt point = partArray[offset+p], closureSize, c;
2496: PetscInt *closure = PETSC_NULL;
2498: /* TODO Include support for height > 0 case */
2499: DMComplexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2500: /* Merge into existing points */
2501: if (partSize+closureSize > maxPartSize) {
2502: PetscInt *tmpPoints;
2504: maxPartSize = PetscMax(partSize+closureSize, 2*maxPartSize);
2505: PetscMalloc(maxPartSize * sizeof(PetscInt), &tmpPoints);
2506: PetscMemcpy(tmpPoints, partPoints, partSize * sizeof(PetscInt));
2507: PetscFree(partPoints);
2508: partPoints = tmpPoints;
2509: }
2510: for(c = 0; c < closureSize; ++c) {
2511: partPoints[partSize+c] = closure[c*2];
2512: }
2513: partSize += closureSize;
2514: PetscSortRemoveDupsInt(&partSize, partPoints);
2515: }
2516: PetscSectionSetDof(*section, rank, partSize);
2517: }
2518: PetscSectionSetUp(*section);
2519: PetscSectionGetStorageSize(*section, &newSize);
2520: PetscMalloc(newSize * sizeof(PetscInt), &allPoints);
2522: for(rank = rStart; rank < rEnd; ++rank) {
2523: PetscInt partSize = 0, newOffset;
2524: PetscInt numPoints, offset, p;
2526: PetscSectionGetDof(pointSection, rank, &numPoints);
2527: PetscSectionGetOffset(pointSection, rank, &offset);
2528: for(p = 0; p < numPoints; ++p) {
2529: PetscInt point = partArray[offset+p], closureSize, c;
2530: PetscInt *closure = PETSC_NULL;
2532: /* TODO Include support for height > 0 case */
2533: DMComplexGetTransitiveClosure(dm, point, PETSC_TRUE, &closureSize, &closure);
2534: /* Merge into existing points */
2535: for(c = 0; c < closureSize; ++c) {
2536: partPoints[partSize+c] = closure[c*2];
2537: }
2538: partSize += closureSize;
2539: PetscSortRemoveDupsInt(&partSize, partPoints);
2540: }
2541: PetscSectionGetOffset(*section, rank, &newOffset);
2542: PetscMemcpy(&allPoints[newOffset], partPoints, partSize * sizeof(PetscInt));
2543: }
2544: ISRestoreIndices(pointPartition, &partArray);
2545: PetscFree(partPoints);
2546: ISCreateGeneral(((PetscObject) dm)->comm, newSize, allPoints, PETSC_OWN_POINTER, partition);
2547: return(0);
2548: }
2552: /*
2553: Input Parameters:
2554: . originalSection
2555: , originalVec
2557: Output Parameters:
2558: . newSection
2559: . newVec
2560: */
2561: PetscErrorCode DMComplexDistributeField(DM dm, PetscSF pointSF, PetscSection originalSection, Vec originalVec, PetscSection newSection, Vec newVec)
2562: {
2563: PetscSF fieldSF;
2564: PetscInt *remoteOffsets, fieldSize;
2565: PetscScalar *originalValues, *newValues;
2566: PetscErrorCode ierr;
2569: PetscSFDistributeSection(pointSF, originalSection, &remoteOffsets, newSection);
2571: PetscSectionGetStorageSize(newSection, &fieldSize);
2572: VecSetSizes(newVec, fieldSize, PETSC_DETERMINE);
2573: VecSetFromOptions(newVec);
2575: VecGetArray(originalVec, &originalValues);
2576: VecGetArray(newVec, &newValues);
2577: PetscSFCreateSectionSF(pointSF, originalSection, remoteOffsets, newSection, &fieldSF);
2578: PetscSFBcastBegin(fieldSF, MPIU_SCALAR, originalValues, newValues);
2579: PetscSFBcastEnd(fieldSF, MPIU_SCALAR, originalValues, newValues);
2580: PetscSFDestroy(&fieldSF);
2581: VecRestoreArray(newVec, &newValues);
2582: VecRestoreArray(originalVec, &originalValues);
2583: return(0);
2584: }
2588: /*@C
2589: DMComplexDistribute - Distributes the mesh and any associated sections.
2591: Not Collective
2593: Input Parameter:
2594: + dm - The original DMComplex object
2595: - partitioner - The partitioning package, or NULL for the default
2597: Output Parameter:
2598: . parallelMesh - The distributed DMComplex object, or PETSC_NULL
2600: Note: If the mesh was not distributed, the return value is PETSC_NULL
2602: Level: intermediate
2604: .keywords: mesh, elements
2605: .seealso: DMComplexCreate(), DMComplexDistributeByFace()
2606: @*/
2607: PetscErrorCode DMComplexDistribute(DM dm, const char partitioner[], DM *dmParallel)
2608: {
2609: DM_Complex *mesh = (DM_Complex *) dm->data, *pmesh;
2610: MPI_Comm comm = ((PetscObject) dm)->comm;
2611: const PetscInt height = 0;
2612: PetscInt dim, numRemoteRanks;
2613: IS cellPart, part;
2614: PetscSection cellPartSection, partSection;
2615: PetscSFNode *remoteRanks;
2616: PetscSF partSF, pointSF, coneSF;
2617: ISLocalToGlobalMapping renumbering;
2618: PetscSection originalConeSection, newConeSection;
2619: PetscInt *remoteOffsets;
2620: PetscInt *cones, *newCones, newConesSize;
2621: PetscBool flg;
2622: PetscMPIInt rank, numProcs, p;
2628: MPI_Comm_rank(comm, &rank);
2629: MPI_Comm_size(comm, &numProcs);
2630: if (numProcs == 1) return(0);
2632: DMComplexGetDimension(dm, &dim);
2633: /* Create cell partition - We need to rewrite to use IS, use the MatPartition stuff */
2634: DMComplexCreatePartition(dm, &cellPartSection, &cellPart, height);
2635: /* Create SF assuming a serial partition for all processes: Could check for IS length here */
2636: if (!rank) {
2637: numRemoteRanks = numProcs;
2638: } else {
2639: numRemoteRanks = 0;
2640: }
2641: PetscMalloc(numRemoteRanks * sizeof(PetscSFNode), &remoteRanks);
2642: for(p = 0; p < numRemoteRanks; ++p) {
2643: remoteRanks[p].rank = p;
2644: remoteRanks[p].index = 0;
2645: }
2646: PetscSFCreate(comm, &partSF);
2647: PetscSFSetGraph(partSF, 1, numRemoteRanks, PETSC_NULL, PETSC_OWN_POINTER, remoteRanks, PETSC_OWN_POINTER);
2648: PetscOptionsHasName(((PetscObject) dm)->prefix, "-partition_view", &flg);
2649: if (flg) {
2650: PetscPrintf(comm, "Cell Partition:\n");
2651: PetscSectionView(cellPartSection, PETSC_VIEWER_STDOUT_WORLD);
2652: ISView(cellPart, PETSC_NULL);
2653: PetscSFView(partSF, PETSC_NULL);
2654: }
2655: /* Close the partition over the mesh */
2656: DMComplexCreatePartitionClosure(dm, cellPartSection, cellPart, &partSection, &part);
2657: ISDestroy(&cellPart);
2658: PetscSectionDestroy(&cellPartSection);
2659: /* Create new mesh */
2660: DMComplexCreate(comm, dmParallel);
2661: DMComplexSetDimension(*dmParallel, dim);
2662: PetscObjectSetName((PetscObject) *dmParallel, "Parallel Mesh");
2663: pmesh = (DM_Complex *) (*dmParallel)->data;
2664: /* Distribute sieve points and the global point numbering (replaces creating remote bases) */
2665: PetscSFConvertPartition(partSF, partSection, part, &renumbering, &pointSF);
2666: if (flg) {
2667: PetscPrintf(comm, "Point Partition:\n");
2668: PetscSectionView(partSection, PETSC_VIEWER_STDOUT_WORLD);
2669: ISView(part, PETSC_NULL);
2670: PetscSFView(pointSF, PETSC_NULL);
2671: PetscPrintf(comm, "Point Renumbering after partition:\n");
2672: ISLocalToGlobalMappingView(renumbering, PETSC_NULL);
2673: }
2674: /* Distribute cone section */
2675: DMComplexGetConeSection(dm, &originalConeSection);
2676: DMComplexGetConeSection(*dmParallel, &newConeSection);
2677: PetscSFDistributeSection(pointSF, originalConeSection, &remoteOffsets, newConeSection);
2678: DMSetUp(*dmParallel);
2679: {
2680: PetscInt pStart, pEnd, p;
2682: PetscSectionGetChart(newConeSection, &pStart, &pEnd);
2683: for(p = pStart; p < pEnd; ++p) {
2684: PetscInt coneSize;
2685: PetscSectionGetDof(newConeSection, p, &coneSize);
2686: pmesh->maxConeSize = PetscMax(pmesh->maxConeSize, coneSize);
2687: }
2688: }
2689: /* Communicate and renumber cones */
2690: PetscSFCreateSectionSF(pointSF, originalConeSection, remoteOffsets, newConeSection, &coneSF);
2691: DMComplexGetCones(dm, &cones);
2692: DMComplexGetCones(*dmParallel, &newCones);
2693: PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones);
2694: PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones);
2695: PetscSectionGetStorageSize(newConeSection, &newConesSize);
2696: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newConesSize, newCones, PETSC_NULL, newCones);
2697: PetscOptionsHasName(((PetscObject) dm)->prefix, "-cones_view", &flg);
2698: if (flg) {
2699: PetscPrintf(comm, "Serial Cone Section:\n");
2700: PetscSectionView(originalConeSection, PETSC_VIEWER_STDOUT_WORLD);
2701: PetscPrintf(comm, "Parallel Cone Section:\n");
2702: PetscSectionView(newConeSection, PETSC_VIEWER_STDOUT_WORLD);
2703: PetscSFView(coneSF, PETSC_NULL);
2704: }
2705: DMComplexGetConeOrientations(dm, &cones);
2706: DMComplexGetConeOrientations(*dmParallel, &newCones);
2707: PetscSFBcastBegin(coneSF, MPIU_INT, cones, newCones);
2708: PetscSFBcastEnd(coneSF, MPIU_INT, cones, newCones);
2709: PetscSFDestroy(&coneSF);
2710: /* Create supports and stratify sieve */
2711: DMComplexSymmetrize(*dmParallel);
2712: DMComplexStratify(*dmParallel);
2713: /* Distribute Coordinates */
2714: {
2715: PetscSection originalCoordSection, newCoordSection;
2716: Vec originalCoordinates, newCoordinates;
2718: DMComplexGetCoordinateSection(dm, &originalCoordSection);
2719: DMComplexGetCoordinateSection(*dmParallel, &newCoordSection);
2720: DMComplexGetCoordinateVec(dm, &originalCoordinates);
2721: DMComplexGetCoordinateVec(*dmParallel, &newCoordinates);
2723: DMComplexDistributeField(dm, pointSF, originalCoordSection, originalCoordinates, newCoordSection, newCoordinates);
2724: }
2725: /* Distribute labels */
2726: {
2727: SieveLabel next = mesh->labels, newNext = PETSC_NULL;
2728: PetscInt numLabels = 0, l;
2730: /* Bcast number of labels */
2731: while(next) {++numLabels; next = next->next;}
2732: MPI_Bcast(&numLabels, 1, MPIU_INT, 0, comm);
2733: next = mesh->labels;
2734: for(l = 0; l < numLabels; ++l) {
2735: SieveLabel newLabel;
2736: const PetscInt *partArray;
2737: PetscInt *stratumSizes = PETSC_NULL, *points = PETSC_NULL;
2738: PetscMPIInt *sendcnts = PETSC_NULL, *offsets = PETSC_NULL, *displs = PETSC_NULL;
2739: PetscInt nameSize, s, p;
2740: size_t len = 0;
2742: PetscNew(struct Sieve_Label, &newLabel);
2743: /* Bcast name (could filter for no points) */
2744: if (!rank) {PetscStrlen(next->name, &len);}
2745: nameSize = len;
2746: MPI_Bcast(&nameSize, 1, MPIU_INT, 0, comm);
2747: PetscMalloc(nameSize+1, &newLabel->name);
2748: if (!rank) {PetscMemcpy(newLabel->name, next->name, nameSize+1);}
2749: MPI_Bcast(newLabel->name, nameSize+1, MPI_CHAR, 0, comm);
2750: /* Bcast numStrata (could filter for no points in stratum) */
2751: if (!rank) {newLabel->numStrata = next->numStrata;}
2752: MPI_Bcast(&newLabel->numStrata, 1, MPIU_INT, 0, comm);
2753: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumValues);
2754: PetscMalloc(newLabel->numStrata * sizeof(PetscInt), &newLabel->stratumSizes);
2755: PetscMalloc((newLabel->numStrata+1) * sizeof(PetscInt), &newLabel->stratumOffsets);
2756: /* Bcast stratumValues (could filter for no points in stratum) */
2757: if (!rank) {PetscMemcpy(newLabel->stratumValues, next->stratumValues, next->numStrata * sizeof(PetscInt));}
2758: MPI_Bcast(newLabel->stratumValues, newLabel->numStrata, MPIU_INT, 0, comm);
2759: /* Find size on each process and Scatter */
2760: if (!rank) {
2761: ISGetIndices(part, &partArray);
2762: PetscMalloc(numProcs*next->numStrata * sizeof(PetscInt), &stratumSizes);
2763: PetscMemzero(stratumSizes, numProcs*next->numStrata * sizeof(PetscInt));
2764: for(s = 0; s < next->numStrata; ++s) {
2765: for(p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
2766: const PetscInt point = next->points[p];
2767: PetscInt proc;
2769: for(proc = 0; proc < numProcs; ++proc) {
2770: PetscInt dof, off, pPart;
2772: PetscSectionGetDof(partSection, proc, &dof);
2773: PetscSectionGetOffset(partSection, proc, &off);
2774: for(pPart = off; pPart < off+dof; ++pPart) {
2775: if (partArray[pPart] == point) {
2776: ++stratumSizes[proc*next->numStrata+s];
2777: break;
2778: }
2779: }
2780: }
2781: }
2782: }
2783: ISRestoreIndices(part, &partArray);
2784: }
2785: MPI_Scatter(stratumSizes, newLabel->numStrata, MPI_INT, newLabel->stratumSizes, newLabel->numStrata, MPI_INT, 0, comm);
2786: /* Calculate stratumOffsets */
2787: newLabel->stratumOffsets[0] = 0;
2788: for(s = 0; s < newLabel->numStrata; ++s) {
2789: newLabel->stratumOffsets[s+1] = newLabel->stratumSizes[s] + newLabel->stratumOffsets[s];
2790: }
2791: /* Pack points and Scatter */
2792: if (!rank) {
2793: PetscMalloc3(numProcs,PetscMPIInt,&sendcnts,numProcs,PetscMPIInt,&offsets,numProcs+1,PetscMPIInt,&displs);
2794: displs[0] = 0;
2795: for(p = 0; p < numProcs; ++p) {
2796: sendcnts[p] = 0;
2797: for(s = 0; s < next->numStrata; ++s) {
2798: sendcnts[p] += stratumSizes[p*next->numStrata+s];
2799: }
2800: offsets[p] = displs[p];
2801: displs[p+1] = displs[p] + sendcnts[p];
2802: }
2803: PetscMalloc(displs[numProcs] * sizeof(PetscInt), &points);
2804: for(s = 0; s < next->numStrata; ++s) {
2805: for(p = next->stratumOffsets[s]; p < next->stratumOffsets[s]+next->stratumSizes[s]; ++p) {
2806: const PetscInt point = next->points[p];
2807: PetscInt proc;
2809: for(proc = 0; proc < numProcs; ++proc) {
2810: PetscInt dof, off, pPart;
2812: PetscSectionGetDof(partSection, proc, &dof);
2813: PetscSectionGetOffset(partSection, proc, &off);
2814: for(pPart = off; pPart < off+dof; ++pPart) {
2815: if (partArray[pPart] == point) {
2816: points[offsets[proc]++] = point;
2817: break;
2818: }
2819: }
2820: }
2821: }
2822: }
2823: }
2824: PetscMalloc(newLabel->stratumOffsets[newLabel->numStrata] * sizeof(PetscInt), &newLabel->points);
2825: MPI_Scatterv(points, sendcnts, displs, MPIU_INT, newLabel->points, newLabel->stratumOffsets[newLabel->numStrata], MPIU_INT, 0, comm);
2826: PetscFree(points);
2827: PetscFree3(sendcnts,offsets,displs);
2828: PetscFree(stratumSizes);
2829: /* Renumber points */
2830: ISGlobalToLocalMappingApply(renumbering, IS_GTOLM_MASK, newLabel->stratumOffsets[newLabel->numStrata], newLabel->points, PETSC_NULL, newLabel->points);
2831: /* Sort points */
2832: for(s = 0; s < newLabel->numStrata; ++s) {
2833: PetscSortInt(newLabel->stratumSizes[s], &newLabel->points[newLabel->stratumOffsets[s]]);
2834: }
2835: /* Insert into list */
2836: if (newNext) {
2837: newNext->next = newLabel;
2838: } else {
2839: pmesh->labels = newLabel;
2840: }
2841: newNext = newLabel;
2842: if (!rank) {next = next->next;}
2843: }
2844: }
2845: /* Cleanup Partition */
2846: ISLocalToGlobalMappingDestroy(&renumbering);
2847: PetscSFDestroy(&partSF);
2848: PetscSectionDestroy(&partSection);
2849: ISDestroy(&part);
2850: /* Create point SF for parallel mesh */
2851: {
2852: const PetscInt *leaves;
2853: PetscSFNode *remotePoints;
2854: PetscInt *rowners, *lowners, *ghostPoints;
2855: PetscInt numRoots, numLeaves, numGhostPoints = 0, p, gp;
2856: PetscInt pStart, pEnd;
2858: DMComplexGetChart(*dmParallel, &pStart, &pEnd);
2859: PetscSFGetGraph(pointSF, &numRoots, &numLeaves, &leaves, PETSC_NULL);
2860: PetscMalloc2(numRoots*2,PetscInt,&rowners,numLeaves*2,PetscInt,&lowners);
2861: for(p = 0; p < numRoots*2; ++p) {
2862: rowners[p] = 0;
2863: }
2864: for(p = 0; p < numLeaves; ++p) {
2865: lowners[p*2+0] = rank;
2866: lowners[p*2+1] = leaves ? leaves[p] : p;
2867: }
2868: #if 0 /* Why doesn't this datatype work */
2869: PetscSFFetchAndOpBegin(pointSF, MPIU_2INT, rowners, lowners, lowners, MPI_MAXLOC);
2870: PetscSFFetchAndOpEnd(pointSF, MPIU_2INT, rowners, lowners, lowners, MPI_MAXLOC);
2871: #endif
2872: PetscSFFetchAndOpBegin(pointSF, MPI_2INT, rowners, lowners, lowners, MPI_MAXLOC);
2873: PetscSFFetchAndOpEnd(pointSF, MPI_2INT, rowners, lowners, lowners, MPI_MAXLOC);
2874: PetscSFBcastBegin(pointSF, MPIU_2INT, rowners, lowners);
2875: PetscSFBcastEnd(pointSF, MPIU_2INT, rowners, lowners);
2876: for(p = 0; p < numLeaves; ++p) {
2877: if (lowners[p*2+0] != rank) ++numGhostPoints;
2878: }
2879: PetscMalloc(numGhostPoints * sizeof(PetscInt), &ghostPoints);
2880: PetscMalloc(numGhostPoints * sizeof(PetscSFNode), &remotePoints);
2881: for(p = 0, gp = 0; p < numLeaves; ++p) {
2882: if (lowners[p*2+0] != rank) {
2883: ghostPoints[gp] = leaves ? leaves[p] : p;
2884: remotePoints[gp].rank = lowners[p*2+0];
2885: remotePoints[gp].index = lowners[p*2+1];
2886: ++gp;
2887: }
2888: }
2889: PetscFree2(rowners,lowners);
2890: PetscSFSetGraph((*dmParallel)->sf, pEnd - pStart, numGhostPoints, ghostPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);
2891: PetscSFSetFromOptions((*dmParallel)->sf);
2892: }
2893: /* Cleanup */
2894: PetscSFDestroy(&pointSF);
2895: DMSetFromOptions(*dmParallel);
2896: return(0);
2897: }
2899: #ifdef PETSC_HAVE_TRIANGLE
2900: #include <triangle.h>
2904: PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx) {
2906: inputCtx->numberofpoints = 0;
2907: inputCtx->numberofpointattributes = 0;
2908: inputCtx->pointlist = PETSC_NULL;
2909: inputCtx->pointattributelist = PETSC_NULL;
2910: inputCtx->pointmarkerlist = PETSC_NULL;
2911: inputCtx->numberofsegments = 0;
2912: inputCtx->segmentlist = PETSC_NULL;
2913: inputCtx->segmentmarkerlist = PETSC_NULL;
2914: inputCtx->numberoftriangleattributes = 0;
2915: inputCtx->trianglelist = PETSC_NULL;
2916: inputCtx->numberofholes = 0;
2917: inputCtx->holelist = PETSC_NULL;
2918: inputCtx->numberofregions = 0;
2919: inputCtx->regionlist = PETSC_NULL;
2920: return(0);
2921: }
2925: PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx) {
2927: outputCtx->numberofpoints = 0;
2928: outputCtx->pointlist = PETSC_NULL;
2929: outputCtx->pointattributelist = PETSC_NULL;
2930: outputCtx->pointmarkerlist = PETSC_NULL;
2931: outputCtx->numberoftriangles = 0;
2932: outputCtx->trianglelist = PETSC_NULL;
2933: outputCtx->triangleattributelist = PETSC_NULL;
2934: outputCtx->neighborlist = PETSC_NULL;
2935: outputCtx->segmentlist = PETSC_NULL;
2936: outputCtx->segmentmarkerlist = PETSC_NULL;
2937: outputCtx->numberofedges = 0;
2938: outputCtx->edgelist = PETSC_NULL;
2939: outputCtx->edgemarkerlist = PETSC_NULL;
2940: return(0);
2941: }
2945: PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx) {
2947: free(outputCtx->pointmarkerlist);
2948: free(outputCtx->edgelist);
2949: free(outputCtx->edgemarkerlist);
2950: free(outputCtx->trianglelist);
2951: free(outputCtx->neighborlist);
2952: return(0);
2953: }
2957: PetscErrorCode DMComplexInterpolate_2D(DM dm, DM *dmInt)
2958: {
2959: DM idm;
2960: DM_Complex *mesh;
2961: PetscInt *off;
2962: const PetscInt numCorners = 3;
2963: PetscInt dim, numCells, cStart, cEnd, c, numVertices, vStart, vEnd;
2964: PetscInt numEdges, firstEdge, edge, e;
2968: DMComplexGetDimension(dm, &dim);
2969: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
2970: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
2971: numCells = cEnd - cStart;
2972: numVertices = vEnd - vStart;
2973: firstEdge = numCells + numVertices;
2974: numEdges = 0 ;
2975: /* Count edges using algorithm from CreateNeighborCSR */
2976: DMComplexCreateNeighborCSR(dm, PETSC_NULL, &off, PETSC_NULL);
2977: if (off) {
2978: numEdges = off[numCells]/2;
2979: /* Account for boundary edges: \sum_c 3 - neighbors = 3*numCells - totalNeighbors */
2980: numEdges += 3*numCells - off[numCells];
2981: }
2982: /* Create interpolated mesh */
2983: DMCreate(((PetscObject) dm)->comm, &idm);
2984: DMSetType(idm, DMCOMPLEX);
2985: DMComplexSetDimension(idm, dim);
2986: DMComplexSetChart(idm, 0, numCells+numVertices+numEdges);
2987: for(c = 0; c < numCells; ++c) {
2988: DMComplexSetConeSize(idm, c, numCorners);
2989: }
2990: for(e = firstEdge; e < firstEdge+numEdges; ++e) {
2991: DMComplexSetConeSize(idm, e, 2);
2992: }
2993: DMSetUp(idm);
2994: for(c = 0, edge = firstEdge; c < numCells; ++c) {
2995: const PetscInt *faces;
2996: PetscInt numFaces, faceSize, f;
2998: DMComplexGetFaces(dm, c, &numFaces, &faceSize, &faces);
2999: if (faceSize != 2) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize);
3000: for(f = 0; f < numFaces; ++f) {
3001: PetscBool found = PETSC_FALSE;
3003: /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */
3004: for(e = firstEdge; e < edge; ++e) {
3005: const PetscInt *cone;
3007: DMComplexGetCone(idm, e, &cone);
3008: if (((faces[f*faceSize+0] == cone[0]) && (faces[f*faceSize+1] == cone[1])) ||
3009: ((faces[f*faceSize+0] == cone[1]) && (faces[f*faceSize+1] == cone[0]))) {
3010: found = PETSC_TRUE;
3011: break;
3012: }
3013: }
3014: if (!found) {
3015: DMComplexSetCone(idm, edge, &faces[f*faceSize]);
3016: ++edge;
3017: }
3018: DMComplexInsertCone(idm, c, f, e);
3019: }
3020: }
3021: if (edge != firstEdge+numEdges) SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges);
3022: PetscFree(off);
3023: DMComplexSymmetrize(idm);
3024: DMComplexStratify(idm);
3025: mesh = (DM_Complex *) (idm)->data;
3026: for(c = 0; c < numCells; ++c) {
3027: const PetscInt *cone, *faces;
3028: PetscInt coneSize, coff, numFaces, faceSize, f;
3030: DMComplexGetConeSize(idm, c, &coneSize);
3031: DMComplexGetCone(idm, c, &cone);
3032: PetscSectionGetOffset(mesh->coneSection, c, &coff);
3033: DMComplexGetFaces(dm, c, &numFaces, &faceSize, &faces);
3034: if (coneSize != numFaces) SETERRQ3(((PetscObject) idm)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numFaces);
3035: for(f = 0; f < numFaces; ++f) {
3036: const PetscInt *econe;
3037: PetscInt esize;
3039: DMComplexGetConeSize(idm, cone[f], &esize);
3040: DMComplexGetCone(idm, cone[f], &econe);
3041: if (esize != 2) SETERRQ2(((PetscObject) idm)->comm, PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[f]);
3042: if ((faces[f*faceSize+0] == econe[0]) && (faces[f*faceSize+1] == econe[1])) {
3043: /* Correctly oriented */
3044: mesh->coneOrientations[coff+f] = 0;
3045: } else if ((faces[f*faceSize+0] == econe[1]) && (faces[f*faceSize+1] == econe[0])) {
3046: /* Start at index 1, and reverse orientation */
3047: mesh->coneOrientations[coff+f] = -(1+1);
3048: }
3049: }
3050: }
3051: *dmInt = idm;
3052: return(0);
3053: }
3057: PetscErrorCode DMComplexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm)
3058: {
3059: MPI_Comm comm = ((PetscObject) boundary)->comm;
3060: DM_Complex *bd = (DM_Complex *) boundary->data;
3061: PetscInt dim = 2;
3062: const PetscBool createConvexHull = PETSC_FALSE;
3063: const PetscBool constrained = PETSC_FALSE;
3064: struct triangulateio in;
3065: struct triangulateio out;
3066: PetscInt vStart, vEnd, v, eStart, eEnd, e;
3067: PetscMPIInt rank;
3068: PetscErrorCode ierr;
3071: MPI_Comm_rank(comm, &rank);
3072: InitInput_Triangle(&in);
3073: InitOutput_Triangle(&out);
3074: DMComplexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3075: in.numberofpoints = vEnd - vStart;
3076: if (in.numberofpoints > 0) {
3077: PetscScalar *array;
3079: PetscMalloc(in.numberofpoints*dim * sizeof(double), &in.pointlist);
3080: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
3081: VecGetArray(bd->coordinates, &array);
3082: for(v = vStart; v < vEnd; ++v) {
3083: const PetscInt idx = v - vStart;
3084: PetscInt off, d;
3086: PetscSectionGetOffset(bd->coordSection, v, &off);
3087: for(d = 0; d < dim; ++d) {
3088: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3089: }
3090: DMComplexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3091: }
3092: VecRestoreArray(bd->coordinates, &array);
3093: }
3094: DMComplexGetHeightStratum(boundary, 0, &eStart, &eEnd);
3095: in.numberofsegments = eEnd - eStart;
3096: if (in.numberofsegments > 0) {
3097: PetscMalloc(in.numberofsegments*2 * sizeof(int), &in.segmentlist);
3098: PetscMalloc(in.numberofsegments * sizeof(int), &in.segmentmarkerlist);
3099: for(e = eStart; e < eEnd; ++e) {
3100: const PetscInt idx = e - eStart;
3101: const PetscInt *cone;
3103: DMComplexGetCone(boundary, e, &cone);
3104: in.segmentlist[idx*2+0] = cone[0] - vStart;
3105: in.segmentlist[idx*2+1] = cone[1] - vStart;
3106: DMComplexGetLabelValue(boundary, "marker", e, &in.segmentmarkerlist[idx]);
3107: }
3108: }
3109: #if 0 /* Do not currently support holes */
3110: PetscReal *holeCoords;
3111: PetscInt h, d;
3113: DMComplexGetHoles(boundary, &in.numberofholes, &holeCords);
3114: if (in.numberofholes > 0) {
3115: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
3116: for(h = 0; h < in.numberofholes; ++h) {
3117: for(d = 0; d < dim; ++d) {
3118: in.holelist[h*dim+d] = holeCoords[h*dim+d];
3119: }
3120: }
3121: }
3122: #endif
3123: if (!rank) {
3124: char args[32];
3126: /* Take away 'Q' for verbose output */
3127: PetscStrcpy(args, "pqezQ");
3128: if (createConvexHull) {
3129: PetscStrcat(args, "c");
3130: }
3131: if (constrained) {
3132: PetscStrcpy(args, "zepDQ");
3133: }
3134: triangulate(args, &in, &out, PETSC_NULL);
3135: }
3136: PetscFree(in.pointlist);
3137: PetscFree(in.pointmarkerlist);
3138: PetscFree(in.segmentlist);
3139: PetscFree(in.segmentmarkerlist);
3140: PetscFree(in.holelist);
3142: DMCreate(comm, dm);
3143: DMSetType(*dm, DMCOMPLEX);
3144: DMComplexSetDimension(*dm, dim);
3145: {
3146: DM_Complex *mesh = (DM_Complex *) (*dm)->data;
3147: const PetscInt numCorners = 3;
3148: const PetscInt numCells = out.numberoftriangles;
3149: const PetscInt numVertices = out.numberofpoints;
3150: int *cells = out.trianglelist;
3151: double *meshCoords = out.pointlist;
3152: PetscInt coordSize, c;
3153: PetscScalar *coords;
3155: DMComplexSetChart(*dm, 0, numCells+numVertices);
3156: for(c = 0; c < numCells; ++c) {
3157: DMComplexSetConeSize(*dm, c, numCorners);
3158: }
3159: DMSetUp(*dm);
3160: for(c = 0; c < numCells; ++c) {
3161: /* Should be numCorners, but c89 sucks shit */
3162: PetscInt cone[3] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells};
3164: DMComplexSetCone(*dm, c, cone);
3165: }
3166: DMComplexSymmetrize(*dm);
3167: DMComplexStratify(*dm);
3168: if (interpolate) {
3169: DM idm;
3171: DMComplexInterpolate_2D(*dm, &idm);
3172: DMDestroy(dm);
3173: *dm = idm;
3174: mesh = (DM_Complex *) (idm)->data;
3175: }
3176: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3177: for(v = numCells; v < numCells+numVertices; ++v) {
3178: PetscSectionSetDof(mesh->coordSection, v, dim);
3179: }
3180: PetscSectionSetUp(mesh->coordSection);
3181: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3182: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3183: VecSetFromOptions(mesh->coordinates);
3184: VecGetArray(mesh->coordinates, &coords);
3185: for(v = 0; v < numVertices; ++v) {
3186: coords[v*dim+0] = meshCoords[v*dim+0];
3187: coords[v*dim+1] = meshCoords[v*dim+1];
3188: }
3189: VecRestoreArray(mesh->coordinates, &coords);
3190: for(v = 0; v < numVertices; ++v) {
3191: if (out.pointmarkerlist[v]) {
3192: DMComplexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3193: }
3194: }
3195: if (interpolate) {
3196: for(e = 0; e < out.numberofedges; e++) {
3197: if (out.edgemarkerlist[e]) {
3198: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3199: const PetscInt *edges;
3200: PetscInt numEdges;
3202: DMComplexJoinPoints(*dm, 2, vertices, &numEdges, &edges);
3203: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3204: DMComplexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
3205: }
3206: }
3207: }
3208: }
3209: #if 0 /* Do not currently support holes */
3210: DMComplexCopyHoles(*dm, boundary);
3211: #endif
3212: FiniOutput_Triangle(&out);
3213: return(0);
3214: }
3218: PetscErrorCode DMComplexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined)
3219: {
3220: MPI_Comm comm = ((PetscObject) dm)->comm;
3221: DM_Complex *mesh = (DM_Complex *) dm->data;
3222: PetscInt dim = 2;
3223: struct triangulateio in;
3224: struct triangulateio out;
3225: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3226: PetscMPIInt rank;
3227: PetscErrorCode ierr;
3230: MPI_Comm_rank(comm, &rank);
3231: InitInput_Triangle(&in);
3232: InitOutput_Triangle(&out);
3233: DMComplexGetDepth(dm, &depth);
3234: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3235: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
3236: in.numberofpoints = vEnd - vStart;
3237: if (in.numberofpoints > 0) {
3238: PetscScalar *array;
3240: PetscMalloc(in.numberofpoints*dim * sizeof(double), &in.pointlist);
3241: PetscMalloc(in.numberofpoints * sizeof(int), &in.pointmarkerlist);
3242: VecGetArray(mesh->coordinates, &array);
3243: for(v = vStart; v < vEnd; ++v) {
3244: const PetscInt idx = v - vStart;
3245: PetscInt off, d;
3247: PetscSectionGetOffset(mesh->coordSection, v, &off);
3248: for(d = 0; d < dim; ++d) {
3249: in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3250: }
3251: DMComplexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
3252: }
3253: VecRestoreArray(mesh->coordinates, &array);
3254: }
3255: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
3256: in.numberofcorners = 3;
3257: in.numberoftriangles = cEnd - cStart;
3258: in.trianglearealist = (double *) maxVolumes;
3259: if (in.numberoftriangles > 0) {
3260: PetscMalloc(in.numberoftriangles*in.numberofcorners * sizeof(int), &in.trianglelist);
3261: for(c = cStart; c < cEnd; ++c) {
3262: const PetscInt idx = c - cStart;
3263: PetscInt *closure = PETSC_NULL;
3264: PetscInt closureSize;
3266: DMComplexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3267: if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize);
3268: for(v = 0; v < 3; ++v) {
3269: in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart;
3270: }
3271: }
3272: }
3273: #if 0 /* Do not currently support holes */
3274: PetscReal *holeCoords;
3275: PetscInt h, d;
3277: DMComplexGetHoles(boundary, &in.numberofholes, &holeCords);
3278: if (in.numberofholes > 0) {
3279: PetscMalloc(in.numberofholes*dim * sizeof(double), &in.holelist);
3280: for(h = 0; h < in.numberofholes; ++h) {
3281: for(d = 0; d < dim; ++d) {
3282: in.holelist[h*dim+d] = holeCoords[h*dim+d];
3283: }
3284: }
3285: }
3286: #endif
3287: if (!rank) {
3288: char args[32];
3290: /* Take away 'Q' for verbose output */
3291: PetscStrcpy(args, "pqezQra");
3292: triangulate(args, &in, &out, PETSC_NULL);
3293: }
3294: PetscFree(in.pointlist);
3295: PetscFree(in.pointmarkerlist);
3296: PetscFree(in.segmentlist);
3297: PetscFree(in.segmentmarkerlist);
3298: PetscFree(in.trianglelist);
3300: DMCreate(comm, dmRefined);
3301: DMSetType(*dmRefined, DMCOMPLEX);
3302: DMComplexSetDimension(*dmRefined, dim);
3303: {
3304: DM_Complex *mesh = (DM_Complex *) (*dmRefined)->data;
3305: const PetscInt numCorners = 3;
3306: const PetscInt numCells = out.numberoftriangles;
3307: const PetscInt numVertices = out.numberofpoints;
3308: int *cells = out.trianglelist;
3309: double *meshCoords = out.pointlist;
3310: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
3311: PetscInt coordSize, c, e;
3312: PetscScalar *coords;
3314: DMComplexSetChart(*dmRefined, 0, numCells+numVertices);
3315: for(c = 0; c < numCells; ++c) {
3316: DMComplexSetConeSize(*dmRefined, c, numCorners);
3317: }
3318: DMSetUp(*dmRefined);
3319: for(c = 0; c < numCells; ++c) {
3320: /* Should be numCorners, but c89 sucks shit */
3321: PetscInt cone[3] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells};
3323: DMComplexSetCone(*dmRefined, c, cone);
3324: }
3325: DMComplexSymmetrize(*dmRefined);
3326: DMComplexStratify(*dmRefined);
3328: if (interpolate) {
3329: DM idm;
3331: DMComplexInterpolate_2D(*dmRefined, &idm);
3332: DMDestroy(dmRefined);
3333: *dmRefined = idm;
3334: mesh = (DM_Complex *) (idm)->data;
3335: }
3336: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3337: for(v = numCells; v < numCells+numVertices; ++v) {
3338: PetscSectionSetDof(mesh->coordSection, v, dim);
3339: }
3340: PetscSectionSetUp(mesh->coordSection);
3341: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3342: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3343: VecSetFromOptions(mesh->coordinates);
3344: VecGetArray(mesh->coordinates, &coords);
3345: for(v = 0; v < numVertices; ++v) {
3346: coords[v*dim+0] = meshCoords[v*dim+0];
3347: coords[v*dim+1] = meshCoords[v*dim+1];
3348: }
3349: VecRestoreArray(mesh->coordinates, &coords);
3350: for(v = 0; v < numVertices; ++v) {
3351: if (out.pointmarkerlist[v]) {
3352: DMComplexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3353: }
3354: }
3355: if (interpolate) {
3356: for(e = 0; e < out.numberofedges; e++) {
3357: if (out.edgemarkerlist[e]) {
3358: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3359: const PetscInt *edges;
3360: PetscInt numEdges;
3362: DMComplexJoinPoints(*dmRefined, 2, vertices, &numEdges, &edges);
3363: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3364: DMComplexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3365: }
3366: }
3367: }
3368: }
3369: #if 0 /* Do not currently support holes */
3370: DMComplexCopyHoles(*dm, boundary);
3371: #endif
3372: FiniOutput_Triangle(&out);
3373: return(0);
3374: }
3375: #endif
3377: #ifdef PETSC_HAVE_TETGEN
3378: #include <tetgen.h>
3381: PetscErrorCode DMComplexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
3382: {
3383: MPI_Comm comm = ((PetscObject) boundary)->comm;
3384: DM_Complex *bd = (DM_Complex *) boundary->data;
3385: const PetscInt dim = 3;
3386: ::tetgenio in;
3387: ::tetgenio out;
3388: PetscInt vStart, vEnd, v, fStart, fEnd, f;
3389: PetscMPIInt rank;
3393: MPI_Comm_rank(comm, &rank);
3394: DMComplexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3395: in.numberofpoints = vEnd - vStart;
3396: if (in.numberofpoints > 0) {
3397: PetscScalar *array;
3399: in.pointlist = new double[in.numberofpoints*dim];
3400: in.pointmarkerlist = new int[in.numberofpoints];
3401: VecGetArray(bd->coordinates, &array);
3402: for(v = vStart; v < vEnd; ++v) {
3403: const PetscInt idx = v - vStart;
3404: PetscInt off, d;
3406: PetscSectionGetOffset(bd->coordSection, v, &off);
3407: for(d = 0; d < dim; ++d) {
3408: in.pointlist[idx*dim + d] = array[off+d];
3409: }
3410: DMComplexGetLabelValue(boundary, "marker", v, &in.pointmarkerlist[idx]);
3411: }
3412: VecRestoreArray(bd->coordinates, &array);
3413: }
3414: DMComplexGetHeightStratum(boundary, 0, &fStart, &fEnd);
3415: in.numberoffacets = fEnd - fStart;
3416: if (in.numberoffacets > 0) {
3417: in.facetlist = new tetgenio::facet[in.numberoffacets];
3418: in.facetmarkerlist = new int[in.numberoffacets];
3419: for(f = fStart; f < fEnd; ++f) {
3420: const PetscInt idx = f - fStart;
3421: PetscInt *points = PETSC_NULL, numPoints, p, numVertices = 0, v;
3423: in.facetlist[idx].numberofpolygons = 1;
3424: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
3425: in.facetlist[idx].numberofholes = 0;
3426: in.facetlist[idx].holelist = NULL;
3428: DMComplexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3429: for(p = 0; p < numPoints*2; p += 2) {
3430: const PetscInt point = points[p];
3431: if ((point >= vStart) && (point < vEnd)) {
3432: points[numVertices++] = point;
3433: }
3434: }
3436: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
3437: poly->numberofvertices = numVertices;
3438: poly->vertexlist = new int[poly->numberofvertices];
3439: for(v = 0; v < numVertices; ++v) {
3440: const PetscInt vIdx = points[v] - vStart;
3441: poly->vertexlist[v] = vIdx;
3442: }
3443: DMComplexGetLabelValue(boundary, "marker", f, &in.facetmarkerlist[idx]);
3444: }
3445: }
3446: if (!rank) {
3447: char args[32];
3449: /* Take away 'Q' for verbose output */
3450: PetscStrcpy(args, "pqezQ");
3451: ::tetrahedralize(args, &in, &out);
3452: }
3453: DMCreate(comm, dm);
3454: DMSetType(*dm, DMCOMPLEX);
3455: DMComplexSetDimension(*dm, dim);
3456: {
3457: DM_Complex *mesh = (DM_Complex *) (*dm)->data;
3458: const PetscInt numCorners = 4;
3459: const PetscInt numCells = out.numberoftetrahedra;
3460: const PetscInt numVertices = out.numberofpoints;
3461: int *cells = out.tetrahedronlist;
3462: double *meshCoords = out.pointlist;
3463: PetscInt coordSize, c;
3464: PetscScalar *coords;
3466: DMComplexSetChart(*dm, 0, numCells+numVertices);
3467: for(c = 0; c < numCells; ++c) {
3468: DMComplexSetConeSize(*dm, c, numCorners);
3469: }
3470: DMSetUp(*dm);
3471: for(c = 0; c < numCells; ++c) {
3472: /* Should be numCorners, but c89 sucks shit */
3473: PetscInt cone[4] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells, cells[c*numCorners+3]+numCells};
3475: DMComplexSetCone(*dm, c, cone);
3476: }
3477: DMComplexSymmetrize(*dm);
3478: DMComplexStratify(*dm);
3479: if (interpolate) {
3480: DM imesh;
3481: PetscInt *off;
3482: PetscInt firstFace = numCells+numVertices, numFaces = 0, face, f, firstEdge, numEdges = 0, edge, e;
3484: SETERRQ(comm, PETSC_ERR_SUP, "Interpolation is not yet implemented in 3D");
3485: /* TODO: Rewrite algorithm here to do all meets with neighboring cells and return counts */
3486: /* Count faces using algorithm from CreateNeighborCSR */
3487: DMComplexCreateNeighborCSR(*dm, PETSC_NULL, &off, PETSC_NULL);
3488: if (off) {
3489: numFaces = off[numCells]/2;
3490: /* Account for boundary faces: \sum_c 4 - neighbors = 4*numCells - totalNeighbors */
3491: numFaces += 4*numCells - off[numCells];
3492: }
3493: firstEdge = firstFace+numFaces;
3494: /* Create interpolated mesh */
3495: DMCreate(comm, &imesh);
3496: DMSetType(imesh, DMCOMPLEX);
3497: DMComplexSetDimension(imesh, dim);
3498: DMComplexSetChart(imesh, 0, numCells+numVertices+numEdges);
3499: for(c = 0; c < numCells; ++c) {
3500: DMComplexSetConeSize(imesh, c, numCorners);
3501: }
3502: for(f = firstFace; f < firstFace+numFaces; ++f) {
3503: DMComplexSetConeSize(imesh, f, 3);
3504: }
3505: for(e = firstEdge; e < firstEdge+numEdges; ++e) {
3506: DMComplexSetConeSize(imesh, e, 2);
3507: }
3508: DMSetUp(imesh);
3509: for(c = 0, face = firstFace; c < numCells; ++c) {
3510: const PetscInt *faces;
3511: PetscInt numFaces, faceSize, f;
3513: DMComplexGetFaces(*dm, c, &numFaces, &faceSize, &faces);
3514: if (faceSize != 2) SETERRQ1(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize);
3515: for(f = 0; f < numFaces; ++f) {
3516: PetscBool found = PETSC_FALSE;
3518: /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */
3519: for(e = firstEdge; e < edge; ++e) {
3520: const PetscInt *cone;
3522: DMComplexGetCone(imesh, e, &cone);
3523: if (((faces[f*faceSize+0] == cone[0]) && (faces[f*faceSize+1] == cone[1])) ||
3524: ((faces[f*faceSize+0] == cone[1]) && (faces[f*faceSize+1] == cone[0]))) {
3525: found = PETSC_TRUE;
3526: break;
3527: }
3528: }
3529: if (!found) {
3530: DMComplexSetCone(imesh, edge, &faces[f*faceSize]);
3531: ++edge;
3532: }
3533: DMComplexInsertCone(imesh, c, f, e);
3534: }
3535: }
3536: if (edge != firstEdge+numEdges) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges);
3537: PetscFree(off);
3538: DMComplexSymmetrize(imesh);
3539: DMComplexStratify(imesh);
3540: mesh = (DM_Complex *) (imesh)->data;
3541: for(c = 0; c < numCells; ++c) {
3542: const PetscInt *cone, *faces;
3543: PetscInt coneSize, coff, numFaces, faceSize, f;
3545: DMComplexGetConeSize(imesh, c, &coneSize);
3546: DMComplexGetCone(imesh, c, &cone);
3547: PetscSectionGetOffset(mesh->coneSection, c, &coff);
3548: DMComplexGetFaces(*dm, c, &numFaces, &faceSize, &faces);
3549: if (coneSize != numFaces) SETERRQ3(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numFaces);
3550: for(f = 0; f < numFaces; ++f) {
3551: const PetscInt *econe;
3552: PetscInt esize;
3554: DMComplexGetConeSize(imesh, cone[f], &esize);
3555: DMComplexGetCone(imesh, cone[f], &econe);
3556: if (esize != 2) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[f]);
3557: if ((faces[f*faceSize+0] == econe[0]) && (faces[f*faceSize+1] == econe[1])) {
3558: /* Correctly oriented */
3559: mesh->coneOrientations[coff+f] = 0;
3560: } else if ((faces[f*faceSize+0] == econe[1]) && (faces[f*faceSize+1] == econe[0])) {
3561: /* Start at index 1, and reverse orientation */
3562: mesh->coneOrientations[coff+f] = -(1+1);
3563: }
3564: }
3565: }
3566: DMDestroy(dm);
3567: *dm = imesh;
3568: }
3569: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3570: for(v = numCells; v < numCells+numVertices; ++v) {
3571: PetscSectionSetDof(mesh->coordSection, v, dim);
3572: }
3573: PetscSectionSetUp(mesh->coordSection);
3574: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3575: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3576: VecSetFromOptions(mesh->coordinates);
3577: VecGetArray(mesh->coordinates, &coords);
3578: for(v = 0; v < numVertices; ++v) {
3579: coords[v*dim+0] = meshCoords[v*dim+0];
3580: coords[v*dim+1] = meshCoords[v*dim+1];
3581: coords[v*dim+2] = meshCoords[v*dim+2];
3582: }
3583: VecRestoreArray(mesh->coordinates, &coords);
3584: for(v = 0; v < numVertices; ++v) {
3585: if (out.pointmarkerlist[v]) {
3586: DMComplexSetLabelValue(*dm, "marker", v+numCells, out.pointmarkerlist[v]);
3587: }
3588: }
3589: if (interpolate) {
3590: PetscInt e;
3592: for(e = 0; e < out.numberofedges; e++) {
3593: if (out.edgemarkerlist[e]) {
3594: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3595: const PetscInt *edges;
3596: PetscInt numEdges;
3598: DMComplexJoinPoints(*dm, 2, vertices, &numEdges, &edges);
3599: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3600: DMComplexSetLabelValue(*dm, "marker", edges[0], out.edgemarkerlist[e]);
3601: }
3602: }
3603: for(f = 0; f < out.numberoftrifaces; f++) {
3604: if (out.trifacemarkerlist[f]) {
3605: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3606: const PetscInt *faces;
3607: PetscInt numFaces;
3609: DMComplexJoinPoints(*dm, 3, vertices, &numFaces, &faces);
3610: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3611: DMComplexSetLabelValue(*dm, "marker", faces[0], out.trifacemarkerlist[f]);
3612: }
3613: }
3614: }
3615: }
3616: return(0);
3617: }
3621: PetscErrorCode DMComplexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
3622: {
3623: MPI_Comm comm = ((PetscObject) dm)->comm;
3624: DM_Complex *mesh = (DM_Complex *) dm->data;
3625: const PetscInt dim = 3;
3626: ::tetgenio in;
3627: ::tetgenio out;
3628: PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
3629: PetscMPIInt rank;
3633: MPI_Comm_rank(comm, &rank);
3634: DMComplexGetDepth(dm, &depth);
3635: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
3636: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
3637: in.numberofpoints = vEnd - vStart;
3638: if (in.numberofpoints > 0) {
3639: PetscScalar *array;
3641: in.pointlist = new double[in.numberofpoints*dim];
3642: in.pointmarkerlist = new int[in.numberofpoints];
3643: VecGetArray(mesh->coordinates, &array);
3644: for(v = vStart; v < vEnd; ++v) {
3645: const PetscInt idx = v - vStart;
3646: PetscInt off, d;
3648: PetscSectionGetOffset(mesh->coordSection, v, &off);
3649: for(d = 0; d < dim; ++d) {
3650: in.pointlist[idx*dim + d] = array[off+d];
3651: }
3652: DMComplexGetLabelValue(dm, "marker", v, &in.pointmarkerlist[idx]);
3653: }
3654: VecRestoreArray(mesh->coordinates, &array);
3655: }
3656: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
3657: in.numberofcorners = 4;
3658: in.numberoftetrahedra = cEnd - cStart;
3659: in.tetrahedronvolumelist = (double *) maxVolumes;
3660: if (in.numberoftetrahedra > 0) {
3661: in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
3662: for(c = cStart; c < cEnd; ++c) {
3663: const PetscInt idx = c - cStart;
3664: PetscInt *closure = PETSC_NULL;
3665: PetscInt closureSize;
3667: DMComplexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
3668: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
3669: for(v = 0; v < 4; ++v) {
3670: in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
3671: }
3672: }
3673: }
3674: // TODO: Put in boundary faces with markers
3675: if (!rank) {
3676: char args[32];
3678: /* Take away 'Q' for verbose output */
3679: //PetscStrcpy(args, "qezQra");
3680: PetscStrcpy(args, "qezraVVVV");
3681: ::tetrahedralize(args, &in, &out);
3682: }
3683: in.tetrahedronvolumelist = NULL;
3685: DMCreate(comm, dmRefined);
3686: DMSetType(*dmRefined, DMCOMPLEX);
3687: DMComplexSetDimension(*dmRefined, dim);
3688: {
3689: DM_Complex *mesh = (DM_Complex *) (*dmRefined)->data;
3690: const PetscInt numCorners = 4;
3691: const PetscInt numCells = out.numberoftetrahedra;
3692: const PetscInt numVertices = out.numberofpoints;
3693: int *cells = out.tetrahedronlist;
3694: double *meshCoords = out.pointlist;
3695: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
3696: PetscInt coordSize, c, e;
3697: PetscScalar *coords;
3699: DMComplexSetChart(*dmRefined, 0, numCells+numVertices);
3700: for(c = 0; c < numCells; ++c) {
3701: DMComplexSetConeSize(*dmRefined, c, numCorners);
3702: }
3703: DMSetUp(*dmRefined);
3704: for(c = 0; c < numCells; ++c) {
3705: /* Should be numCorners, but c89 sucks shit */
3706: PetscInt cone[4] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells, cells[c*numCorners+3]+numCells};
3708: DMComplexSetCone(*dmRefined, c, cone);
3709: }
3710: DMComplexSymmetrize(*dmRefined);
3711: DMComplexStratify(*dmRefined);
3713: if (interpolate) {
3714: DM imesh;
3715: PetscInt *off;
3716: PetscInt firstEdge = numCells+numVertices, numEdges, edge, e;
3718: SETERRQ(comm, PETSC_ERR_SUP, "Interpolation is not yet implemented in 3D");
3719: /* Count edges using algorithm from CreateNeighborCSR */
3720: DMComplexCreateNeighborCSR(*dmRefined, PETSC_NULL, &off, PETSC_NULL);
3721: if (off) {
3722: numEdges = off[numCells]/2;
3723: /* Account for boundary edges: \sum_c 3 - neighbors = 3*numCells - totalNeighbors */
3724: numEdges += 3*numCells - off[numCells];
3725: }
3726: /* Create interpolated mesh */
3727: DMCreate(comm, &imesh);
3728: DMSetType(imesh, DMCOMPLEX);
3729: DMComplexSetDimension(imesh, dim);
3730: DMComplexSetChart(imesh, 0, numCells+numVertices+numEdges);
3731: for(c = 0; c < numCells; ++c) {
3732: DMComplexSetConeSize(imesh, c, numCorners);
3733: }
3734: for(e = firstEdge; e < firstEdge+numEdges; ++e) {
3735: DMComplexSetConeSize(imesh, e, 2);
3736: }
3737: DMSetUp(imesh);
3738: for(c = 0, edge = firstEdge; c < numCells; ++c) {
3739: const PetscInt *faces;
3740: PetscInt numFaces, faceSize, f;
3742: DMComplexGetFaces(*dmRefined, c, &numFaces, &faceSize, &faces);
3743: if (faceSize != 2) SETERRQ1(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize);
3744: for(f = 0; f < numFaces; ++f) {
3745: PetscBool found = PETSC_FALSE;
3747: /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */
3748: for(e = firstEdge; e < edge; ++e) {
3749: const PetscInt *cone;
3751: DMComplexGetCone(imesh, e, &cone);
3752: if (((faces[f*faceSize+0] == cone[0]) && (faces[f*faceSize+1] == cone[1])) ||
3753: ((faces[f*faceSize+0] == cone[1]) && (faces[f*faceSize+1] == cone[0]))) {
3754: found = PETSC_TRUE;
3755: break;
3756: }
3757: }
3758: if (!found) {
3759: DMComplexSetCone(imesh, edge, &faces[f*faceSize]);
3760: ++edge;
3761: }
3762: DMComplexInsertCone(imesh, c, f, e);
3763: }
3764: }
3765: if (edge != firstEdge+numEdges) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges);
3766: PetscFree(off);
3767: DMComplexSymmetrize(imesh);
3768: DMComplexStratify(imesh);
3769: mesh = (DM_Complex *) (imesh)->data;
3770: for(c = 0; c < numCells; ++c) {
3771: const PetscInt *cone, *faces;
3772: PetscInt coneSize, coff, numFaces, faceSize, f;
3774: DMComplexGetConeSize(imesh, c, &coneSize);
3775: DMComplexGetCone(imesh, c, &cone);
3776: PetscSectionGetOffset(mesh->coneSection, c, &coff);
3777: DMComplexGetFaces(*dmRefined, c, &numFaces, &faceSize, &faces);
3778: if (coneSize != numFaces) SETERRQ3(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numFaces);
3779: for(f = 0; f < numFaces; ++f) {
3780: const PetscInt *econe;
3781: PetscInt esize;
3783: DMComplexGetConeSize(imesh, cone[f], &esize);
3784: DMComplexGetCone(imesh, cone[f], &econe);
3785: if (esize != 2) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[f]);
3786: if ((faces[f*faceSize+0] == econe[0]) && (faces[f*faceSize+1] == econe[1])) {
3787: /* Correctly oriented */
3788: mesh->coneOrientations[coff+f] = 0;
3789: } else if ((faces[f*faceSize+0] == econe[1]) && (faces[f*faceSize+1] == econe[0])) {
3790: /* Start at index 1, and reverse orientation */
3791: mesh->coneOrientations[coff+f] = -(1+1);
3792: }
3793: }
3794: }
3795: DMDestroy(dmRefined);
3796: *dmRefined = imesh;
3797: }
3798: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
3799: for(v = numCells; v < numCells+numVertices; ++v) {
3800: PetscSectionSetDof(mesh->coordSection, v, dim);
3801: }
3802: PetscSectionSetUp(mesh->coordSection);
3803: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
3804: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
3805: VecSetFromOptions(mesh->coordinates);
3806: VecGetArray(mesh->coordinates, &coords);
3807: for(v = 0; v < numVertices; ++v) {
3808: coords[v*dim+0] = meshCoords[v*dim+0];
3809: coords[v*dim+1] = meshCoords[v*dim+1];
3810: }
3811: VecRestoreArray(mesh->coordinates, &coords);
3812: for(v = 0; v < numVertices; ++v) {
3813: if (out.pointmarkerlist[v]) {
3814: DMComplexSetLabelValue(*dmRefined, "marker", v+numCells, out.pointmarkerlist[v]);
3815: }
3816: }
3817: if (interpolate) {
3818: PetscInt f;
3820: for(e = 0; e < out.numberofedges; e++) {
3821: if (out.edgemarkerlist[e]) {
3822: const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
3823: const PetscInt *edges;
3824: PetscInt numEdges;
3826: DMComplexJoinPoints(*dmRefined, 2, vertices, &numEdges, &edges);
3827: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
3828: DMComplexSetLabelValue(*dmRefined, "marker", edges[0], out.edgemarkerlist[e]);
3829: }
3830: }
3831: for(f = 0; f < out.numberoftrifaces; f++) {
3832: if (out.trifacemarkerlist[f]) {
3833: const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
3834: const PetscInt *faces;
3835: PetscInt numFaces;
3837: DMComplexJoinPoints(*dmRefined, 3, vertices, &numFaces, &faces);
3838: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
3839: DMComplexSetLabelValue(*dmRefined, "marker", faces[0], out.trifacemarkerlist[f]);
3840: }
3841: }
3842: }
3843: }
3844: return(0);
3845: }
3846: #endif
3848: #ifdef PETSC_HAVE_CTETGEN
3849: #include "ctetgen.h"
3853: PetscErrorCode DMComplexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
3854: {
3855: MPI_Comm comm = ((PetscObject) boundary)->comm;
3856: DM_Complex *bd = (DM_Complex *) boundary->data;
3857: const PetscInt dim = 3;
3858: PLC *in, *out;
3859: PetscInt verbose = 0, vStart, vEnd, v, fStart, fEnd, f;
3860: PetscMPIInt rank;
3864: PetscOptionsGetInt(((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, PETSC_NULL);
3865: MPI_Comm_rank(comm, &rank);
3866: DMComplexGetDepthStratum(boundary, 0, &vStart, &vEnd);
3867: PLCCreate(&in);
3868: PLCCreate(&out);
3869: in->numberofpoints = vEnd - vStart;
3870: if (in->numberofpoints > 0) {
3871: PetscScalar *array;
3873: PetscMalloc(in->numberofpoints*dim * sizeof(PetscReal), &in->pointlist);
3874: PetscMalloc(in->numberofpoints * sizeof(int), &in->pointmarkerlist);
3875: VecGetArray(bd->coordinates, &array);
3876: for(v = vStart; v < vEnd; ++v) {
3877: const PetscInt idx = v - vStart;
3878: PetscInt off, d, m;
3880: PetscSectionGetOffset(bd->coordSection, v, &off);
3881: for(d = 0; d < dim; ++d) {
3882: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
3883: }
3884: DMComplexGetLabelValue(boundary, "marker", v, &m);
3885: in->pointmarkerlist[idx] = (int) m;
3886: }
3887: VecRestoreArray(bd->coordinates, &array);
3888: }
3889: DMComplexGetHeightStratum(boundary, 0, &fStart, &fEnd);
3890: in->numberoffacets = fEnd - fStart;
3891: if (in->numberoffacets > 0) {
3892: PetscMalloc(in->numberoffacets * sizeof(facet), &in->facetlist);
3893: PetscMalloc(in->numberoffacets * sizeof(int), &in->facetmarkerlist);
3894: for(f = fStart; f < fEnd; ++f) {
3895: const PetscInt idx = f - fStart;
3896: PetscInt *points = PETSC_NULL, numPoints, p, numVertices = 0, v, m;
3897: polygon *poly;
3899: in->facetlist[idx].numberofpolygons = 1;
3900: PetscMalloc(in->facetlist[idx].numberofpolygons * sizeof(polygon), &in->facetlist[idx].polygonlist);
3901: in->facetlist[idx].numberofholes = 0;
3902: in->facetlist[idx].holelist = PETSC_NULL;
3904: DMComplexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
3905: for(p = 0; p < numPoints*2; p += 2) {
3906: const PetscInt point = points[p];
3907: if ((point >= vStart) && (point < vEnd)) {
3908: points[numVertices++] = point;
3909: }
3910: }
3912: poly = in->facetlist[idx].polygonlist;
3913: poly->numberofvertices = numVertices;
3914: PetscMalloc(poly->numberofvertices * sizeof(int), &poly->vertexlist);
3915: for(v = 0; v < numVertices; ++v) {
3916: const PetscInt vIdx = points[v] - vStart;
3917: poly->vertexlist[v] = vIdx;
3918: }
3919: DMComplexGetLabelValue(boundary, "marker", f, &m);
3920: in->facetmarkerlist[idx] = (int) m;
3921: }
3922: }
3923: if (!rank) {
3924: TetGenOpts t;
3926: TetGenOptsInitialize(&t);
3927: t.in = boundary; /* Should go away */
3928: t.plc = 1;
3929: t.quality = 1;
3930: t.edgesout = 1;
3931: t.zeroindex = 1;
3932: t.quiet = 1;
3933: t.verbose = verbose;
3934: TetGenCheckOpts(&t);
3935: TetGenTetrahedralize(&t, in, out);
3936: }
3937: DMCreate(comm, dm);
3938: DMSetType(*dm, DMCOMPLEX);
3939: DMComplexSetDimension(*dm, dim);
3940: {
3941: DM_Complex *mesh = (DM_Complex *) (*dm)->data;
3942: const PetscInt numCorners = 4;
3943: const PetscInt numCells = out->numberoftetrahedra;
3944: const PetscInt numVertices = out->numberofpoints;
3945: int *cells = out->tetrahedronlist;
3946: PetscReal *meshCoords = out->pointlist;
3947: PetscInt coordSize, c;
3948: PetscScalar *coords;
3950: DMComplexSetChart(*dm, 0, numCells+numVertices);
3951: for(c = 0; c < numCells; ++c) {
3952: DMComplexSetConeSize(*dm, c, numCorners);
3953: }
3954: DMSetUp(*dm);
3955: for(c = 0; c < numCells; ++c) {
3956: /* Should be numCorners, but c89 sucks shit */
3957: PetscInt cone[4] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells, cells[c*numCorners+3]+numCells};
3959: DMComplexSetCone(*dm, c, cone);
3960: }
3961: DMComplexSymmetrize(*dm);
3962: DMComplexStratify(*dm);
3963: if (interpolate) {
3964: DM imesh;
3965: PetscInt *off;
3966: PetscInt firstFace = numCells+numVertices, numFaces = 0, f, firstEdge, numEdges = 0, edge, e;
3968: SETERRQ(comm, PETSC_ERR_SUP, "Interpolation is not yet implemented in 3D");
3969: /* TODO: Rewrite algorithm here to do all meets with neighboring cells and return counts */
3970: /* Count faces using algorithm from CreateNeighborCSR */
3971: DMComplexCreateNeighborCSR(*dm, PETSC_NULL, &off, PETSC_NULL);
3972: if (off) {
3973: numFaces = off[numCells]/2;
3974: /* Account for boundary faces: \sum_c 4 - neighbors = 4*numCells - totalNeighbors */
3975: numFaces += 4*numCells - off[numCells];
3976: }
3977: firstEdge = firstFace+numFaces;
3978: /* Create interpolated mesh */
3979: DMCreate(comm, &imesh);
3980: DMSetType(imesh, DMCOMPLEX);
3981: DMComplexSetDimension(imesh, dim);
3982: DMComplexSetChart(imesh, 0, numCells+numVertices+numEdges);
3983: for(c = 0; c < numCells; ++c) {
3984: DMComplexSetConeSize(imesh, c, numCorners);
3985: }
3986: for(f = firstFace; f < firstFace+numFaces; ++f) {
3987: DMComplexSetConeSize(imesh, f, 3);
3988: }
3989: for(e = firstEdge; e < firstEdge+numEdges; ++e) {
3990: DMComplexSetConeSize(imesh, e, 2);
3991: }
3992: DMSetUp(imesh);
3993: for(c = 0; c < numCells; ++c) {
3994: const PetscInt *faces;
3995: PetscInt numFaces, faceSize, f;
3997: DMComplexGetFaces(*dm, c, &numFaces, &faceSize, &faces);
3998: if (faceSize != 2) SETERRQ1(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize);
3999: for(f = 0; f < numFaces; ++f) {
4000: PetscBool found = PETSC_FALSE;
4002: /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */
4003: for(e = firstEdge; e < edge; ++e) {
4004: const PetscInt *cone;
4006: DMComplexGetCone(imesh, e, &cone);
4007: if (((faces[f*faceSize+0] == cone[0]) && (faces[f*faceSize+1] == cone[1])) ||
4008: ((faces[f*faceSize+0] == cone[1]) && (faces[f*faceSize+1] == cone[0]))) {
4009: found = PETSC_TRUE;
4010: break;
4011: }
4012: }
4013: if (!found) {
4014: DMComplexSetCone(imesh, edge, &faces[f*faceSize]);
4015: ++edge;
4016: }
4017: DMComplexInsertCone(imesh, c, f, e);
4018: }
4019: }
4020: if (edge != firstEdge+numEdges) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges);
4021: PetscFree(off);
4022: DMComplexSymmetrize(imesh);
4023: DMComplexStratify(imesh);
4024: mesh = (DM_Complex *) (imesh)->data;
4025: for(c = 0; c < numCells; ++c) {
4026: const PetscInt *cone, *faces;
4027: PetscInt coneSize, coff, numFaces, faceSize, f;
4029: DMComplexGetConeSize(imesh, c, &coneSize);
4030: DMComplexGetCone(imesh, c, &cone);
4031: PetscSectionGetOffset(mesh->coneSection, c, &coff);
4032: DMComplexGetFaces(*dm, c, &numFaces, &faceSize, &faces);
4033: if (coneSize != numFaces) SETERRQ3(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numFaces);
4034: for(f = 0; f < numFaces; ++f) {
4035: const PetscInt *econe;
4036: PetscInt esize;
4038: DMComplexGetConeSize(imesh, cone[f], &esize);
4039: DMComplexGetCone(imesh, cone[f], &econe);
4040: if (esize != 2) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[f]);
4041: if ((faces[f*faceSize+0] == econe[0]) && (faces[f*faceSize+1] == econe[1])) {
4042: /* Correctly oriented */
4043: mesh->coneOrientations[coff+f] = 0;
4044: } else if ((faces[f*faceSize+0] == econe[1]) && (faces[f*faceSize+1] == econe[0])) {
4045: /* Start at index 1, and reverse orientation */
4046: mesh->coneOrientations[coff+f] = -(1+1);
4047: }
4048: }
4049: }
4050: DMDestroy(dm);
4051: *dm = imesh;
4052: }
4053: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
4054: for(v = numCells; v < numCells+numVertices; ++v) {
4055: PetscSectionSetDof(mesh->coordSection, v, dim);
4056: }
4057: PetscSectionSetUp(mesh->coordSection);
4058: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
4059: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
4060: VecSetFromOptions(mesh->coordinates);
4061: VecGetArray(mesh->coordinates, &coords);
4062: for(v = 0; v < numVertices; ++v) {
4063: coords[v*dim+0] = meshCoords[v*dim+0];
4064: coords[v*dim+1] = meshCoords[v*dim+1];
4065: coords[v*dim+2] = meshCoords[v*dim+2];
4066: }
4067: VecRestoreArray(mesh->coordinates, &coords);
4068: for(v = 0; v < numVertices; ++v) {
4069: if (out->pointmarkerlist[v]) {
4070: DMComplexSetLabelValue(*dm, "marker", v+numCells, out->pointmarkerlist[v]);
4071: }
4072: }
4073: if (interpolate) {
4074: PetscInt e;
4076: for(e = 0; e < out->numberofedges; e++) {
4077: if (out->edgemarkerlist[e]) {
4078: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
4079: const PetscInt *edges;
4080: PetscInt numEdges;
4082: DMComplexJoinPoints(*dm, 2, vertices, &numEdges, &edges);
4083: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
4084: DMComplexSetLabelValue(*dm, "marker", edges[0], out->edgemarkerlist[e]);
4085: }
4086: }
4087: for(f = 0; f < out->numberoftrifaces; f++) {
4088: if (out->trifacemarkerlist[f]) {
4089: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
4090: const PetscInt *faces;
4091: PetscInt numFaces;
4093: DMComplexJoinPoints(*dm, 3, vertices, &numFaces, &faces);
4094: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
4095: DMComplexSetLabelValue(*dm, "marker", faces[0], out->trifacemarkerlist[f]);
4096: }
4097: }
4098: }
4099: }
4101: PLCDestroy(&in);
4102: PLCDestroy(&out);
4103: return(0);
4104: }
4108: PetscErrorCode DMComplexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
4109: {
4110: MPI_Comm comm = ((PetscObject) dm)->comm;
4111: DM_Complex *mesh = (DM_Complex *) dm->data;
4112: const PetscInt dim = 3;
4113: PLC *in, *out;
4114: PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal;
4115: PetscMPIInt rank;
4119: PetscOptionsGetInt(((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, PETSC_NULL);
4120: MPI_Comm_rank(comm, &rank);
4121: DMComplexGetDepth(dm, &depth);
4122: MPI_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);
4123: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
4124: PLCCreate(&in);
4125: PLCCreate(&out);
4126: in->numberofpoints = vEnd - vStart;
4127: if (in->numberofpoints > 0) {
4128: PetscScalar *array;
4130: PetscMalloc(in->numberofpoints*dim * sizeof(PetscReal), &in->pointlist);
4131: PetscMalloc(in->numberofpoints * sizeof(int), &in->pointmarkerlist);
4132: VecGetArray(mesh->coordinates, &array);
4133: for(v = vStart; v < vEnd; ++v) {
4134: const PetscInt idx = v - vStart;
4135: PetscInt off, d, m;
4137: PetscSectionGetOffset(mesh->coordSection, v, &off);
4138: for(d = 0; d < dim; ++d) {
4139: in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
4140: }
4141: DMComplexGetLabelValue(dm, "marker", v, &m);
4142: in->pointmarkerlist[idx] = (int) m;
4143: }
4144: VecRestoreArray(mesh->coordinates, &array);
4145: }
4146: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
4147: in->numberofcorners = 4;
4148: in->numberoftetrahedra = cEnd - cStart;
4149: in->tetrahedronvolumelist = maxVolumes;
4150: if (in->numberoftetrahedra > 0) {
4151: PetscMalloc(in->numberoftetrahedra*in->numberofcorners * sizeof(int), &in->tetrahedronlist);
4152: for(c = cStart; c < cEnd; ++c) {
4153: const PetscInt idx = c - cStart;
4154: PetscInt *closure = PETSC_NULL;
4155: PetscInt closureSize;
4157: DMComplexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
4158: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
4159: for(v = 0; v < 4; ++v) {
4160: in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
4161: }
4162: }
4163: }
4164: if (!rank) {
4165: TetGenOpts t;
4167: TetGenOptsInitialize(&t);
4168: t.in = dm; /* Should go away */
4169: t.refine = 1;
4170: t.varvolume = 1;
4171: t.quality = 1;
4172: t.edgesout = 1;
4173: t.zeroindex = 1;
4174: t.quiet = 1;
4175: t.verbose = verbose; /* Change this */
4176: TetGenCheckOpts(&t);
4177: TetGenTetrahedralize(&t, in, out);
4178: }
4180: DMCreate(comm, dmRefined);
4181: DMSetType(*dmRefined, DMCOMPLEX);
4182: DMComplexSetDimension(*dmRefined, dim);
4183: {
4184: DM_Complex *mesh = (DM_Complex *) (*dmRefined)->data;
4185: const PetscInt numCorners = 4;
4186: const PetscInt numCells = out->numberoftetrahedra;
4187: const PetscInt numVertices = out->numberofpoints;
4188: int *cells = out->tetrahedronlist;
4189: PetscReal *meshCoords = out->pointlist;
4190: PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE;
4191: PetscInt coordSize, c, e;
4192: PetscScalar *coords;
4194: DMComplexSetChart(*dmRefined, 0, numCells+numVertices);
4195: for(c = 0; c < numCells; ++c) {
4196: DMComplexSetConeSize(*dmRefined, c, numCorners);
4197: }
4198: DMSetUp(*dmRefined);
4199: for(c = 0; c < numCells; ++c) {
4200: /* Should be numCorners, but c89 sucks shit */
4201: PetscInt cone[4] = {cells[c*numCorners+0]+numCells, cells[c*numCorners+1]+numCells, cells[c*numCorners+2]+numCells, cells[c*numCorners+3]+numCells};
4203: DMComplexSetCone(*dmRefined, c, cone);
4204: }
4205: DMComplexSymmetrize(*dmRefined);
4206: DMComplexStratify(*dmRefined);
4208: if (interpolate) {
4209: DM imesh;
4210: PetscInt *off;
4211: PetscInt firstEdge = numCells+numVertices, numEdges, edge, e;
4213: SETERRQ(comm, PETSC_ERR_SUP, "Interpolation is not yet implemented in 3D");
4214: /* Count edges using algorithm from CreateNeighborCSR */
4215: DMComplexCreateNeighborCSR(*dmRefined, PETSC_NULL, &off, PETSC_NULL);
4216: if (off) {
4217: numEdges = off[numCells]/2;
4218: /* Account for boundary edges: \sum_c 3 - neighbors = 3*numCells - totalNeighbors */
4219: numEdges += 3*numCells - off[numCells];
4220: }
4221: /* Create interpolated mesh */
4222: DMCreate(comm, &imesh);
4223: DMSetType(imesh, DMCOMPLEX);
4224: DMComplexSetDimension(imesh, dim);
4225: DMComplexSetChart(imesh, 0, numCells+numVertices+numEdges);
4226: for(c = 0; c < numCells; ++c) {
4227: DMComplexSetConeSize(imesh, c, numCorners);
4228: }
4229: for(e = firstEdge; e < firstEdge+numEdges; ++e) {
4230: DMComplexSetConeSize(imesh, e, 2);
4231: }
4232: DMSetUp(imesh);
4233: for(c = 0, edge = firstEdge; c < numCells; ++c) {
4234: const PetscInt *faces;
4235: PetscInt numFaces, faceSize, f;
4237: DMComplexGetFaces(*dmRefined, c, &numFaces, &faceSize, &faces);
4238: if (faceSize != 2) SETERRQ1(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Triangles cannot have face of size %D", faceSize);
4239: for(f = 0; f < numFaces; ++f) {
4240: PetscBool found = PETSC_FALSE;
4242: /* TODO Need join of vertices to check for existence of edges, which needs support (could set edge support), so just brute force for now */
4243: for(e = firstEdge; e < edge; ++e) {
4244: const PetscInt *cone;
4246: DMComplexGetCone(imesh, e, &cone);
4247: if (((faces[f*faceSize+0] == cone[0]) && (faces[f*faceSize+1] == cone[1])) ||
4248: ((faces[f*faceSize+0] == cone[1]) && (faces[f*faceSize+1] == cone[0]))) {
4249: found = PETSC_TRUE;
4250: break;
4251: }
4252: }
4253: if (!found) {
4254: DMComplexSetCone(imesh, edge, &faces[f*faceSize]);
4255: ++edge;
4256: }
4257: DMComplexInsertCone(imesh, c, f, e);
4258: }
4259: }
4260: if (edge != firstEdge+numEdges) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D should be %D", edge-firstEdge, numEdges);
4261: PetscFree(off);
4262: DMComplexSymmetrize(imesh);
4263: DMComplexStratify(imesh);
4264: mesh = (DM_Complex *) (imesh)->data;
4265: for(c = 0; c < numCells; ++c) {
4266: const PetscInt *cone, *faces;
4267: PetscInt coneSize, coff, numFaces, faceSize, f;
4269: DMComplexGetConeSize(imesh, c, &coneSize);
4270: DMComplexGetCone(imesh, c, &cone);
4271: PetscSectionGetOffset(mesh->coneSection, c, &coff);
4272: DMComplexGetFaces(*dmRefined, c, &numFaces, &faceSize, &faces);
4273: if (coneSize != numFaces) SETERRQ3(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edges %D for cell %D should be %D", coneSize, c, numFaces);
4274: for(f = 0; f < numFaces; ++f) {
4275: const PetscInt *econe;
4276: PetscInt esize;
4278: DMComplexGetConeSize(imesh, cone[f], &esize);
4279: DMComplexGetCone(imesh, cone[f], &econe);
4280: if (esize != 2) SETERRQ2(((PetscObject) imesh)->comm, PETSC_ERR_PLIB, "Invalid number of edge endpoints %D for edge %D should be 2", esize, cone[f]);
4281: if ((faces[f*faceSize+0] == econe[0]) && (faces[f*faceSize+1] == econe[1])) {
4282: /* Correctly oriented */
4283: mesh->coneOrientations[coff+f] = 0;
4284: } else if ((faces[f*faceSize+0] == econe[1]) && (faces[f*faceSize+1] == econe[0])) {
4285: /* Start at index 1, and reverse orientation */
4286: mesh->coneOrientations[coff+f] = -(1+1);
4287: }
4288: }
4289: }
4290: DMDestroy(dmRefined);
4291: *dmRefined = imesh;
4292: }
4293: PetscSectionSetChart(mesh->coordSection, numCells, numCells + numVertices);
4294: for(v = numCells; v < numCells+numVertices; ++v) {
4295: PetscSectionSetDof(mesh->coordSection, v, dim);
4296: }
4297: PetscSectionSetUp(mesh->coordSection);
4298: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
4299: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
4300: VecSetFromOptions(mesh->coordinates);
4301: VecGetArray(mesh->coordinates, &coords);
4302: for(v = 0; v < numVertices; ++v) {
4303: coords[v*dim+0] = meshCoords[v*dim+0];
4304: coords[v*dim+1] = meshCoords[v*dim+1];
4305: coords[v*dim+2] = meshCoords[v*dim+2];
4306: }
4307: VecRestoreArray(mesh->coordinates, &coords);
4308: for(v = 0; v < numVertices; ++v) {
4309: if (out->pointmarkerlist[v]) {
4310: DMComplexSetLabelValue(*dmRefined, "marker", v+numCells, out->pointmarkerlist[v]);
4311: }
4312: }
4313: if (interpolate) {
4314: PetscInt f;
4316: for(e = 0; e < out->numberofedges; e++) {
4317: if (out->edgemarkerlist[e]) {
4318: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
4319: const PetscInt *edges;
4320: PetscInt numEdges;
4322: DMComplexJoinPoints(*dmRefined, 2, vertices, &numEdges, &edges);
4323: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
4324: DMComplexSetLabelValue(*dmRefined, "marker", edges[0], out->edgemarkerlist[e]);
4325: }
4326: }
4327: for(f = 0; f < out->numberoftrifaces; f++) {
4328: if (out->trifacemarkerlist[f]) {
4329: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
4330: const PetscInt *faces;
4331: PetscInt numFaces;
4333: DMComplexJoinPoints(*dmRefined, 3, vertices, &numFaces, &faces);
4334: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
4335: DMComplexSetLabelValue(*dmRefined, "marker", faces[0], out->trifacemarkerlist[f]);
4336: }
4337: }
4338: }
4339: }
4340: PLCDestroy(&in);
4341: PLCDestroy(&out);
4342: return(0);
4343: }
4344: #endif
4348: /*@C
4349: DMComplexGenerate - Generates a mesh.
4351: Not Collective
4353: Input Parameters:
4354: + boundary - The DMComplex boundary object
4355: . name - The mesh generation package name
4356: - interpolate - Flag to create intermediate mesh elements
4358: Output Parameter:
4359: . mesh - The DMComplex object
4361: Level: intermediate
4363: .keywords: mesh, elements
4364: .seealso: DMComplexCreate(), DMRefine()
4365: @*/
4366: PetscErrorCode DMComplexGenerate(DM boundary, const char name[], PetscBool interpolate, DM *mesh)
4367: {
4368: PetscInt dim;
4369: char genname[1024];
4370: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
4376: DMComplexGetDimension(boundary, &dim);
4377: PetscOptionsGetString(((PetscObject) boundary)->prefix, "-dm_complex_generator", genname, 1024, &flg);
4378: if (flg) {name = genname;}
4379: if (name) {
4380: PetscStrcmp(name, "triangle", &isTriangle);
4381: PetscStrcmp(name, "tetgen", &isTetgen);
4382: PetscStrcmp(name, "ctetgen", &isCTetgen);
4383: }
4384: switch(dim) {
4385: case 1:
4386: if (!name || isTriangle) {
4387: #ifdef PETSC_HAVE_TRIANGLE
4388: DMComplexGenerate_Triangle(boundary, interpolate, mesh);
4389: #else
4390: SETERRQ(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Mesh generation needs external package support.\nPlease reconfigure with --download-triangle.");
4391: #endif
4392: } else SETERRQ1(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
4393: break;
4394: case 2:
4395: if (!name || isCTetgen) {
4396: #ifdef PETSC_HAVE_CTETGEN
4397: DMComplexGenerate_CTetgen(boundary, interpolate, mesh);
4398: #else
4399: SETERRQ(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --download-ctetgen.");
4400: #endif
4401: } else if (isTetgen) {
4402: #ifdef PETSC_HAVE_TETGEN
4403: DMComplexGenerate_Tetgen(boundary, interpolate, mesh);
4404: #else
4405: SETERRQ(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
4406: #endif
4407: } else SETERRQ1(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
4408: break;
4409: default:
4410: SETERRQ1(((PetscObject) boundary)->comm, PETSC_ERR_SUP, "Mesh generation for a dimension %d boundary is not supported.", dim);
4411: }
4412: return(0);
4413: }
4417: PetscErrorCode DMComplexSetRefinementLimit(DM dm, PetscReal refinementLimit)
4418: {
4419: DM_Complex *mesh = (DM_Complex *) dm->data;
4423: mesh->refinementLimit = refinementLimit;
4424: return(0);
4425: }
4429: PetscErrorCode DMComplexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
4430: {
4431: DM_Complex *mesh = (DM_Complex *) dm->data;
4437: /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
4438: *refinementLimit = mesh->refinementLimit;
4439: return(0);
4440: }
4444: PetscErrorCode DMRefine_Complex(DM dm, MPI_Comm comm, DM *dmRefined)
4445: {
4446: PetscReal refinementLimit;
4447: PetscInt dim, cStart, cEnd;
4448: char genname[1024], *name = PETSC_NULL;
4449: PetscBool isTriangle = PETSC_FALSE, isTetgen = PETSC_FALSE, isCTetgen = PETSC_FALSE, flg;
4453: DMComplexGetRefinementLimit(dm, &refinementLimit);
4454: if (refinementLimit == 0.0) return(0);
4455: DMComplexGetDimension(dm, &dim);
4456: DMComplexGetHeightStratum(dm, 0, &cStart, &cEnd);
4457: PetscOptionsGetString(((PetscObject) dm)->prefix, "-dm_complex_generator", genname, 1024, &flg);
4458: if (flg) {name = genname;}
4459: if (name) {
4460: PetscStrcmp(name, "triangle", &isTriangle);
4461: PetscStrcmp(name, "tetgen", &isTetgen);
4462: PetscStrcmp(name, "ctetgen", &isCTetgen);
4463: }
4464: switch(dim) {
4465: case 2:
4466: if (!name || isTriangle) {
4467: #ifdef PETSC_HAVE_TRIANGLE
4468: double *maxVolumes;
4469: PetscInt c;
4471: PetscMalloc((cEnd - cStart) * sizeof(double), &maxVolumes);
4472: for(c = 0; c < cEnd-cStart; ++c) {
4473: maxVolumes[c] = refinementLimit;
4474: }
4475: DMComplexRefine_Triangle(dm, maxVolumes, dmRefined);
4476: #else
4477: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Mesh refinement needs external package support.\nPlease reconfigure with --download-triangle.");
4478: #endif
4479: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Unknown 2D mesh generation package %s", name);
4480: break;
4481: case 3:
4482: if (!name || isCTetgen) {
4483: #ifdef PETSC_HAVE_CTETGEN
4484: PetscReal *maxVolumes;
4485: PetscInt c;
4487: PetscMalloc((cEnd - cStart) * sizeof(PetscReal), &maxVolumes);
4488: for(c = 0; c < cEnd-cStart; ++c) {
4489: maxVolumes[c] = refinementLimit;
4490: }
4491: DMComplexRefine_CTetgen(dm, maxVolumes, dmRefined);
4492: #else
4493: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "CTetgen needs external package support.\nPlease reconfigure with --with-ctetgen.");
4494: #endif
4495: } else if (isTetgen) {
4496: #ifdef PETSC_HAVE_TETGEN
4497: double *maxVolumes;
4498: PetscInt c;
4500: PetscMalloc((cEnd - cStart) * sizeof(double), &maxVolumes);
4501: for(c = 0; c < cEnd-cStart; ++c) {
4502: maxVolumes[c] = refinementLimit;
4503: }
4504: DMComplexRefine_Tetgen(dm, maxVolumes, dmRefined);
4505: #else
4506: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Tetgen needs external package support.\nPlease reconfigure with --with-c-language=cxx --download-tetgen.");
4507: #endif
4508: } else SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Unknown 3D mesh generation package %s", name);
4509: break;
4510: default:
4511: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Mesh refinement in dimension %d is not supported.", dim);
4512: }
4513: return(0);
4514: }
4518: PetscErrorCode DMComplexGetDepth(DM dm, PetscInt *depth) {
4519: PetscInt d;
4525: DMComplexGetLabelSize(dm, "depth", &d);
4526: *depth = d-1;
4527: return(0);
4528: }
4532: PetscErrorCode DMComplexGetDepthStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) {
4533: DM_Complex *mesh = (DM_Complex *) dm->data;
4534: SieveLabel next = mesh->labels;
4535: PetscBool flg = PETSC_FALSE;
4536: PetscInt depth;
4541: if (stratumValue < 0) {
4542: DMComplexGetChart(dm, start, end);
4543: return(0);
4544: } else {
4545: PetscInt pStart, pEnd;
4547: if (start) {*start = 0;}
4548: if (end) {*end = 0;}
4549: DMComplexGetChart(dm, &pStart, &pEnd);
4550: if (pStart == pEnd) {return(0);}
4551: }
4552: DMComplexHasLabel(dm, "depth", &flg);
4553: if (!flg) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");
4554: /* We should have a generic GetLabel() and a Label class */
4555: while(next) {
4556: PetscStrcmp("depth", next->name, &flg);
4557: if (flg) break;
4558: next = next->next;
4559: }
4560: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
4561: depth = stratumValue;
4562: if ((depth < 0) || (depth >= next->numStrata)) {
4563: if (start) {*start = 0;}
4564: if (end) {*end = 0;}
4565: } else {
4566: if (start) {*start = next->points[next->stratumOffsets[depth]];}
4567: if (end) {*end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;}
4568: }
4569: return(0);
4570: }
4574: PetscErrorCode DMComplexGetHeightStratum(DM dm, PetscInt stratumValue, PetscInt *start, PetscInt *end) {
4575: DM_Complex *mesh = (DM_Complex *) dm->data;
4576: SieveLabel next = mesh->labels;
4577: PetscBool flg = PETSC_FALSE;
4578: PetscInt depth;
4583: if (stratumValue < 0) {
4584: DMComplexGetChart(dm, start, end);
4585: } else {
4586: PetscInt pStart, pEnd;
4588: if (start) {*start = 0;}
4589: if (end) {*end = 0;}
4590: DMComplexGetChart(dm, &pStart, &pEnd);
4591: if (pStart == pEnd) {return(0);}
4592: }
4593: DMComplexHasLabel(dm, "depth", &flg);
4594: if (!flg) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "No label named depth was found");
4595: /* We should have a generic GetLabel() and a Label class */
4596: while(next) {
4597: PetscStrcmp("depth", next->name, &flg);
4598: if (flg) break;
4599: next = next->next;
4600: }
4601: /* Strata are sorted and contiguous -- In addition, depth/height is either full or 1-level */
4602: depth = next->stratumValues[next->numStrata-1] - stratumValue;
4603: if ((depth < 0) || (depth >= next->numStrata)) {
4604: if (start) {*start = 0;}
4605: if (end) {*end = 0;}
4606: } else {
4607: if (start) {*start = next->points[next->stratumOffsets[depth]];}
4608: if (end) {*end = next->points[next->stratumOffsets[depth]+next->stratumSizes[depth]-1]+1;}
4609: }
4610: return(0);
4611: }
4615: PetscErrorCode DMComplexCreateSection(DM dm, PetscInt dim, PetscInt numFields, PetscInt numComp[], PetscInt numDof[], PetscInt numBC, PetscInt bcField[], IS bcPoints[], PetscSection *section) {
4616: PetscInt *numDofTot, *maxConstraints;
4617: PetscInt pStart = 0, pEnd = 0;
4618: PetscInt p, d, f, bc;
4622: PetscMalloc2(dim+1,PetscInt,&numDofTot,numFields+1,PetscInt,&maxConstraints);
4623: for(f = 0; f <= numFields; ++f) {maxConstraints[f] = 0;}
4624: for(d = 0; d <= dim; ++d) {
4625: numDofTot[d] = 0;
4626: for(f = 0; f < numFields; ++f) {
4627: numDofTot[d] += numDof[f*(dim+1)+d];
4628: }
4629: }
4630: PetscSectionCreate(((PetscObject) dm)->comm, section);
4631: if (numFields > 1) {
4632: PetscSectionSetNumFields(*section, numFields);
4633: if (numComp) {
4634: for(f = 0; f < numFields; ++f) {
4635: PetscSectionSetFieldComponents(*section, f, numComp[f]);
4636: }
4637: }
4638: } else {
4639: numFields = 0;
4640: }
4641: DMComplexGetChart(dm, &pStart, &pEnd);
4642: PetscSectionSetChart(*section, pStart, pEnd);
4643: for(d = 0; d <= dim; ++d) {
4644: DMComplexGetDepthStratum(dm, d, &pStart, &pEnd);
4645: for(p = pStart; p < pEnd; ++p) {
4646: for(f = 0; f < numFields; ++f) {
4647: PetscSectionSetFieldDof(*section, p, f, numDof[f*(dim+1)+d]);
4648: }
4649: PetscSectionSetDof(*section, p, numDofTot[d]);
4650: }
4651: }
4652: if (numBC) {
4653: for(bc = 0; bc < numBC; ++bc) {
4654: PetscInt field = 0;
4655: const PetscInt *idx;
4656: PetscInt n, i;
4658: if (numFields) {field = bcField[bc];}
4659: ISGetLocalSize(bcPoints[bc], &n);
4660: ISGetIndices(bcPoints[bc], &idx);
4661: for(i = 0; i < n; ++i) {
4662: const PetscInt p = idx[i];
4663: PetscInt depth, numConst;
4665: DMComplexGetLabelValue(dm, "depth", p, &depth);
4666: numConst = numDof[field*(dim+1)+depth];
4667: maxConstraints[field] = PetscMax(maxConstraints[field], numConst);
4668: if (numFields) {
4669: PetscSectionSetFieldConstraintDof(*section, p, field, numConst);
4670: }
4671: PetscSectionAddConstraintDof(*section, p, numConst);
4672: }
4673: ISRestoreIndices(bcPoints[bc], &idx);
4674: }
4675: for(f = 0; f < numFields; ++f) {
4676: maxConstraints[numFields] += maxConstraints[f];
4677: }
4678: }
4679: PetscSectionSetUp(*section);
4680: if (maxConstraints[numFields]) {
4681: PetscInt *indices;
4683: PetscMalloc(maxConstraints[numFields] * sizeof(PetscInt), &indices);
4684: PetscSectionGetChart(*section, &pStart, &pEnd);
4685: for(p = pStart; p < pEnd; ++p) {
4686: PetscInt cDof;
4688: PetscSectionGetConstraintDof(*section, p, &cDof);
4689: if (cDof) {
4690: if (cDof > maxConstraints[numFields]) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_LIB, "Likely memory corruption, point %D cDof %D > maxConstraints %D", p, cDof, maxConstraints);
4691: if (numFields) {
4692: PetscInt numConst = 0, fOff = 0;
4694: for(f = 0; f < numFields; ++f) {
4695: PetscInt cfDof, fDof;
4697: PetscSectionGetFieldDof(*section, p, f, &fDof);
4698: PetscSectionGetFieldConstraintDof(*section, p, f, &cfDof);
4699: for(d = 0; d < cfDof; ++d) {
4700: indices[numConst+d] = fOff+d;
4701: }
4702: PetscSectionSetFieldConstraintIndices(*section, p, f, &indices[numConst]);
4703: numConst += cfDof;
4704: fOff += fDof;
4705: }
4706: if (cDof != numConst) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_LIB, "Total number of field constraints %D should be %D", numConst, cDof);
4707: } else {
4708: for(d = 0; d < cDof; ++d) {
4709: indices[d] = d;
4710: }
4711: }
4712: PetscSectionSetConstraintIndices(*section, p, indices);
4713: }
4714: }
4715: PetscFree(indices);
4716: }
4717: PetscFree2(numDofTot,maxConstraints);
4718: {
4719: PetscBool view = PETSC_FALSE;
4721: PetscOptionsHasName(((PetscObject) dm)->prefix, "-section_view", &view);
4722: if (view) {
4723: PetscSectionView(*section, PETSC_VIEWER_STDOUT_WORLD);
4724: }
4725: }
4726: return(0);
4727: }
4731: PetscErrorCode DMComplexGetCoordinateSection(DM dm, PetscSection *section) {
4732: DM_Complex *mesh = (DM_Complex *) dm->data;
4737: *section = mesh->coordSection;
4738: return(0);
4739: }
4743: PetscErrorCode DMComplexSetCoordinateSection(DM dm, PetscSection section) {
4744: DM_Complex *mesh = (DM_Complex *) dm->data;
4749: PetscSectionDestroy(&mesh->coordSection);
4750: mesh->coordSection = section;
4751: return(0);
4752: }
4756: PetscErrorCode DMComplexGetConeSection(DM dm, PetscSection *section) {
4757: DM_Complex *mesh = (DM_Complex *) dm->data;
4761: if (section) *section = mesh->coneSection;
4762: return(0);
4763: }
4767: PetscErrorCode DMComplexGetCones(DM dm, PetscInt *cones[]) {
4768: DM_Complex *mesh = (DM_Complex *) dm->data;
4772: if (cones) *cones = mesh->cones;
4773: return(0);
4774: }
4778: PetscErrorCode DMComplexGetConeOrientations(DM dm, PetscInt *coneOrientations[]) {
4779: DM_Complex *mesh = (DM_Complex *) dm->data;
4783: if (coneOrientations) *coneOrientations = mesh->coneOrientations;
4784: return(0);
4785: }
4789: PetscErrorCode DMComplexGetCoordinateVec(DM dm, Vec *coordinates) {
4790: DM_Complex *mesh = (DM_Complex *) dm->data;
4795: *coordinates = mesh->coordinates;
4796: return(0);
4797: }
4799: /******************************** FEM Support **********************************/
4803: /*@C
4804: DMComplexVecGetClosure - Get an array of the values on the closure of 'point'
4806: Not collective
4808: Input Parameters:
4809: + dm - The DM
4810: . section - The section describing the layout in v, or PETSC_NULL to use the default section
4811: . v - The local vector
4812: - point - The sieve point in the DM
4814: Output Parameters:
4815: . values - The array of values, which is a borrowed array and should not be freed
4817: Level: intermediate
4819: .seealso DMComplexVecSetClosure(), DMComplexMatSetClosure()
4820: @*/
4821: PetscErrorCode DMComplexVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar *values[]) {
4822: PetscScalar *array, *vArray;
4823: PetscInt *points = PETSC_NULL;
4824: PetscInt offsets[32];
4825: PetscInt numFields, size, numPoints, pStart, pEnd, p, q, f;
4826: PetscErrorCode ierr;
4831: if (!section) {
4832: DMGetDefaultSection(dm, §ion);
4833: }
4834: PetscSectionGetNumFields(section, &numFields);
4835: if (numFields > 32) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 32", numFields);
4836: PetscMemzero(offsets, 32 * sizeof(PetscInt));
4837: DMComplexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
4838: /* Compress out points not in the section */
4839: PetscSectionGetChart(section, &pStart, &pEnd);
4840: for(p = 0, q = 0; p < numPoints*2; p += 2) {
4841: if ((points[p] >= pStart) && (points[p] < pEnd)) {
4842: points[q*2] = points[p];
4843: points[q*2+1] = points[p+1];
4844: ++q;
4845: }
4846: }
4847: numPoints = q;
4848: for(p = 0, size = 0; p < numPoints*2; p += 2) {
4849: PetscInt dof, fdof;
4851: PetscSectionGetDof(section, points[p], &dof);
4852: for(f = 0; f < numFields; ++f) {
4853: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4854: offsets[f+1] += fdof;
4855: }
4856: size += dof;
4857: }
4858: DMGetWorkArray(dm, 2*size, &array);
4859: VecGetArray(v, &vArray);
4860: for(p = 0; p < numPoints*2; p += 2) {
4861: PetscInt o = points[p+1];
4862: PetscInt dof, off, d;
4863: PetscScalar *varr;
4865: PetscSectionGetDof(section, points[p], &dof);
4866: PetscSectionGetOffset(section, points[p], &off);
4867: varr = &vArray[off];
4868: if (numFields) {
4869: PetscInt fdof, foff, fcomp, f, c;
4871: for(f = 0, foff = 0; f < numFields; ++f) {
4872: PetscSectionGetFieldDof(section, points[p], f, &fdof);
4873: if (o >= 0) {
4874: for(d = 0; d < fdof; ++d, ++offsets[f]) {
4875: array[offsets[f]] = varr[foff+d];
4876: }
4877: } else {
4878: PetscSectionGetFieldComponents(section, f, &fcomp);
4879: for(d = fdof/fcomp-1; d >= 0; --d) {
4880: for(c = 0; c < fcomp; ++c, ++offsets[f]) {
4881: array[offsets[f]] = varr[foff+d*fcomp+c];
4882: }
4883: }
4884: }
4885: foff += fdof;
4886: }
4887: } else {
4888: if (o >= 0) {
4889: for(d = 0; d < dof; ++d, ++offsets[0]) {
4890: array[offsets[0]] = varr[d];
4891: }
4892: } else {
4893: for(d = dof-1; d >= 0; --d, ++offsets[0]) {
4894: array[offsets[0]] = varr[d];
4895: }
4896: }
4897: }
4898: }
4899: VecRestoreArray(v, &vArray);
4900: *values = array;
4901: return(0);
4902: }
4904: PETSC_STATIC_INLINE void add (PetscScalar *x, PetscScalar y) {*x += y;}
4905: PETSC_STATIC_INLINE void insert(PetscScalar *x, PetscScalar y) {*x = y;}
4909: PetscErrorCode updatePoint_private(PetscSection section, PetscInt point, PetscInt dof, void (*fuse)(PetscScalar *, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[])
4910: {
4911: PetscInt cdof; /* The number of constraints on this point */
4912: PetscInt *cdofs; /* The indices of the constrained dofs on this point */
4913: PetscScalar *a;
4914: PetscInt off, cind = 0, k;
4918: PetscSectionGetConstraintDof(section, point, &cdof);
4919: PetscSectionGetOffset(section, point, &off);
4920: a = &array[off];
4921: if (!cdof || setBC) {
4922: if (orientation >= 0) {
4923: for(k = 0; k < dof; ++k) {
4924: fuse(&a[k], values[k]);
4925: }
4926: } else {
4927: for(k = 0; k < dof; ++k) {
4928: fuse(&a[k], values[dof-k-1]);
4929: }
4930: }
4931: } else {
4932: PetscSectionGetConstraintIndices(section, point, &cdofs);
4933: if (orientation >= 0) {
4934: for(k = 0; k < dof; ++k) {
4935: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4936: fuse(&a[k], values[k]);
4937: }
4938: } else {
4939: for(k = 0; k < dof; ++k) {
4940: if ((cind < cdof) && (k == cdofs[cind])) {++cind; continue;}
4941: fuse(&a[k], values[dof-k-1]);
4942: }
4943: }
4944: }
4945: return(0);
4946: }
4950: PetscErrorCode updatePointFields_private(PetscSection section, PetscInt point, PetscInt foffs[], void (*fuse)(PetscScalar *, PetscScalar), PetscBool setBC, PetscInt orientation, const PetscScalar values[], PetscScalar array[]) {
4951: PetscScalar *a;
4952: PetscInt numFields, off, foff, f;
4956: PetscSectionGetNumFields(section, &numFields);
4957: PetscSectionGetOffset(section, point, &off);
4958: a = &array[off];
4959: for(f = 0, foff = 0; f < numFields; ++f) {
4960: PetscInt fdof, fcomp, fcdof;
4961: PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
4962: PetscInt cind = 0, k, c;
4964: PetscSectionGetFieldComponents(section, f, &fcomp);
4965: PetscSectionGetFieldDof(section, point, f, &fdof);
4966: PetscSectionGetFieldConstraintDof(section, point, f, &fcdof);
4967: if (!fcdof || setBC) {
4968: if (orientation >= 0) {
4969: for(k = 0; k < fdof; ++k) {
4970: fuse(&a[foff+k], values[foffs[f]+k]);
4971: }
4972: } else {
4973: for(k = fdof/fcomp-1; k >= 0; --k) {
4974: for(c = 0; c < fcomp; ++c) {
4975: fuse(&a[foff+(fdof/fcomp-1-k)*fcomp+c], values[foffs[f]+k*fcomp+c]);
4976: }
4977: }
4978: }
4979: } else {
4980: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
4981: if (orientation >= 0) {
4982: for(k = 0; k < fdof; ++k) {
4983: if ((cind < fcdof) && (k == fcdofs[cind])) {++cind; continue;}
4984: fuse(&a[foff+k], values[foffs[f]+k]);
4985: }
4986: } else {
4987: for(k = fdof/fcomp-1; k >= 0; --k) {
4988: for(c = 0; c < fcomp; ++c) {
4989: if ((cind < fcdof) && (k*fcomp+c == fcdofs[cind])) {++cind; continue;}
4990: fuse(&a[foff+(fdof/fcomp-1-k)*fcomp+c], values[foffs[f]+k*fcomp+c]);
4991: }
4992: }
4993: }
4994: }
4995: foff += fdof;
4996: foffs[f] += fdof;
4997: }
4998: return(0);
4999: }
5003: /*@C
5004: DMComplexVecSetClosure - Set an array of the values on the closure of 'point'
5006: Not collective
5008: Input Parameters:
5009: + dm - The DM
5010: . section - The section describing the layout in v, or PETSC_NULL to use the default sectionw
5011: . v - The local vector
5012: . point - The sieve point in the DM
5013: . values - The array of values, which is a borrowed array and should not be freed
5014: - mode - The insert mode, where INSERT_ALL_VALUES and ADD_ALL_VALUES also overwrite boundary conditions
5016: Level: intermediate
5018: .seealso DMComplexVecGetClosure(), DMComplexMatSetClosure()
5019: @*/
5020: PetscErrorCode DMComplexVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt point, const PetscScalar values[], InsertMode mode) {
5021: PetscScalar *array;
5022: PetscInt *points = PETSC_NULL;
5023: PetscInt offsets[32];
5024: PetscInt numFields, numPoints, off, dof, pStart, pEnd, p, q, f;
5025: PetscErrorCode ierr;
5030: if (!section) {
5031: DMGetDefaultSection(dm, §ion);
5032: }
5033: PetscSectionGetNumFields(section, &numFields);
5034: if (numFields > 32) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 32", numFields);
5035: PetscMemzero(offsets, 32 * sizeof(PetscInt));
5036: DMComplexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5037: /* Compress out points not in the section */
5038: PetscSectionGetChart(section, &pStart, &pEnd);
5039: for(p = 0, q = 0; p < numPoints*2; p += 2) {
5040: if ((points[p] >= pStart) && (points[p] < pEnd)) {
5041: points[q*2] = points[p];
5042: points[q*2+1] = points[p+1];
5043: ++q;
5044: }
5045: }
5046: numPoints = q;
5047: for(p = 0; p < numPoints*2; p += 2) {
5048: PetscInt fdof;
5050: for(f = 0; f < numFields; ++f) {
5051: PetscSectionGetFieldDof(section, points[p], f, &fdof);
5052: offsets[f+1] += fdof;
5053: }
5054: }
5055: VecGetArray(v, &array);
5056: if (numFields) {
5057: switch(mode) {
5058: case INSERT_VALUES:
5059: for(p = 0; p < numPoints*2; p += 2) {
5060: PetscInt o = points[p+1];
5061: updatePointFields_private(section, points[p], offsets, insert, PETSC_FALSE, o, values, array);
5062: } break;
5063: case INSERT_ALL_VALUES:
5064: for(p = 0; p < numPoints*2; p += 2) {
5065: PetscInt o = points[p+1];
5066: updatePointFields_private(section, points[p], offsets, insert, PETSC_TRUE, o, values, array);
5067: } break;
5068: case ADD_VALUES:
5069: for(p = 0; p < numPoints*2; p += 2) {
5070: PetscInt o = points[p+1];
5071: updatePointFields_private(section, points[p], offsets, add, PETSC_FALSE, o, values, array);
5072: } break;
5073: case ADD_ALL_VALUES:
5074: for(p = 0; p < numPoints*2; p += 2) {
5075: PetscInt o = points[p+1];
5076: updatePointFields_private(section, points[p], offsets, add, PETSC_TRUE, o, values, array);
5077: } break;
5078: default:
5079: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5080: }
5081: } else {
5082: switch(mode) {
5083: case INSERT_VALUES:
5084: for(p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5085: PetscInt o = points[p+1];
5086: PetscSectionGetDof(section, points[p], &dof);
5087: updatePoint_private(section, points[p], dof, insert, PETSC_FALSE, o, &values[off], array);
5088: } break;
5089: case INSERT_ALL_VALUES:
5090: for(p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5091: PetscInt o = points[p+1];
5092: PetscSectionGetDof(section, points[p], &dof);
5093: updatePoint_private(section, points[p], dof, insert, PETSC_TRUE, o, &values[off], array);
5094: } break;
5095: case ADD_VALUES:
5096: for(p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5097: PetscInt o = points[p+1];
5098: PetscSectionGetDof(section, points[p], &dof);
5099: updatePoint_private(section, points[p], dof, add, PETSC_FALSE, o, &values[off], array);
5100: } break;
5101: case ADD_ALL_VALUES:
5102: for(p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5103: PetscInt o = points[p+1];
5104: PetscSectionGetDof(section, points[p], &dof);
5105: updatePoint_private(section, points[p], dof, add, PETSC_TRUE, o, &values[off], array);
5106: } break;
5107: default:
5108: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insert mode %D", mode);
5109: }
5110: }
5111: VecRestoreArray(v, &array);
5112: return(0);
5113: }
5117: PetscErrorCode DMComplexPrintMatSetValues(Mat A, PetscInt point, PetscInt numIndices, const PetscInt indices[], PetscScalar values[])
5118: {
5119: PetscMPIInt rank;
5120: PetscInt i, j;
5124: MPI_Comm_rank(((PetscObject) A)->comm, &rank);
5125: PetscPrintf(PETSC_COMM_SELF, "[%D]mat for sieve point %D\n", rank, point);
5126: for(i = 0; i < numIndices; i++) {
5127: PetscPrintf(PETSC_COMM_SELF, "[%D]mat indices[%D] = %D\n", rank, i, indices[i]);
5128: }
5129: for(i = 0; i < numIndices; i++) {
5130: PetscPrintf(PETSC_COMM_SELF, "[%D]", rank);
5131: for(j = 0; j < numIndices; j++) {
5132: #ifdef PETSC_USE_COMPLEX
5133: PetscPrintf(PETSC_COMM_SELF, " (%G,%G)", PetscRealPart(values[i*numIndices+j]), PetscImaginaryPart(values[i*numIndices+j]));
5134: #else
5135: PetscPrintf(PETSC_COMM_SELF, " %G", values[i*numIndices+j]);
5136: #endif
5137: }
5138: PetscPrintf(PETSC_COMM_SELF, "\n");
5139: }
5140: return(0);
5141: }
5145: /* . off - The global offset of this point */
5146: PetscErrorCode indicesPoint_private(PetscSection section, PetscInt point, PetscInt dof, PetscInt off, PetscBool setBC, PetscInt orientation, PetscInt indices[]) {
5147: PetscInt cdof; /* The number of constraints on this point */
5148: PetscInt *cdofs; /* The indices of the constrained dofs on this point */
5149: PetscInt cind = 0, k;
5153: PetscSectionGetDof(section, point, &dof);
5154: PetscSectionGetConstraintDof(section, point, &cdof);
5155: if (!cdof || setBC) {
5156: if (orientation >= 0) {
5157: for(k = 0; k < dof; ++k) {
5158: indices[k] = off+k;
5159: }
5160: } else {
5161: for(k = 0; k < dof; ++k) {
5162: indices[dof-k-1] = off+k;
5163: }
5164: }
5165: } else {
5166: PetscSectionGetConstraintIndices(section, point, &cdofs);
5167: if (orientation >= 0) {
5168: for(k = 0; k < dof; ++k) {
5169: if ((cind < cdof) && (k == cdofs[cind])) {
5170: /* Insert check for returning constrained indices */
5171: indices[k] = -(off+k+1);
5172: ++cind;
5173: } else {
5174: indices[k] = off+k;
5175: }
5176: }
5177: } else {
5178: for(k = 0; k < dof; ++k) {
5179: if ((cind < cdof) && (k == cdofs[cind])) {
5180: /* Insert check for returning constrained indices */
5181: indices[dof-k-1] = -(off+k+1);
5182: ++cind;
5183: } else {
5184: indices[dof-k-1] = off+k;
5185: }
5186: }
5187: }
5188: }
5189: return(0);
5190: }
5194: /* . off - The global offset of this point */
5195: PetscErrorCode indicesPointFields_private(PetscSection section, PetscInt point, PetscInt off, PetscInt foffs[], PetscBool setBC, PetscInt orientation, PetscInt indices[]) {
5196: PetscInt numFields, foff, f;
5200: PetscSectionGetNumFields(section, &numFields);
5201: for(f = 0, foff = 0; f < numFields; ++f) {
5202: PetscInt fdof, fcomp, cfdof;
5203: PetscInt *fcdofs; /* The indices of the constrained dofs for field f on this point */
5204: PetscInt cind = 0, k, c;
5206: PetscSectionGetFieldComponents(section, f, &fcomp);
5207: PetscSectionGetFieldDof(section, point, f, &fdof);
5208: PetscSectionGetFieldConstraintDof(section, point, f, &cfdof);
5209: if (!cfdof || setBC) {
5210: if (orientation >= 0) {
5211: for(k = 0; k < fdof; ++k) {
5212: indices[foffs[f]+k] = off+foff+k;
5213: }
5214: } else {
5215: for(k = fdof/fcomp-1; k >= 0; --k) {
5216: for(c = 0; c < fcomp; ++c) {
5217: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
5218: }
5219: }
5220: }
5221: } else {
5222: PetscSectionGetFieldConstraintIndices(section, point, f, &fcdofs);
5223: if (orientation >= 0) {
5224: for(k = 0; k < fdof; ++k) {
5225: if ((cind < cfdof) && (k == fcdofs[cind])) {
5226: indices[foffs[f]+k] = -(off+foff+k+1);
5227: ++cind;
5228: } else {
5229: indices[foffs[f]+k] = off+foff+k;
5230: }
5231: }
5232: } else {
5233: for(k = fdof/fcomp-1; k >= 0; --k) {
5234: for(c = 0; c < fcomp; ++c) {
5235: if ((cind < cfdof) && ((fdof/fcomp-1-k)*fcomp+c == fcdofs[cind])) {
5236: indices[foffs[f]+k*fcomp+c] = -(off+foff+(fdof/fcomp-1-k)*fcomp+c+1);
5237: ++cind;
5238: } else {
5239: indices[foffs[f]+k*fcomp+c] = off+foff+(fdof/fcomp-1-k)*fcomp+c;
5240: }
5241: }
5242: }
5243: }
5244: }
5245: foff += fdof - cfdof;
5246: foffs[f] += fdof;
5247: }
5248: return(0);
5249: }
5253: PetscErrorCode DMComplexMatSetClosure(DM dm, PetscSection section, PetscSection globalSection, Mat A, PetscInt point, PetscScalar values[], InsertMode mode)
5254: {
5255: DM_Complex *mesh = (DM_Complex *) dm->data;
5256: PetscInt *points = PETSC_NULL;
5257: PetscInt *indices;
5258: PetscInt offsets[32];
5259: PetscInt numFields, numPoints, numIndices, dof, off, globalOff, pStart, pEnd, p, q, f;
5260: PetscBool useDefault = !section ? PETSC_TRUE : PETSC_FALSE;
5261: PetscBool useGlobalDefault = !globalSection ? PETSC_TRUE : PETSC_FALSE;
5262: PetscErrorCode ierr;
5267: if (useDefault) {
5268: DMGetDefaultSection(dm, §ion);
5269: }
5270: if (useGlobalDefault) {
5271: if (useDefault) {
5272: DMGetDefaultGlobalSection(dm, &globalSection);
5273: } else {
5274: PetscSectionCreateGlobalSection(section, dm->sf, &globalSection);
5275: }
5276: }
5277: PetscSectionGetNumFields(section, &numFields);
5278: if (numFields > 32) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of fields %D limited to 32", numFields);
5279: PetscMemzero(offsets, 32 * sizeof(PetscInt));
5280: DMComplexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);
5281: /* Compress out points not in the section */
5282: PetscSectionGetChart(section, &pStart, &pEnd);
5283: for(p = 0, q = 0; p < numPoints*2; p += 2) {
5284: if ((points[p] >= pStart) && (points[p] < pEnd)) {
5285: points[q*2] = points[p];
5286: points[q*2+1] = points[p+1];
5287: ++q;
5288: }
5289: }
5290: numPoints = q;
5291: for(p = 0, numIndices = 0; p < numPoints*2; p += 2) {
5292: PetscInt fdof;
5294: PetscSectionGetDof(section, points[p], &dof);
5295: for(f = 0; f < numFields; ++f) {
5296: PetscSectionGetFieldDof(section, points[p], f, &fdof);
5297: offsets[f+1] += fdof;
5298: }
5299: numIndices += dof;
5300: }
5301: DMGetWorkArray(dm, numIndices, (PetscScalar **) &indices);
5302: if (numFields) {
5303: for(p = 0; p < numPoints*2; p += 2) {
5304: PetscInt o = points[p+1];
5305: PetscSectionGetOffset(globalSection, points[p], &globalOff);
5306: indicesPointFields_private(section, points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE, o, indices);
5307: }
5308: } else {
5309: for(p = 0, off = 0; p < numPoints*2; p += 2, off += dof) {
5310: PetscInt o = points[p+1];
5311: PetscSectionGetOffset(globalSection, points[p], &globalOff);
5312: indicesPoint_private(section, points[p], dof, globalOff < 0 ? -(globalOff+1) : globalOff, PETSC_FALSE, o, &indices[off]);
5313: }
5314: }
5315: if (useGlobalDefault && !useDefault) {
5316: PetscSectionDestroy(&globalSection);
5317: }
5318: if (mesh->printSetValues) {DMComplexPrintMatSetValues(A, point, numIndices, indices, values);}
5319: MatSetValues(A, numIndices, indices, numIndices, indices, values, mode);
5320: if (ierr) {
5321: PetscMPIInt rank;
5322: PetscErrorCode ierr2;
5324: ierr2 = MPI_Comm_rank(((PetscObject) A)->comm, &rank);CHKERRQ(ierr2);
5325: ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%D]ERROR in DMComplexMatSetClosure\n", rank);CHKERRQ(ierr2);
5326: ierr2 = DMComplexPrintMatSetValues(A, point, numIndices, indices, values);CHKERRQ(ierr2);
5327:
5328: }
5329: return(0);
5330: }
5334: PetscErrorCode DMComplexComputeTriangleGeometry_private(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
5335: {
5336: DM_Complex *mesh = (DM_Complex *) dm->data;
5337: const PetscScalar *coords;
5338: const PetscInt dim = 2;
5339: PetscInt d, f;
5340: PetscErrorCode ierr;
5343: DMComplexVecGetClosure(dm, mesh->coordSection, mesh->coordinates, e, &coords);
5344: if (v0) {
5345: for(d = 0; d < dim; d++) {
5346: v0[d] = PetscRealPart(coords[d]);
5347: }
5348: }
5349: if (J) {
5350: for(d = 0; d < dim; d++) {
5351: for(f = 0; f < dim; f++) {
5352: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5353: }
5354: }
5355: *detJ = J[0]*J[3] - J[1]*J[2];
5356: #if 0
5357: if (detJ < 0.0) {
5358: const PetscReal xLength = mesh->periodicity[0];
5360: if (xLength != 0.0) {
5361: PetscReal v0x = coords[0*dim+0];
5363: if (v0x == 0.0) {
5364: v0x = v0[0] = xLength;
5365: }
5366: for(f = 0; f < dim; f++) {
5367: const PetscReal px = coords[(f+1)*dim+0] == 0.0 ? xLength : coords[(f+1)*dim+0];
5369: J[0*dim+f] = 0.5*(px - v0x);
5370: }
5371: }
5372: detJ = J[0]*J[3] - J[1]*J[2];
5373: }
5374: #endif
5375: PetscLogFlops(8.0 + 3.0);
5376: }
5377: if (invJ) {
5378: const PetscReal invDet = 1.0/(*detJ);
5380: invJ[0] = invDet*J[3];
5381: invJ[1] = -invDet*J[1];
5382: invJ[2] = -invDet*J[2];
5383: invJ[3] = invDet*J[0];
5384: PetscLogFlops(5.0);
5385: }
5386: return(0);
5387: }
5391: PetscErrorCode DMComplexComputeRectangleGeometry_private(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
5392: {
5393: DM_Complex *mesh = (DM_Complex *) dm->data;
5394: const PetscScalar *coords;
5395: const PetscInt dim = 2;
5396: PetscInt d, f;
5397: PetscErrorCode ierr;
5400: DMComplexVecGetClosure(dm, mesh->coordSection, mesh->coordinates, e, &coords);
5401: if (v0) {
5402: for(d = 0; d < dim; d++) {
5403: v0[d] = PetscRealPart(coords[d]);
5404: }
5405: }
5406: if (J) {
5407: for(d = 0; d < dim; d++) {
5408: for(f = 0; f < dim; f++) {
5409: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5410: }
5411: }
5412: *detJ = J[0]*J[3] - J[1]*J[2];
5413: PetscLogFlops(8.0 + 3.0);
5414: }
5415: if (invJ) {
5416: const PetscReal invDet = 1.0/(*detJ);
5418: invJ[0] = invDet*J[3];
5419: invJ[1] = -invDet*J[1];
5420: invJ[2] = -invDet*J[2];
5421: invJ[3] = invDet*J[0];
5422: PetscLogFlopsNoError(5.0);
5423: }
5424: *detJ *= 2.0;
5425: return(0);
5426: }
5430: PetscErrorCode DMComplexComputeTetrahedronGeometry_private(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
5431: {
5432: DM_Complex *mesh = (DM_Complex *) dm->data;
5433: const PetscScalar *coords;
5434: const PetscInt dim = 3;
5435: PetscInt d, f;
5436: PetscErrorCode ierr;
5439: DMComplexVecGetClosure(dm, mesh->coordSection, mesh->coordinates, e, &coords);
5440: if (v0) {
5441: for(d = 0; d < dim; d++) {
5442: v0[d] = PetscRealPart(coords[d]);
5443: }
5444: }
5445: if (J) {
5446: for(d = 0; d < dim; d++) {
5447: for(f = 0; f < dim; f++) {
5448: J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5449: }
5450: }
5451: /* ??? This does not work with CTetGen: The minus sign is here since I orient the first face to get the outward normal */
5452: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
5453: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
5454: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
5455: PetscLogFlops(18.0 + 12.0);
5456: }
5457: if (invJ) {
5458: const PetscReal invDet = -1.0/(*detJ);
5460: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
5461: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
5462: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
5463: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
5464: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
5465: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
5466: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
5467: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
5468: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
5469: PetscLogFlops(37.0);
5470: }
5471: return(0);
5472: }
5476: PetscErrorCode DMComplexComputeHexahedronGeometry_private(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
5477: {
5478: DM_Complex *mesh = (DM_Complex *) dm->data;
5479: const PetscScalar *coords;
5480: const PetscInt dim = 3;
5481: PetscInt d;
5482: PetscErrorCode ierr;
5485: DMComplexVecGetClosure(dm, mesh->coordSection, mesh->coordinates, e, &coords);
5486: if (v0) {
5487: for(d = 0; d < dim; d++) {
5488: v0[d] = PetscRealPart(coords[d]);
5489: }
5490: }
5491: if (J) {
5492: for(d = 0; d < dim; d++) {
5493: J[d*dim+0] = 0.5*(PetscRealPart(coords[(0+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5494: J[d*dim+1] = 0.5*(PetscRealPart(coords[(1+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5495: J[d*dim+2] = 0.5*(PetscRealPart(coords[(3+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
5496: }
5497: *detJ = (J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
5498: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
5499: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]));
5500: PetscLogFlops(18.0 + 12.0);
5501: }
5502: if (invJ) {
5503: const PetscReal invDet = -1.0/(*detJ);
5505: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
5506: invJ[0*3+1] = invDet*(J[0*3+2]*J[2*3+1] - J[0*3+1]*J[2*3+2]);
5507: invJ[0*3+2] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
5508: invJ[1*3+0] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
5509: invJ[1*3+1] = invDet*(J[0*3+0]*J[2*3+2] - J[0*3+2]*J[2*3+0]);
5510: invJ[1*3+2] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
5511: invJ[2*3+0] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
5512: invJ[2*3+1] = invDet*(J[0*3+1]*J[2*3+0] - J[0*3+0]*J[2*3+1]);
5513: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
5514: PetscLogFlops(37.0);
5515: }
5516: *detJ *= 8.0;
5517: return(0);
5518: }
5522: PetscErrorCode DMComplexComputeCellGeometry(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ) {
5523: PetscInt dim, maxConeSize;
5527: DMComplexGetDimension(dm, &dim);
5528: DMComplexGetMaxSizes(dm, &maxConeSize, PETSC_NULL);
5529: switch(dim) {
5530: case 2:
5531: switch(maxConeSize) {
5532: case 3:
5533: DMComplexComputeTriangleGeometry_private(dm, cell, v0, J, invJ, detJ);
5534: break;
5535: case 4:
5536: DMComplexComputeRectangleGeometry_private(dm, cell, v0, J, invJ, detJ);
5537: break;
5538: default:
5539: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Unsupported number of cell vertices %D for element geometry computation", maxConeSize);
5540: }
5541: break;
5542: case 3:
5543: switch(maxConeSize) {
5544: case 4:
5545: DMComplexComputeTetrahedronGeometry_private(dm, cell, v0, J, invJ, detJ);
5546: break;
5547: case 8:
5548: DMComplexComputeHexahedronGeometry_private(dm, cell, v0, J, invJ, detJ);
5549: break;
5550: default:
5551: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Unsupported number of cell vertices %D for element geometry computation", maxConeSize);
5552: }
5553: break;
5554: default:
5555: SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
5556: }
5557: return(0);
5558: }
5562: PetscErrorCode DMComplexGetFaceOrientation(DM dm, PetscInt cell, PetscInt numCorners, PetscInt indices[], PetscInt oppositeVertex, PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented) {
5563: MPI_Comm comm = ((PetscObject) dm)->comm;
5564: PetscBool posOrient = PETSC_FALSE;
5565: const PetscInt debug = 0;
5566: PetscInt cellDim, faceSize, f;
5569: DMComplexGetDimension(dm, &cellDim);
5570: if (debug) {PetscPrintf(comm, "cellDim: %d numCorners: %d\n", cellDim, numCorners);}
5572: if (cellDim == numCorners-1) {
5573: /* Simplices */
5574: faceSize = numCorners-1;
5575: posOrient = !(oppositeVertex%2) ? PETSC_TRUE : PETSC_FALSE;
5576: } else if (cellDim == 1 && numCorners == 3) {
5577: /* Quadratic line */
5578: faceSize = 1;
5579: posOrient = PETSC_TRUE;
5580: } else if (cellDim == 2 && numCorners == 4) {
5581: /* Quads */
5582: faceSize = 2;
5583: if ((indices[1] > indices[0]) && (indices[1] - indices[0] == 1)) {
5584: posOrient = PETSC_TRUE;
5585: } else if ((indices[0] == 3) && (indices[1] == 0)) {
5586: posOrient = PETSC_TRUE;
5587: } else {
5588: if (((indices[0] > indices[1]) && (indices[0] - indices[1] == 1)) || ((indices[0] == 0) && (indices[1] == 3))) {
5589: posOrient = PETSC_FALSE;
5590: } else {
5591: SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossedge");
5592: }
5593: }
5594: } else if (cellDim == 2 && numCorners == 6) {
5595: /* Quadratic triangle (I hate this) */
5596: /* Edges are determined by the first 2 vertices (corners of edges) */
5597: const PetscInt faceSizeTri = 3;
5598: PetscInt sortedIndices[3], i, iFace;
5599: PetscBool found = PETSC_FALSE;
5600: PetscInt faceVerticesTriSorted[9] = {
5601: 0, 3, 4, /* bottom */
5602: 1, 4, 5, /* right */
5603: 2, 3, 5, /* left */
5604: };
5605: PetscInt faceVerticesTri[9] = {
5606: 0, 3, 4, /* bottom */
5607: 1, 4, 5, /* right */
5608: 2, 5, 3, /* left */
5609: };
5611: faceSize = faceSizeTri;
5612: for(i = 0; i < faceSizeTri; ++i) sortedIndices[i] = indices[i];
5613: PetscSortInt(faceSizeTri, sortedIndices);
5614: for(iFace = 0; iFace < 4; ++iFace) {
5615: const PetscInt ii = iFace*faceSizeTri;
5616: PetscInt fVertex, cVertex;
5618: if ((sortedIndices[0] == faceVerticesTriSorted[ii+0]) &&
5619: (sortedIndices[1] == faceVerticesTriSorted[ii+1])) {
5620: for(fVertex = 0; fVertex < faceSizeTri; ++fVertex) {
5621: for(cVertex = 0; cVertex < faceSizeTri; ++cVertex) {
5622: if (indices[cVertex] == faceVerticesTri[ii+fVertex]) {
5623: faceVertices[fVertex] = origVertices[cVertex];
5624: break;
5625: }
5626: }
5627: }
5628: found = PETSC_TRUE;
5629: break;
5630: }
5631: }
5632: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tri crossface");
5633: if (posOriented) {*posOriented = PETSC_TRUE;}
5634: return(0);
5635: } else if (cellDim == 2 && numCorners == 9) {
5636: /* Quadratic quad (I hate this) */
5637: /* Edges are determined by the first 2 vertices (corners of edges) */
5638: const PetscInt faceSizeQuad = 3;
5639: PetscInt sortedIndices[3], i, iFace;
5640: PetscBool found = PETSC_FALSE;
5641: PetscInt faceVerticesQuadSorted[12] = {
5642: 0, 1, 4, /* bottom */
5643: 1, 2, 5, /* right */
5644: 2, 3, 6, /* top */
5645: 0, 3, 7, /* left */
5646: };
5647: PetscInt faceVerticesQuad[12] = {
5648: 0, 1, 4, /* bottom */
5649: 1, 2, 5, /* right */
5650: 2, 3, 6, /* top */
5651: 3, 0, 7, /* left */
5652: };
5654: faceSize = faceSizeQuad;
5655: for(i = 0; i < faceSizeQuad; ++i) sortedIndices[i] = indices[i];
5656: PetscSortInt(faceSizeQuad, sortedIndices);
5657: for(iFace = 0; iFace < 4; ++iFace) {
5658: const PetscInt ii = iFace*faceSizeQuad;
5659: PetscInt fVertex, cVertex;
5661: if ((sortedIndices[0] == faceVerticesQuadSorted[ii+0]) &&
5662: (sortedIndices[1] == faceVerticesQuadSorted[ii+1])) {
5663: for(fVertex = 0; fVertex < faceSizeQuad; ++fVertex) {
5664: for(cVertex = 0; cVertex < faceSizeQuad; ++cVertex) {
5665: if (indices[cVertex] == faceVerticesQuad[ii+fVertex]) {
5666: faceVertices[fVertex] = origVertices[cVertex];
5667: break;
5668: }
5669: }
5670: }
5671: found = PETSC_TRUE;
5672: break;
5673: }
5674: }
5675: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid quad crossface");
5676: if (posOriented) {*posOriented = PETSC_TRUE;}
5677: return(0);
5678: } else if (cellDim == 3 && numCorners == 8) {
5679: /* Hexes
5680: A hex is two oriented quads with the normal of the first
5681: pointing up at the second.
5683: 7---6
5684: /| /|
5685: 4---5 |
5686: | 3-|-2
5687: |/ |/
5688: 0---1
5690: Faces are determined by the first 4 vertices (corners of faces) */
5691: const PetscInt faceSizeHex = 4;
5692: PetscInt sortedIndices[4], i, iFace;
5693: PetscBool found = PETSC_FALSE;
5694: PetscInt faceVerticesHexSorted[24] = {
5695: 0, 1, 2, 3, /* bottom */
5696: 4, 5, 6, 7, /* top */
5697: 0, 1, 4, 5, /* front */
5698: 1, 2, 5, 6, /* right */
5699: 2, 3, 6, 7, /* back */
5700: 0, 3, 4, 7, /* left */
5701: };
5702: PetscInt faceVerticesHex[24] = {
5703: 3, 2, 1, 0, /* bottom */
5704: 4, 5, 6, 7, /* top */
5705: 0, 1, 5, 4, /* front */
5706: 1, 2, 6, 5, /* right */
5707: 2, 3, 7, 6, /* back */
5708: 3, 0, 4, 7, /* left */
5709: };
5711: faceSize = faceSizeHex;
5712: for(i = 0; i < faceSizeHex; ++i) sortedIndices[i] = indices[i];
5713: PetscSortInt(faceSizeHex, sortedIndices);
5714: for(iFace = 0; iFace < 6; ++iFace) {
5715: const PetscInt ii = iFace*faceSizeHex;
5716: PetscInt fVertex, cVertex;
5718: if ((sortedIndices[0] == faceVerticesHexSorted[ii+0]) &&
5719: (sortedIndices[1] == faceVerticesHexSorted[ii+1]) &&
5720: (sortedIndices[2] == faceVerticesHexSorted[ii+2]) &&
5721: (sortedIndices[3] == faceVerticesHexSorted[ii+3])) {
5722: for(fVertex = 0; fVertex < faceSizeHex; ++fVertex) {
5723: for(cVertex = 0; cVertex < faceSizeHex; ++cVertex) {
5724: if (indices[cVertex] == faceVerticesHex[ii+fVertex]) {
5725: faceVertices[fVertex] = origVertices[cVertex];
5726: break;
5727: }
5728: }
5729: }
5730: found = PETSC_TRUE;
5731: break;
5732: }
5733: }
5734: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");
5735: if (posOriented) {*posOriented = PETSC_TRUE;}
5736: return(0);
5737: } else if (cellDim == 3 && numCorners == 10) {
5738: /* Quadratic tet */
5739: /* Faces are determined by the first 3 vertices (corners of faces) */
5740: const PetscInt faceSizeTet = 6;
5741: PetscInt sortedIndices[6], i, iFace;
5742: PetscBool found = PETSC_FALSE;
5743: PetscInt faceVerticesTetSorted[24] = {
5744: 0, 1, 2, 6, 7, 8, /* bottom */
5745: 0, 3, 4, 6, 7, 9, /* front */
5746: 1, 4, 5, 7, 8, 9, /* right */
5747: 2, 3, 5, 6, 8, 9, /* left */
5748: };
5749: PetscInt faceVerticesTet[24] = {
5750: 0, 1, 2, 6, 7, 8, /* bottom */
5751: 0, 4, 3, 6, 7, 9, /* front */
5752: 1, 5, 4, 7, 8, 9, /* right */
5753: 2, 3, 5, 8, 6, 9, /* left */
5754: };
5756: faceSize = faceSizeTet;
5757: for(i = 0; i < faceSizeTet; ++i) sortedIndices[i] = indices[i];
5758: PetscSortInt(faceSizeTet, sortedIndices);
5759: for(iFace=0; iFace < 6; ++iFace) {
5760: const PetscInt ii = iFace*faceSizeTet;
5761: PetscInt fVertex, cVertex;
5763: if ((sortedIndices[0] == faceVerticesTetSorted[ii+0]) &&
5764: (sortedIndices[1] == faceVerticesTetSorted[ii+1]) &&
5765: (sortedIndices[2] == faceVerticesTetSorted[ii+2]) &&
5766: (sortedIndices[3] == faceVerticesTetSorted[ii+3])) {
5767: for(fVertex = 0; fVertex < faceSizeTet; ++fVertex) {
5768: for(cVertex = 0; cVertex < faceSizeTet; ++cVertex) {
5769: if (indices[cVertex] == faceVerticesTet[ii+fVertex]) {
5770: faceVertices[fVertex] = origVertices[cVertex];
5771: break;
5772: }
5773: }
5774: }
5775: found = PETSC_TRUE;
5776: break;
5777: }
5778: }
5779: if (!found) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid tet crossface");
5780: if (posOriented) {*posOriented = PETSC_TRUE;}
5781: return(0);
5782: } else if (cellDim == 3 && numCorners == 27) {
5783: /* Quadratic hexes (I hate this)
5784: A hex is two oriented quads with the normal of the first
5785: pointing up at the second.
5787: 7---6
5788: /| /|
5789: 4---5 |
5790: | 3-|-2
5791: |/ |/
5792: 0---1
5794: Faces are determined by the first 4 vertices (corners of faces) */
5795: const PetscInt faceSizeQuadHex = 9;
5796: PetscInt sortedIndices[9], i, iFace;
5797: PetscBool found = PETSC_FALSE;
5798: PetscInt faceVerticesQuadHexSorted[54] = {
5799: 0, 1, 2, 3, 8, 9, 10, 11, 24, /* bottom */
5800: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
5801: 0, 1, 4, 5, 8, 12, 16, 17, 22, /* front */
5802: 1, 2, 5, 6, 9, 13, 17, 18, 21, /* right */
5803: 2, 3, 6, 7, 10, 14, 18, 19, 23, /* back */
5804: 0, 3, 4, 7, 11, 15, 16, 19, 20, /* left */
5805: };
5806: PetscInt faceVerticesQuadHex[54] = {
5807: 3, 2, 1, 0, 10, 9, 8, 11, 24, /* bottom */
5808: 4, 5, 6, 7, 12, 13, 14, 15, 25, /* top */
5809: 0, 1, 5, 4, 8, 17, 12, 16, 22, /* front */
5810: 1, 2, 6, 5, 9, 18, 13, 17, 21, /* right */
5811: 2, 3, 7, 6, 10, 19, 14, 18, 23, /* back */
5812: 3, 0, 4, 7, 11, 16, 15, 19, 20 /* left */
5813: };
5815: faceSize = faceSizeQuadHex;
5816: for(i = 0; i < faceSizeQuadHex; ++i) sortedIndices[i] = indices[i];
5817: PetscSortInt(faceSizeQuadHex, sortedIndices);
5818: for(iFace = 0; iFace < 6; ++iFace) {
5819: const PetscInt ii = iFace*faceSizeQuadHex;
5820: PetscInt fVertex, cVertex;
5822: if ((sortedIndices[0] == faceVerticesQuadHexSorted[ii+0]) &&
5823: (sortedIndices[1] == faceVerticesQuadHexSorted[ii+1]) &&
5824: (sortedIndices[2] == faceVerticesQuadHexSorted[ii+2]) &&
5825: (sortedIndices[3] == faceVerticesQuadHexSorted[ii+3])) {
5826: for(fVertex = 0; fVertex < faceSizeQuadHex; ++fVertex) {
5827: for(cVertex = 0; cVertex < faceSizeQuadHex; ++cVertex) {
5828: if (indices[cVertex] == faceVerticesQuadHex[ii+fVertex]) {
5829: faceVertices[fVertex] = origVertices[cVertex];
5830: break;
5831: }
5832: }
5833: }
5834: found = PETSC_TRUE;
5835: break;
5836: }
5837: }
5838: if (!found) {SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid hex crossface");}
5839: if (posOriented) {*posOriented = PETSC_TRUE;}
5840: return(0);
5841: } else SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Unknown cell type for faceOrientation().");
5842: if (!posOrient) {
5843: if (debug) {PetscPrintf(comm, " Reversing initial face orientation\n");}
5844: for(f = 0; f < faceSize; ++f) {
5845: faceVertices[f] = origVertices[faceSize-1 - f];
5846: }
5847: } else {
5848: if (debug) {PetscPrintf(comm, " Keeping initial face orientation\n");}
5849: for(f = 0; f < faceSize; ++f) {
5850: faceVertices[f] = origVertices[f];
5851: }
5852: }
5853: if (posOriented) {*posOriented = posOrient;}
5854: return(0);
5855: }
5859: PetscErrorCode DMComplexGetOrientedFace(DM dm, PetscInt cell, PetscInt faceSize, const PetscInt face[], PetscInt numCorners, PetscInt indices[], PetscInt origVertices[], PetscInt faceVertices[], PetscBool *posOriented)
5860: {
5861: const PetscInt *cone;
5862: PetscInt coneSize, v, f, v2;
5863: PetscInt oppositeVertex = -1;
5864: PetscErrorCode ierr;
5867: DMComplexGetConeSize(dm, cell, &coneSize);
5868: DMComplexGetCone(dm, cell, &cone);
5869: for(v = 0, v2 = 0; v < coneSize; ++v) {
5870: PetscBool found = PETSC_FALSE;
5872: for(f = 0; f < faceSize; ++f) {
5873: if (face[f] == cone[v]) {found = PETSC_TRUE; break;}
5874: }
5875: if (found) {
5876: indices[v2] = v;
5877: origVertices[v2] = cone[v];
5878: ++v2;
5879: } else {
5880: oppositeVertex = v;
5881: }
5882: }
5883: DMComplexGetFaceOrientation(dm, cell, numCorners, indices, oppositeVertex, origVertices, faceVertices, posOriented);
5884: return(0);
5885: }
5887: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
5888: {
5889: switch(i) {
5890: case 0:
5891: switch(j) {
5892: case 0: return 0;
5893: case 1:
5894: switch(k) {
5895: case 0: return 0;
5896: case 1: return 0;
5897: case 2: return 1;
5898: }
5899: case 2:
5900: switch(k) {
5901: case 0: return 0;
5902: case 1: return -1;
5903: case 2: return 0;
5904: }
5905: }
5906: case 1:
5907: switch(j) {
5908: case 0:
5909: switch(k) {
5910: case 0: return 0;
5911: case 1: return 0;
5912: case 2: return -1;
5913: }
5914: case 1: return 0;
5915: case 2:
5916: switch(k) {
5917: case 0: return 1;
5918: case 1: return 0;
5919: case 2: return 0;
5920: }
5921: }
5922: case 2:
5923: switch(j) {
5924: case 0:
5925: switch(k) {
5926: case 0: return 0;
5927: case 1: return 1;
5928: case 2: return 0;
5929: }
5930: case 1:
5931: switch(k) {
5932: case 0: return -1;
5933: case 1: return 0;
5934: case 2: return 0;
5935: }
5936: case 2: return 0;
5937: }
5938: }
5939: return 0;
5940: }
5944: /*@C
5945: DMComplexCreateRigidBody - create rigid body modes from coordinates
5947: Collective on DM
5949: Input Arguments:
5950: + dm - the DM
5951: - section - the local section associated with the rigid field, or PETSC_NULL for the default section
5952: - globalSection - the global section associated with the rigid field, or PETSC_NULL for the default section
5954: Output Argument:
5955: . sp - the null space
5957: Note: This is necessary to take account of Dirichlet conditions on the displacements
5959: Level: advanced
5961: .seealso: MatNullSpaceCreate()
5962: @*/
5963: PetscErrorCode DMComplexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
5964: {
5965: MPI_Comm comm = ((PetscObject) dm)->comm;
5966: Vec coordinates, localMode, mode[6];
5967: PetscSection coordSection;
5968: PetscScalar *coords;
5969: PetscInt dim, vStart, vEnd, v, n, m, d, i, j;
5973: DMComplexGetDimension(dm, &dim);
5974: if (dim == 1) {
5975: MatNullSpaceCreate(comm, PETSC_TRUE, 0, PETSC_NULL, sp);
5976: return(0);
5977: }
5978: if (!section) {DMGetDefaultSection(dm, §ion);}
5979: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
5980: PetscSectionGetConstrainedStorageSize(globalSection, &n);
5981: DMComplexGetDepthStratum(dm, 0, &vStart, &vEnd);
5982: DMComplexGetCoordinateSection(dm, &coordSection);
5983: DMComplexGetCoordinateVec(dm, &coordinates);
5984: m = (dim*(dim+1))/2;
5985: VecCreate(comm, &mode[0]);
5986: VecSetSizes(mode[0], n, PETSC_DETERMINE);
5987: VecSetUp(mode[0]);
5988: for(i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
5989: /* Assume P1 */
5990: DMGetLocalVector(dm, &localMode);
5991: for(d = 0; d < dim; ++d) {
5992: PetscScalar values[3] = {0.0, 0.0, 0.0};
5994: values[d] = 1.0;
5995: VecSet(localMode, 0.0);
5996: for(v = vStart; v < vEnd; ++v) {
5997: DMComplexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
5998: }
5999: DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
6000: DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
6001: }
6002: VecGetArray(coordinates, &coords);
6003: for(d = dim; d < dim*(dim+1)/2; ++d) {
6004: PetscInt i, j, k = dim > 2 ? d - dim : d;
6006: VecSet(localMode, 0.0);
6007: for(v = vStart; v < vEnd; ++v) {
6008: PetscScalar values[3] = {0.0, 0.0, 0.0};
6009: PetscInt off;
6011: PetscSectionGetOffset(coordSection, v, &off);
6012: for(i = 0; i < dim; ++i) {
6013: for(j = 0; j < dim; ++j) {
6014: values[j] += (PetscReal)epsilon(i, j, k)*coords[off+i];
6015: }
6016: }
6017: DMComplexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
6018: }
6019: DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
6020: DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
6021: }
6022: VecRestoreArray(coordinates, &coords);
6023: DMRestoreLocalVector(dm, &localMode);
6024: for(i = 0; i < dim; ++i) {VecNormalize(mode[i], PETSC_NULL);}
6025: /* Orthonormalize system */
6026: for(i = dim; i < m; ++i) {
6027: PetscScalar dots[6];
6029: VecMDot(mode[i], i, mode, dots);
6030: for(j = 0; j < i; ++j) dots[j] *= -1.0;
6031: VecMAXPY(mode[i], i, dots, mode);
6032: VecNormalize(mode[i], PETSC_NULL);
6033: }
6034: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
6035: for(i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
6036: return(0);
6037: }