Actual source code: stag.c
1: /*
2: Implementation of DMStag, defining dimension-independent functions in the
3: DM API. stag1d.c, stag2d.c, and stag3d.c may include dimension-specific
4: implementations of DM API functions, and other files here contain additional
5: DMStag-specific API functions, as well as internal functions.
6: */
7: #include <petsc/private/dmstagimpl.h>
8: #include <petscsf.h>
10: static PetscErrorCode DMCreateFieldDecomposition_Stag(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
11: {
12: PetscInt f0,f1,f2,f3,dof0,dof1,dof2,dof3,n_entries,k,d,cnt,n_fields,dim;
13: DMStagStencil *stencil0,*stencil1,*stencil2,*stencil3;
15: DMGetDimension(dm,&dim);
16: DMStagGetDOF(dm,&dof0,&dof1,&dof2,&dof3);
17: DMStagGetEntriesPerElement(dm,&n_entries);
19: f0 = 1;
20: f1 = f2 = f3 = 0;
21: if (dim == 1) {
22: f1 = 1;
23: } else if (dim == 2) {
24: f1 = 2;
25: f2 = 1;
26: } else if (dim == 3) {
27: f1 = 3;
28: f2 = 3;
29: f3 = 1;
30: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dimension %" PetscInt_FMT,dim);
32: PetscCalloc1(f0*dof0,&stencil0);
33: PetscCalloc1(f1*dof1,&stencil1);
34: if (dim >= 2) {
35: PetscCalloc1(f2*dof2,&stencil2);
36: }
37: if (dim >= 3) {
38: PetscCalloc1(f3*dof3,&stencil3);
39: }
40: for (k=0; k<f0; ++k) {
41: for (d=0; d<dof0; ++d) {
42: stencil0[dof0*k + d].i = 0; stencil0[dof0*k + d].j = 0; stencil0[dof0*k + d].j = 0;
43: }
44: }
45: for (k=0; k<f1; ++k) {
46: for (d=0; d<dof1; ++d) {
47: stencil1[dof1*k + d].i = 0; stencil1[dof1*k + d].j = 0; stencil1[dof1*k + d].j = 0;
48: }
49: }
50: if (dim >= 2) {
51: for (k=0; k<f2; ++k) {
52: for (d=0; d<dof2; ++d) {
53: stencil2[dof2*k + d].i = 0; stencil2[dof2*k + d].j = 0; stencil2[dof2*k + d].j = 0;
54: }
55: }
56: }
57: if (dim >= 3) {
58: for (k=0; k<f3; ++k) {
59: for (d=0; d<dof3; ++d) {
60: stencil3[dof3*k + d].i = 0; stencil3[dof3*k + d].j = 0; stencil3[dof3*k + d].j = 0;
61: }
62: }
63: }
65: n_fields = 0;
66: if (dof0 != 0) ++n_fields;
67: if (dof1 != 0) ++n_fields;
68: if (dim >=2 && dof2 != 0) ++n_fields;
69: if (dim >=3 && dof3 != 0) ++n_fields;
70: if (len) *len = n_fields;
72: if (islist) {
73: PetscMalloc1(n_fields,islist);
75: if (dim == 1) {
76: /* face, element */
77: for (d=0; d<dof0; ++d) {
78: stencil0[d].loc = DMSTAG_LEFT;
79: stencil0[d].c = d;
80: }
81: for (d=0; d<dof1; ++d) {
82: stencil1[d].loc = DMSTAG_ELEMENT;
83: stencil1[d].c = d;
84: }
85: } else if (dim == 2) {
86: /* vertex, edge(down,left), element */
87: for (d=0; d<dof0; ++d) {
88: stencil0[d].loc = DMSTAG_DOWN_LEFT;
89: stencil0[d].c = d;
90: }
91: /* edge */
92: cnt = 0;
93: for (d=0; d<dof1; ++d) {
94: stencil1[cnt].loc = DMSTAG_DOWN; stencil1[cnt].c = d;
95: ++cnt;
96: }
97: for (d=0; d<dof1; ++d) {
98: stencil1[cnt].loc = DMSTAG_LEFT; stencil1[cnt].c = d;
99: ++cnt;
100: }
101: /* element */
102: for (d=0; d<dof2; ++d) {
103: stencil2[d].loc = DMSTAG_ELEMENT;
104: stencil2[d].c = d;
105: }
106: } else if (dim == 3) {
107: /* vertex, edge(down,left), face(down,left,back), element */
108: for (d=0; d<dof0; ++d) {
109: stencil0[d].loc = DMSTAG_BACK_DOWN_LEFT;
110: stencil0[d].c = d;
111: }
112: /* edges */
113: cnt = 0;
114: for (d=0; d<dof1; ++d) {
115: stencil1[cnt].loc = DMSTAG_BACK_DOWN; stencil1[cnt].c = d;
116: ++cnt;
117: }
118: for (d=0; d<dof1; ++d) {
119: stencil1[cnt].loc = DMSTAG_BACK_LEFT; stencil1[cnt].c = d;
120: ++cnt;
121: }
122: for (d=0; d<dof1; ++d) {
123: stencil1[cnt].loc = DMSTAG_DOWN_LEFT; stencil1[cnt].c = d;
124: ++cnt;
125: }
126: /* faces */
127: cnt = 0;
128: for (d=0; d<dof2; ++d) {
129: stencil2[cnt].loc = DMSTAG_BACK; stencil2[cnt].c = d;
130: ++cnt;
131: }
132: for (d=0; d<dof2; ++d) {
133: stencil2[cnt].loc = DMSTAG_DOWN; stencil2[cnt].c = d;
134: ++cnt;
135: }
136: for (d=0; d<dof2; ++d) {
137: stencil2[cnt].loc = DMSTAG_LEFT; stencil2[cnt].c = d;
138: ++cnt;
139: }
140: /* elements */
141: for (d=0; d<dof3; ++d) {
142: stencil3[d].loc = DMSTAG_ELEMENT;
143: stencil3[d].c = d;
144: }
145: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dimension %" PetscInt_FMT,dim);
147: cnt = 0;
148: if (dof0 != 0) {
149: DMStagCreateISFromStencils(dm,f0*dof0,stencil0,&(*islist)[cnt]);
150: ++cnt;
151: }
152: if (dof1 != 0) {
153: DMStagCreateISFromStencils(dm,f1*dof1,stencil1,&(*islist)[cnt]);
154: ++cnt;
155: }
156: if (dim >= 2 && dof2 != 0) {
157: DMStagCreateISFromStencils(dm,f2*dof2,stencil2,&(*islist)[cnt]);
158: ++cnt;
159: }
160: if (dim >= 3 && dof3 != 0) {
161: DMStagCreateISFromStencils(dm,f3*dof3,stencil3,&(*islist)[cnt]);
162: ++cnt;
163: }
164: }
166: if (namelist) {
167: PetscMalloc1(n_fields,namelist);
168: cnt = 0;
169: if (dim == 1) {
170: if (dof0 != 0) {
171: PetscStrallocpy("vertex",&(*namelist)[cnt]);
172: ++cnt;
173: }
174: if (dof1 != 0) {
175: PetscStrallocpy("element",&(*namelist)[cnt]);
176: ++cnt;
177: }
178: } else if (dim == 2) {
179: if (dof0 != 0) {
180: PetscStrallocpy("vertex",&(*namelist)[cnt]);
181: ++cnt;
182: }
183: if (dof1 != 0) {
184: PetscStrallocpy("face",&(*namelist)[cnt]);
185: ++cnt;
186: }
187: if (dof2 != 0) {
188: PetscStrallocpy("element",&(*namelist)[cnt]);
189: ++cnt;
190: }
191: } else if (dim == 3) {
192: if (dof0 != 0) {
193: PetscStrallocpy("vertex",&(*namelist)[cnt]);
194: ++cnt;
195: }
196: if (dof1 != 0) {
197: PetscStrallocpy("edge",&(*namelist)[cnt]);
198: ++cnt;
199: }
200: if (dof2 != 0) {
201: PetscStrallocpy("face",&(*namelist)[cnt]);
202: ++cnt;
203: }
204: if (dof3 != 0) {
205: PetscStrallocpy("element",&(*namelist)[cnt]);
206: ++cnt;
207: }
208: }
209: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported dimension %" PetscInt_FMT,dim);
210: if (dmlist) {
211: PetscMalloc1(n_fields,dmlist);
212: cnt = 0;
213: if (dof0 != 0) {
214: DMStagCreateCompatibleDMStag(dm,dof0,0,0,0,&(*dmlist)[cnt]);
215: ++cnt;
216: }
217: if (dof1 != 0) {
218: DMStagCreateCompatibleDMStag(dm,0,dof1,0,0,&(*dmlist)[cnt]);
219: ++cnt;
220: }
221: if (dim >= 2 && dof2 != 0) {
222: DMStagCreateCompatibleDMStag(dm,0,0,dof2,0,&(*dmlist)[cnt]);
223: ++cnt;
224: }
225: if (dim >= 3 && dof3 != 0) {
226: DMStagCreateCompatibleDMStag(dm,0,0,0,dof3,&(*dmlist)[cnt]);
227: ++cnt;
228: }
229: }
230: PetscFree(stencil0);
231: PetscFree(stencil1);
232: if (dim >= 2) {
233: PetscFree(stencil2);
234: }
235: if (dim >= 3) {
236: PetscFree(stencil3);
237: }
238: return 0;
239: }
241: static PetscErrorCode DMClone_Stag(DM dm,DM *newdm)
242: {
243: /* Destroy the DM created by generic logic in DMClone() */
244: if (*newdm) {
245: DMDestroy(newdm);
246: }
247: DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);
248: DMSetUp(*newdm);
249: return 0;
250: }
252: static PetscErrorCode DMDestroy_Stag(DM dm)
253: {
254: DM_Stag *stag;
255: PetscInt i;
257: stag = (DM_Stag*)dm->data;
258: for (i=0; i<DMSTAG_MAX_DIM; ++i) {
259: PetscFree(stag->l[i]);
260: }
261: VecScatterDestroy(&stag->gtol);
262: VecScatterDestroy(&stag->ltog_injective);
263: PetscFree(stag->neighbors);
264: PetscFree(stag->locationOffsets);
265: PetscFree(stag->coordinateDMType);
266: PetscFree(stag);
267: return 0;
268: }
270: static PetscErrorCode DMCreateGlobalVector_Stag(DM dm,Vec *vec)
271: {
272: DM_Stag * const stag = (DM_Stag*)dm->data;
275: VecCreateMPI(PetscObjectComm((PetscObject)dm),stag->entries,PETSC_DECIDE,vec);
276: VecSetDM(*vec,dm);
277: /* Could set some ops, as DMDA does */
278: VecSetLocalToGlobalMapping(*vec,dm->ltogmap);
279: return 0;
280: }
282: static PetscErrorCode DMCreateLocalVector_Stag(DM dm,Vec *vec)
283: {
284: DM_Stag * const stag = (DM_Stag*)dm->data;
287: VecCreateSeq(PETSC_COMM_SELF,stag->entriesGhost,vec);
288: VecSetBlockSize(*vec,stag->entriesPerElement);
289: VecSetDM(*vec,dm);
290: return 0;
291: }
293: static PetscErrorCode DMCreateMatrix_Stag(DM dm,Mat *mat)
294: {
295: MatType mat_type;
296: PetscBool is_shell,is_aij;
297: PetscInt dim,entries;
298: ISLocalToGlobalMapping ltogmap;
301: DMGetDimension(dm,&dim);
302: DMGetMatType(dm,&mat_type);
303: DMStagGetEntries(dm,&entries);
304: MatCreate(PetscObjectComm((PetscObject)dm),mat);
305: MatSetSizes(*mat,entries,entries,PETSC_DETERMINE,PETSC_DETERMINE);
306: MatSetType(*mat,mat_type);
307: MatSetUp(*mat);
308: DMGetLocalToGlobalMapping(dm,<ogmap);
309: MatSetLocalToGlobalMapping(*mat,ltogmap,ltogmap);
310: MatSetDM(*mat,dm);
312: /* Compare to similar and perhaps superior logic in DMCreateMatrix_DA, which creates
313: the matrix first and then performs this logic by checking for preallocation functions */
314: PetscStrcmp(mat_type,MATAIJ,&is_aij);
315: if (!is_aij) {
316: PetscStrcmp(mat_type,MATSEQAIJ,&is_aij);
317: } else if (!is_aij) {
318: PetscStrcmp(mat_type,MATMPIAIJ,&is_aij);
319: }
320: PetscStrcmp(mat_type,MATSHELL,&is_shell);
321: if (is_aij) {
322: Mat preallocator;
323: PetscInt m,n;
324: const PetscBool fill_with_zeros = PETSC_FALSE;
326: MatCreate(PetscObjectComm((PetscObject)dm),&preallocator);
327: MatSetType(preallocator,MATPREALLOCATOR);
328: MatGetLocalSize(*mat,&m,&n);
329: MatSetSizes(preallocator,m,n,PETSC_DECIDE,PETSC_DECIDE);
330: MatSetLocalToGlobalMapping(preallocator,ltogmap,ltogmap);
331: MatSetUp(preallocator);
332: switch (dim) {
333: case 1:
334: DMCreateMatrix_Stag_1D_AIJ_Assemble(dm,preallocator);
335: break;
336: case 2:
337: DMCreateMatrix_Stag_2D_AIJ_Assemble(dm,preallocator);
338: break;
339: case 3:
340: DMCreateMatrix_Stag_3D_AIJ_Assemble(dm,preallocator);
341: break;
342: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %" PetscInt_FMT,dim);
343: }
344: MatPreallocatorPreallocate(preallocator,fill_with_zeros,*mat);
345: MatDestroy(&preallocator);
347: if (!dm->prealloc_only) {
348: switch (dim) {
349: case 1:
350: DMCreateMatrix_Stag_1D_AIJ_Assemble(dm,*mat);
351: break;
352: case 2:
353: DMCreateMatrix_Stag_2D_AIJ_Assemble(dm,*mat);
354: break;
355: case 3:
356: DMCreateMatrix_Stag_3D_AIJ_Assemble(dm,*mat);
357: break;
358: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %" PetscInt_FMT,dim);
359: }
360: }
361: /* Note: GPU-related logic, e.g. at the end of DMCreateMatrix_DA_1d_MPIAIJ, is not included here
362: but might be desirable */
363: } else if (is_shell) {
364: /* nothing more to do */
365: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for Mattype %s",mat_type);
366: return 0;
367: }
369: static PetscErrorCode DMGetCompatibility_Stag(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
370: {
371: const DM_Stag * const stag = (DM_Stag*)dm->data;
372: const DM_Stag * const stag2 = (DM_Stag*)dm2->data;
373: PetscInt dim,dim2,i;
374: MPI_Comm comm;
375: PetscMPIInt sameComm;
376: DMType type2;
377: PetscBool sameType;
379: DMGetType(dm2,&type2);
380: PetscStrcmp(DMSTAG,type2,&sameType);
381: if (!sameType) {
382: PetscInfo((PetscObject)dm,"DMStag compatibility check not implemented with DM of type %s\n",type2);
383: *set = PETSC_FALSE;
384: return 0;
385: }
387: PetscObjectGetComm((PetscObject)dm,&comm);
388: MPI_Comm_compare(comm,PetscObjectComm((PetscObject)dm2),&sameComm);
389: if (sameComm != MPI_IDENT) {
390: PetscInfo((PetscObject)dm,"DMStag objects have different communicators: %d != %d\n",comm,PetscObjectComm((PetscObject)dm2));
391: *set = PETSC_FALSE;
392: return 0;
393: }
394: DMGetDimension(dm ,&dim);
395: DMGetDimension(dm2,&dim2);
396: if (dim != dim2) {
397: PetscInfo((PetscObject)dm,"DMStag objects have different dimensions");
398: *set = PETSC_TRUE;
399: *compatible = PETSC_FALSE;
400: return 0;
401: }
402: for (i=0; i<dim; ++i) {
403: if (stag->N[i] != stag2->N[i]) {
404: PetscInfo((PetscObject)dm,"DMStag objects have different global numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
405: *set = PETSC_TRUE;
406: *compatible = PETSC_FALSE;
407: return 0;
408: }
409: if (stag->n[i] != stag2->n[i]) {
410: PetscInfo((PetscObject)dm,"DMStag objects have different local numbers of elements in dimension %D: %D != %D\n",i,stag->n[i],stag2->n[i]);
411: *set = PETSC_TRUE;
412: *compatible = PETSC_FALSE;
413: return 0;
414: }
415: if (stag->boundaryType[i] != stag2->boundaryType[i]) {
416: PetscInfo((PetscObject)dm,"DMStag objects have different boundary types in dimension %d: %s != %s\n",i,stag->boundaryType[i],stag2->boundaryType[i]);
417: *set = PETSC_TRUE;
418: *compatible = PETSC_FALSE;
419: return 0;
420: }
421: }
422: /* Note: we include stencil type and width in the notion of compatibility, as this affects
423: the "atlas" (local subdomains). This might be irritating in legitimate cases
424: of wanting to transfer between two other-wise compatible DMs with different
425: stencil characteristics. */
426: if (stag->stencilType != stag2->stencilType) {
427: PetscInfo((PetscObject)dm,"DMStag objects have different ghost stencil types: %s != %s\n",DMStagStencilTypes[stag->stencilType],DMStagStencilTypes[stag2->stencilType]);
428: *set = PETSC_TRUE;
429: *compatible = PETSC_FALSE;
430: return 0;
431: }
432: if (stag->stencilWidth != stag2->stencilWidth) {
433: PetscInfo((PetscObject)dm,"DMStag objects have different ghost stencil widths: %D != %D\n",stag->stencilWidth,stag->stencilWidth);
434: *set = PETSC_TRUE;
435: *compatible = PETSC_FALSE;
436: return 0;
437: }
438: *set = PETSC_TRUE;
439: *compatible = PETSC_TRUE;
440: return 0;
441: }
443: /*
444: Note there are several orderings in play here.
445: In all cases, non-element dof are associated with the element that they are below/left/behind, and the order in 2D proceeds vertex/bottom edge/left edge/element (with all dof on each together).
446: Also in all cases, only subdomains which are the last in their dimension have partial elements.
448: 1) "Natural" Ordering (not used). Number adding each full or partial (on the right or top) element, starting at the bottom left (i=0,j=0) and proceeding across the entire domain, row by row to get a global numbering.
449: 2) Global ("PETSc") ordering. The same as natural, but restricted to each domain. So, traverse all elements (again starting at the bottom left and going row-by-row) on rank 0, then continue numbering with rank 1, and so on.
450: 3) Local ordering. Including ghost elements (both interior and on the right/top/front to complete partial elements), use the same convention to create a local numbering.
451: */
453: static PetscErrorCode DMLocalToGlobalBegin_Stag(DM dm,Vec l,InsertMode mode,Vec g)
454: {
455: DM_Stag * const stag = (DM_Stag*)dm->data;
457: if (mode == ADD_VALUES) {
458: VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE);
459: } else if (mode == INSERT_VALUES) {
460: if (stag->ltog_injective) {
461: VecScatterBegin(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
462: } else {
463: VecScatterBegin(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
464: }
465: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
466: return 0;
467: }
469: static PetscErrorCode DMLocalToGlobalEnd_Stag(DM dm,Vec l,InsertMode mode,Vec g)
470: {
471: DM_Stag * const stag = (DM_Stag*)dm->data;
473: if (mode == ADD_VALUES) {
474: VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE);
475: } else if (mode == INSERT_VALUES) {
476: if (stag->ltog_injective) {
477: VecScatterEnd(stag->ltog_injective,l,g,mode,SCATTER_FORWARD);
478: } else {
479: VecScatterEnd(stag->gtol,l,g,mode,SCATTER_REVERSE_LOCAL);
480: }
481: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported InsertMode");
482: return 0;
483: }
485: static PetscErrorCode DMGlobalToLocalBegin_Stag(DM dm,Vec g,InsertMode mode,Vec l)
486: {
487: DM_Stag * const stag = (DM_Stag*)dm->data;
489: VecScatterBegin(stag->gtol,g,l,mode,SCATTER_FORWARD);
490: return 0;
491: }
493: static PetscErrorCode DMGlobalToLocalEnd_Stag(DM dm,Vec g,InsertMode mode,Vec l)
494: {
495: DM_Stag * const stag = (DM_Stag*)dm->data;
497: VecScatterEnd(stag->gtol,g,l,mode,SCATTER_FORWARD);
498: return 0;
499: }
501: /*
502: If a stratum is active (non-zero dof), make it active in the coordinate DM.
503: */
504: static PetscErrorCode DMCreateCoordinateDM_Stag(DM dm,DM *dmc)
505: {
506: DM_Stag * const stag = (DM_Stag*)dm->data;
507: PetscInt dim;
508: PetscBool isstag,isproduct;
513: DMGetDimension(dm,&dim);
514: PetscStrcmp(stag->coordinateDMType,DMSTAG,&isstag);
515: PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&isproduct);
516: if (isstag) {
517: PetscCall(DMStagCreateCompatibleDMStag(dm,
518: stag->dof[0] > 0 ? dim : 0,
519: stag->dof[1] > 0 ? dim : 0,
520: stag->dof[2] > 0 ? dim : 0,
521: stag->dof[3] > 0 ? dim : 0,
522: dmc));
523: } else if (isproduct) {
524: DMCreate(PETSC_COMM_WORLD,dmc);
525: DMSetType(*dmc,DMPRODUCT);
526: DMSetDimension(*dmc,dim);
527: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported coordinate DM type %s",stag->coordinateDMType);
528: return 0;
529: }
531: static PetscErrorCode DMGetNeighbors_Stag(DM dm,PetscInt *nRanks,const PetscMPIInt *ranks[])
532: {
533: DM_Stag * const stag = (DM_Stag*)dm->data;
534: PetscInt dim;
536: DMGetDimension(dm,&dim);
537: switch (dim) {
538: case 1: *nRanks = 3; break;
539: case 2: *nRanks = 9; break;
540: case 3: *nRanks = 27; break;
541: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Get neighbors not implemented for dim = %D",dim);
542: }
543: *ranks = stag->neighbors;
544: return 0;
545: }
547: static PetscErrorCode DMView_Stag(DM dm,PetscViewer viewer)
548: {
549: DM_Stag * const stag = (DM_Stag*)dm->data;
550: PetscBool isascii,viewAllRanks;
551: PetscMPIInt rank,size;
552: PetscInt dim,maxRanksToView,i;
554: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
555: MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
556: DMGetDimension(dm,&dim);
557: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
558: if (isascii) {
559: PetscViewerASCIIPrintf(viewer,"Dimension: %D\n",dim);
560: switch (dim) {
561: case 1:
562: PetscViewerASCIIPrintf(viewer,"Global size: %D\n",stag->N[0]);
563: break;
564: case 2:
565: PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D\n",stag->N[0],stag->N[1]);
566: PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1]);
567: break;
568: case 3:
569: PetscViewerASCIIPrintf(viewer,"Global sizes: %D x %D x %D\n",stag->N[0],stag->N[1],stag->N[2]);
570: PetscViewerASCIIPrintf(viewer,"Parallel decomposition: %D x %D x %D ranks\n",stag->nRanks[0],stag->nRanks[1],stag->nRanks[2]);
571: break;
572: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
573: }
574: PetscViewerASCIIPrintf(viewer,"Boundary ghosting:");
575: for (i=0; i<dim; ++i) {
576: PetscViewerASCIIPrintf(viewer," %s",DMBoundaryTypes[stag->boundaryType[i]]);
577: }
578: PetscViewerASCIIPrintf(viewer,"\n");
579: PetscViewerASCIIPrintf(viewer,"Elementwise ghost stencil: %s",DMStagStencilTypes[stag->stencilType]);
580: if (stag->stencilType != DMSTAG_STENCIL_NONE) {
581: PetscViewerASCIIPrintf(viewer,", width %D\n",stag->stencilWidth);
582: } else {
583: PetscViewerASCIIPrintf(viewer,"\n");
584: }
585: PetscViewerASCIIPrintf(viewer,"%D DOF per vertex (0D)\n",stag->dof[0]);
586: if (dim == 3) {
587: PetscViewerASCIIPrintf(viewer,"%D DOF per edge (1D)\n",stag->dof[1]);
588: }
589: if (dim > 1) {
590: PetscViewerASCIIPrintf(viewer,"%D DOF per face (%DD)\n",stag->dof[dim-1],dim-1);
591: }
592: PetscViewerASCIIPrintf(viewer,"%D DOF per element (%DD)\n",stag->dof[dim],dim);
593: if (dm->coordinateDM) {
594: PetscViewerASCIIPrintf(viewer,"Has coordinate DM\n");
595: }
596: maxRanksToView = 16;
597: viewAllRanks = (PetscBool)(size <= maxRanksToView);
598: if (viewAllRanks) {
599: PetscViewerASCIIPushSynchronized(viewer);
600: switch (dim) {
601: case 1:
602: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D (%D with ghosts)\n",rank,stag->n[0],stag->nGhost[0]);
603: break;
604: case 2:
605: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d)\n",rank,stag->rank[0],stag->rank[1]);
606: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D (%D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->nGhost[0],stag->nGhost[1]);
607: break;
608: case 3:
609: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Rank coordinates (%d,%d,%d)\n",rank,stag->rank[0],stag->rank[1],stag->rank[2]);
610: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local elements : %D x %D x %D (%D x %D x %D with ghosts)\n",rank,stag->n[0],stag->n[1],stag->n[2],stag->nGhost[0],stag->nGhost[1],stag->nGhost[2]);
611: break;
612: default: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"not implemented for dim==%D",dim);
613: }
614: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local native entries: %d\n",rank,stag->entries);
615: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local entries total : %d\n",rank,stag->entriesGhost);
616: PetscViewerFlush(viewer);
617: PetscViewerASCIIPopSynchronized(viewer);
618: } else {
619: PetscViewerASCIIPrintf(viewer,"(Per-rank information omitted since >%D ranks used)\n",maxRanksToView);
620: }
621: }
622: return 0;
623: }
625: static PetscErrorCode DMSetFromOptions_Stag(PetscOptionItems *PetscOptionsObject,DM dm)
626: {
627: DM_Stag * const stag = (DM_Stag*)dm->data;
628: PetscInt dim;
630: DMGetDimension(dm,&dim);
631: PetscOptionsHead(PetscOptionsObject,"DMStag Options");
632: PetscOptionsInt("-stag_grid_x","Number of grid points in x direction","DMStagSetGlobalSizes",stag->N[0],&stag->N[0],NULL);
633: if (dim > 1) PetscOptionsInt("-stag_grid_y","Number of grid points in y direction","DMStagSetGlobalSizes",stag->N[1],&stag->N[1],NULL);
634: if (dim > 2) PetscOptionsInt("-stag_grid_z","Number of grid points in z direction","DMStagSetGlobalSizes",stag->N[2],&stag->N[2],NULL);
635: PetscOptionsInt("-stag_ranks_x","Number of ranks in x direction","DMStagSetNumRanks",stag->nRanks[0],&stag->nRanks[0],NULL);
636: if (dim > 1) PetscOptionsInt("-stag_ranks_y","Number of ranks in y direction","DMStagSetNumRanks",stag->nRanks[1],&stag->nRanks[1],NULL);
637: if (dim > 2) PetscOptionsInt("-stag_ranks_z","Number of ranks in z direction","DMStagSetNumRanks",stag->nRanks[2],&stag->nRanks[2],NULL);
638: PetscOptionsInt("-stag_stencil_width","Elementwise stencil width","DMStagSetStencilWidth",stag->stencilWidth,&stag->stencilWidth,NULL);
639: PetscOptionsEnum("-stag_stencil_type","Elementwise stencil stype","DMStagSetStencilType",DMStagStencilTypes,(PetscEnum)stag->stencilType,(PetscEnum*)&stag->stencilType,NULL);
640: PetscOptionsEnum("-stag_boundary_type_x","Treatment of (physical) boundaries in x direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[0],(PetscEnum*)&stag->boundaryType[0],NULL);
641: PetscOptionsEnum("-stag_boundary_type_y","Treatment of (physical) boundaries in y direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[1],(PetscEnum*)&stag->boundaryType[1],NULL);
642: PetscOptionsEnum("-stag_boundary_type_z","Treatment of (physical) boundaries in z direction","DMStagSetBoundaryTypes",DMBoundaryTypes,(PetscEnum)stag->boundaryType[2],(PetscEnum*)&stag->boundaryType[2],NULL);
643: PetscOptionsInt("-stag_dof_0","Number of dof per 0-cell (vertex)","DMStagSetDOF",stag->dof[0],&stag->dof[0],NULL);
644: PetscOptionsInt("-stag_dof_1","Number of dof per 1-cell (element in 1D, face in 2D, edge in 3D)","DMStagSetDOF",stag->dof[1],&stag->dof[1],NULL);
645: PetscOptionsInt("-stag_dof_2","Number of dof per 2-cell (element in 2D, face in 3D)","DMStagSetDOF",stag->dof[2],&stag->dof[2],NULL);
646: PetscOptionsInt("-stag_dof_3","Number of dof per 3-cell (element in 3D)","DMStagSetDOF",stag->dof[3],&stag->dof[3],NULL);
647: PetscOptionsTail();
648: return 0;
649: }
651: /*MC
652: DMSTAG = "stag" - A DM object representing a "staggered grid" or a structured cell complex.
654: This implementation parallels the DMDA implementation in many ways, but allows degrees of freedom
655: to be associated with all "strata" in a logically-rectangular grid.
657: Each stratum can be characterized by the dimension of the entities ("points", to borrow the DMPLEX
658: terminology), from 0- to 3-dimensional.
660: In some cases this numbering is used directly, for example with DMStagGetDOF().
661: To allow easier reading and to some extent more similar code between different-dimensional implementations
662: of the same problem, we associate canonical names for each type of point, for each dimension of DMStag.
664: 1-dimensional DMStag objects have vertices (0D) and elements (1D).
666: 2-dimensional DMStag objects have vertices (0D), faces (1D), and elements (2D).
668: 3-dimensional DMStag objects have vertices (0D), edges (1D), faces (2D), and elements (3D).
670: This naming is reflected when viewing a DMStag object with DMView() , and in forming
671: convenient options prefixes when creating a decomposition with DMCreateFieldDecomposition().
673: Level: beginner
675: .seealso: DM, DMPRODUCT, DMDA, DMPLEX, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMType, DMCreate(), DMSetType()
676: M*/
678: PETSC_EXTERN PetscErrorCode DMCreate_Stag(DM dm)
679: {
680: DM_Stag *stag;
681: PetscInt i,dim;
684: PetscNewLog(dm,&stag);
685: dm->data = stag;
687: stag->gtol = NULL;
688: stag->ltog_injective = NULL;
689: for (i=0; i<DMSTAG_MAX_STRATA; ++i) stag->dof[i] = 0;
690: for (i=0; i<DMSTAG_MAX_DIM; ++i) stag->l[i] = NULL;
691: stag->stencilType = DMSTAG_STENCIL_NONE;
692: stag->stencilWidth = 0;
693: for (i=0; i<DMSTAG_MAX_DIM; ++i) stag->nRanks[i] = -1;
694: stag->coordinateDMType = NULL;
696: DMGetDimension(dm,&dim);
699: PetscMemzero(dm->ops,sizeof(*(dm->ops)));
700: dm->ops->createcoordinatedm = DMCreateCoordinateDM_Stag;
701: dm->ops->createglobalvector = DMCreateGlobalVector_Stag;
702: dm->ops->createinterpolation = NULL;
703: dm->ops->createlocalvector = DMCreateLocalVector_Stag;
704: dm->ops->creatematrix = DMCreateMatrix_Stag;
705: dm->ops->destroy = DMDestroy_Stag;
706: dm->ops->getneighbors = DMGetNeighbors_Stag;
707: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Stag;
708: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Stag;
709: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Stag;
710: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Stag;
711: dm->ops->setfromoptions = DMSetFromOptions_Stag;
712: switch (dim) {
713: case 1: dm->ops->setup = DMSetUp_Stag_1d; break;
714: case 2: dm->ops->setup = DMSetUp_Stag_2d; break;
715: case 3: dm->ops->setup = DMSetUp_Stag_3d; break;
716: default : SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %d",dim);
717: }
718: dm->ops->clone = DMClone_Stag;
719: dm->ops->view = DMView_Stag;
720: dm->ops->getcompatibility = DMGetCompatibility_Stag;
721: dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Stag;
722: return 0;
723: }