Actual source code: complexcreate.c
petsc-3.3-p7 2013-05-11
1: #define PETSCDM_DLL
2: #include <petsc-private/compleximpl.h> /*I "petscdmcomplex.h" I*/
3: #include <petscdmda.h>
7: PetscErrorCode DMSetFromOptions_Complex(DM dm)
8: {
9: DM_Complex *mesh = (DM_Complex *) dm->data;
14: PetscOptionsHead("DMComplex Options");
15: /* Handle DMComplex refinement */
16: /* Handle associated vectors */
17: /* Handle viewing */
18: PetscOptionsBool("-dm_complex_print_set_values", "Output all set values info", "DMView", PETSC_FALSE, &mesh->printSetValues, PETSC_NULL);
19: PetscOptionsTail();
20: return(0);
21: }
25: /*
26: Simple square boundary:
28: 18--5-17--4--16
29: | | |
30: 6 10 3
31: | | |
32: 19-11-20--9--15
33: | | |
34: 7 8 2
35: | | |
36: 12--0-13--1--14
37: */
38: PetscErrorCode DMComplexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
39: {
40: DM_Complex *mesh = (DM_Complex *) dm->data;
41: PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
42: PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
43: PetscScalar *coords;
44: PetscInt coordSize;
45: PetscMPIInt rank;
46: PetscInt v, vx, vy;
50: MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
51: if (!rank) {
52: PetscInt e, ex, ey;
54: DMComplexSetChart(dm, 0, numEdges+numVertices);
55: for(e = 0; e < numEdges; ++e) {
56: DMComplexSetConeSize(dm, e, 2);
57: }
58: DMSetUp(dm); /* Allocate space for cones */
59: for(vy = 0; vy <= edges[1]; vy++) {
60: for(ex = 0; ex < edges[0]; ex++) {
61: PetscInt edge = vy*edges[0] + ex;
62: PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
63: PetscInt cone[2] = {vertex, vertex+1};
65: DMComplexSetCone(dm, edge, cone);
66: if ((vy == 0) || (vy == edges[1])) {
67: DMComplexSetLabelValue(dm, "marker", edge, 1);
68: DMComplexSetLabelValue(dm, "marker", cone[0], 1);
69: if (ex == edges[0]-1) {
70: DMComplexSetLabelValue(dm, "marker", cone[1], 1);
71: }
72: }
73: }
74: }
75: for(vx = 0; vx <= edges[0]; vx++) {
76: for(ey = 0; ey < edges[1]; ey++) {
77: PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1);
78: PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
79: PetscInt cone[2] = {vertex, vertex+edges[0]+1};
81: DMComplexSetCone(dm, edge, cone);
82: if ((vx == 0) || (vx == edges[0])) {
83: DMComplexSetLabelValue(dm, "marker", edge, 1);
84: DMComplexSetLabelValue(dm, "marker", cone[0], 1);
85: if (ey == edges[1]-1) {
86: DMComplexSetLabelValue(dm, "marker", cone[1], 1);
87: }
88: }
89: }
90: }
91: }
92: DMComplexSymmetrize(dm);
93: DMComplexStratify(dm);
94: /* Build coordinates */
95: PetscSectionSetChart(mesh->coordSection, numEdges, numEdges + numVertices);
96: for(v = numEdges; v < numEdges+numVertices; ++v) {
97: PetscSectionSetDof(mesh->coordSection, v, 2);
98: }
99: PetscSectionSetUp(mesh->coordSection);
100: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
101: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
102: VecSetFromOptions(mesh->coordinates);
103: VecGetArray(mesh->coordinates, &coords);
104: for(vy = 0; vy <= edges[1]; ++vy) {
105: for(vx = 0; vx <= edges[0]; ++vx) {
106: coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
107: coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
108: }
109: }
110: VecRestoreArray(mesh->coordinates, &coords);
111: return(0);
112: }
116: /*
117: Simple cubic boundary:
119: 2-------3
120: /| /|
121: 6-------7 |
122: | | | |
123: | 0-----|-1
124: |/ |/
125: 4-------5
126: */
127: PetscErrorCode DMComplexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
128: {
129: DM_Complex *mesh = (DM_Complex *) dm->data;
130: PetscInt numVertices = (faces[0]+1)*(faces[1]+1)*(faces[2]+1);
131: PetscInt numFaces = 6;
132: PetscScalar *coords;
133: PetscInt coordSize;
134: PetscMPIInt rank;
135: PetscInt v, vx, vy, vz;
139: if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Must have at least 1 face per side");
140: if ((faces[0] > 1) || (faces[1] > 1) || (faces[2] > 1)) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Currently can't handle more than 1 face per side");
141: PetscMalloc(numVertices*2 * sizeof(PetscReal), &coords);
142: MPI_Comm_rank(((PetscObject) dm)->comm, &rank);
143: if (!rank) {
144: PetscInt f;
146: DMComplexSetChart(dm, 0, numFaces+numVertices);
147: for(f = 0; f < numFaces; ++f) {
148: DMComplexSetConeSize(dm, f, 4);
149: }
150: DMSetUp(dm); /* Allocate space for cones */
151: for(v = 0; v < numFaces+numVertices; ++v) {
152: DMComplexSetLabelValue(dm, "marker", v, 1);
153: }
154: { /* Side 0 (Front) */
155: PetscInt cone[4] = {numFaces+4, numFaces+5, numFaces+7, numFaces+6};
156: DMComplexSetCone(dm, 0, cone);
157: }
158: { /* Side 1 (Back) */
159: PetscInt cone[4] = {numFaces+1, numFaces+0, numFaces+2, numFaces+3};
160: DMComplexSetCone(dm, 1, cone);
161: }
162: { /* Side 2 (Bottom) */
163: PetscInt cone[4] = {numFaces+0, numFaces+1, numFaces+5, numFaces+4};
164: DMComplexSetCone(dm, 2, cone);
165: }
166: { /* Side 3 (Top) */
167: PetscInt cone[4] = {numFaces+6, numFaces+7, numFaces+3, numFaces+2};
168: DMComplexSetCone(dm, 3, cone);
169: }
170: { /* Side 4 (Left) */
171: PetscInt cone[4] = {numFaces+0, numFaces+4, numFaces+6, numFaces+2};
172: DMComplexSetCone(dm, 4, cone);
173: }
174: { /* Side 5 (Right) */
175: PetscInt cone[4] = {numFaces+5, numFaces+1, numFaces+3, numFaces+7};
176: DMComplexSetCone(dm, 5, cone);
177: }
178: }
179: DMComplexSymmetrize(dm);
180: DMComplexStratify(dm);
181: /* Build coordinates */
182: PetscSectionSetChart(mesh->coordSection, numFaces, numFaces + numVertices);
183: for(v = numFaces; v < numFaces+numVertices; ++v) {
184: PetscSectionSetDof(mesh->coordSection, v, 3);
185: }
186: PetscSectionSetUp(mesh->coordSection);
187: PetscSectionGetStorageSize(mesh->coordSection, &coordSize);
188: VecSetSizes(mesh->coordinates, coordSize, PETSC_DETERMINE);
189: VecSetFromOptions(mesh->coordinates);
190: VecGetArray(mesh->coordinates, &coords);
191: for(vz = 0; vz <= faces[2]; ++vz) {
192: for(vy = 0; vy <= faces[1]; ++vy) {
193: for(vx = 0; vx <= faces[0]; ++vx) {
194: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
195: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
196: coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
197: }
198: }
199: }
200: VecRestoreArray(mesh->coordinates, &coords);
201: return(0);
202: }
206: PetscErrorCode DMComplexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool interpolate, DM *dm) {
207: DM boundary;
212: DMCreate(comm, &boundary);
214: DMSetType(boundary, DMCOMPLEX);
215: DMComplexSetDimension(boundary, dim-1);
216: switch(dim) {
217: case 2:
218: {
219: PetscReal lower[2] = {0.0, 0.0};
220: PetscReal upper[2] = {1.0, 1.0};
221: PetscInt edges[2] = {2, 2};
223: DMComplexCreateSquareBoundary(boundary, lower, upper, edges);
224: break;
225: }
226: case 3:
227: {
228: PetscReal lower[3] = {0.0, 0.0, 0.0};
229: PetscReal upper[3] = {1.0, 1.0, 1.0};
230: PetscInt faces[3] = {1, 1, 1};
232: DMComplexCreateCubeBoundary(boundary, lower, upper, faces);
233: break;
234: }
235: default:
236: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %d", dim);
237: }
238: DMComplexGenerate(boundary, PETSC_NULL, interpolate, dm);
239: DMDestroy(&boundary);
240: return(0);
241: }
243: /* External function declarations here */
244: extern PetscErrorCode DMCreateInterpolation_Complex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
245: extern PetscErrorCode DMCreateMatrix_Complex(DM dm, const MatType mtype, Mat *J);
246: extern PetscErrorCode DMRefine_Complex(DM dm, MPI_Comm comm, DM *dmRefined);
247: extern PetscErrorCode DMSetUp_Complex(DM dm);
248: extern PetscErrorCode DMDestroy_Complex(DM dm);
249: extern PetscErrorCode DMView_Complex(DM dm, PetscViewer viewer);
251: EXTERN_C_BEGIN
254: PetscErrorCode DMCreate_Complex(DM dm)
255: {
256: DM_Complex *mesh;
261: PetscNewLog(dm, DM_Complex, &mesh);
262: dm->data = mesh;
264: mesh->dim = 0;
265: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coneSection);
266: mesh->maxConeSize = 0;
267: mesh->cones = PETSC_NULL;
268: mesh->coneOrientations = PETSC_NULL;
269: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->supportSection);
270: mesh->maxSupportSize = 0;
271: mesh->supports = PETSC_NULL;
272: PetscSectionCreate(((PetscObject) dm)->comm, &mesh->coordSection);
273: VecCreate(((PetscObject) dm)->comm, &mesh->coordinates);
274: PetscObjectSetName((PetscObject) mesh->coordinates, "coordinates");
275: mesh->refinementLimit = -1.0;
277: mesh->meetTmpA = PETSC_NULL;
278: mesh->meetTmpB = PETSC_NULL;
279: mesh->joinTmpA = PETSC_NULL;
280: mesh->joinTmpB = PETSC_NULL;
281: mesh->closureTmpA = PETSC_NULL;
282: mesh->closureTmpB = PETSC_NULL;
283: mesh->facesTmp = PETSC_NULL;
285: mesh->labels = PETSC_NULL;
286: mesh->printSetValues = PETSC_FALSE;
288: PetscStrallocpy(VECSTANDARD, &dm->vectype);
289: dm->ops->view = DMView_Complex;
290: dm->ops->setfromoptions = DMSetFromOptions_Complex;
291: dm->ops->setup = DMSetUp_Complex;
292: dm->ops->createglobalvector = PETSC_NULL;
293: dm->ops->createlocalvector = PETSC_NULL;
294: dm->ops->createlocaltoglobalmapping = PETSC_NULL;
295: dm->ops->createlocaltoglobalmappingblock = PETSC_NULL;
296: dm->ops->createfieldis = PETSC_NULL;
297: dm->ops->getcoloring = 0;
298: dm->ops->creatematrix = DMCreateMatrix_Complex;
299: dm->ops->createinterpolation= 0;
300: dm->ops->getaggregates = 0;
301: dm->ops->getinjection = 0;
302: dm->ops->refine = DMRefine_Complex;
303: dm->ops->coarsen = 0;
304: dm->ops->refinehierarchy = 0;
305: dm->ops->coarsenhierarchy = 0;
306: dm->ops->globaltolocalbegin = PETSC_NULL;
307: dm->ops->globaltolocalend = PETSC_NULL;
308: dm->ops->localtoglobalbegin = PETSC_NULL;
309: dm->ops->localtoglobalend = PETSC_NULL;
310: dm->ops->destroy = DMDestroy_Complex;
311: return(0);
312: }
313: EXTERN_C_END
317: /*@
318: DMComplexCreate - Creates a DMComplex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
320: Collective on MPI_Comm
322: Input Parameter:
323: . comm - The communicator for the DMComplex object
325: Output Parameter:
326: . mesh - The DMComplex object
328: Level: beginner
330: .keywords: DMComplex, create
331: @*/
332: PetscErrorCode DMComplexCreate(MPI_Comm comm, DM *mesh)
333: {
338: DMCreate(comm, mesh);
339: DMSetType(*mesh, DMCOMPLEX);
340: return(0);
341: }