1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
2
3 /*@
4 DMDASetSizes - Sets the number of grid points in the three dimensional directions
5
6 Logically Collective on da
7
8 Input Parameters:
9 + da - the DMDA
10 . M - the global X size
11 . N - the global Y size
12 - P - the global Z size
13
14 Level: intermediate
15
16 Developer Notes:
17 Since the dimension may not yet have been set the code cannot error check for non-positive Y and Z number of grid points
18
19 .seealso: PetscSplitOwnership()
20 @*/
DMDASetSizes(DM da,PetscInt M,PetscInt N,PetscInt P)21 PetscErrorCode DMDASetSizes(DM da, PetscInt M, PetscInt N, PetscInt P)
22 {
23 DM_DA *dd = (DM_DA*)da->data;
24
25 PetscFunctionBegin;
26 PetscValidHeaderSpecificType(da, DM_CLASSID, 1,DMDA);
27 PetscValidLogicalCollectiveInt(da,M,2);
28 PetscValidLogicalCollectiveInt(da,N,3);
29 PetscValidLogicalCollectiveInt(da,P,4);
30 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
31 if (M < 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in X direction must be positive");
32 if (N < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Y direction must be positive");
33 if (P < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Number of grid points in Z direction must be positive");
34
35 dd->M = M;
36 dd->N = N;
37 dd->P = P;
38 PetscFunctionReturn(0);
39 }
40
41 /*@
42 DMDASetNumProcs - Sets the number of processes in each dimension
43
44 Logically Collective on da
45
46 Input Parameters:
47 + da - the DMDA
48 . m - the number of X procs (or PETSC_DECIDE)
49 . n - the number of Y procs (or PETSC_DECIDE)
50 - p - the number of Z procs (or PETSC_DECIDE)
51
52 Level: intermediate
53
54 .seealso: DMDASetSizes(), PetscSplitOwnership()
55 @*/
DMDASetNumProcs(DM da,PetscInt m,PetscInt n,PetscInt p)56 PetscErrorCode DMDASetNumProcs(DM da, PetscInt m, PetscInt n, PetscInt p)
57 {
58 DM_DA *dd = (DM_DA*)da->data;
59 PetscErrorCode ierr;
60
61 PetscFunctionBegin;
62 PetscValidHeaderSpecificType(da, DM_CLASSID, 1,DMDA);
63 PetscValidLogicalCollectiveInt(da,m,2);
64 PetscValidLogicalCollectiveInt(da,n,3);
65 PetscValidLogicalCollectiveInt(da,p,4);
66 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
67 dd->m = m;
68 dd->n = n;
69 dd->p = p;
70 if (da->dim == 2) {
71 PetscMPIInt size;
72 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr);
73 if ((dd->m > 0) && (dd->n < 0)) {
74 dd->n = size/dd->m;
75 if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in X direction not divisible into comm size %d",m,size);
76 }
77 if ((dd->n > 0) && (dd->m < 0)) {
78 dd->m = size/dd->n;
79 if (dd->n*dd->m != size) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"%D processes in Y direction not divisible into comm size %d",n,size);
80 }
81 }
82 PetscFunctionReturn(0);
83 }
84
85 /*@
86 DMDASetBoundaryType - Sets the type of ghost nodes on domain boundaries.
87
88 Not collective
89
90 Input Parameter:
91 + da - The DMDA
92 - bx,by,bz - One of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC
93
94 Level: intermediate
95
96 .seealso: DMDACreate(), DMDestroy(), DMDA, DMBoundaryType
97 @*/
DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)98 PetscErrorCode DMDASetBoundaryType(DM da,DMBoundaryType bx,DMBoundaryType by,DMBoundaryType bz)
99 {
100 DM_DA *dd = (DM_DA*)da->data;
101
102 PetscFunctionBegin;
103 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
104 PetscValidLogicalCollectiveEnum(da,bx,2);
105 PetscValidLogicalCollectiveEnum(da,by,3);
106 PetscValidLogicalCollectiveEnum(da,bz,4);
107 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
108 dd->bx = bx;
109 dd->by = by;
110 dd->bz = bz;
111 PetscFunctionReturn(0);
112 }
113
114 /*@
115 DMDASetDof - Sets the number of degrees of freedom per vertex
116
117 Not collective
118
119 Input Parameters:
120 + da - The DMDA
121 - dof - Number of degrees of freedom
122
123 Level: intermediate
124
125 .seealso: DMDAGetDof(), DMDACreate(), DMDestroy(), DMDA
126 @*/
DMDASetDof(DM da,PetscInt dof)127 PetscErrorCode DMDASetDof(DM da, PetscInt dof)
128 {
129 DM_DA *dd = (DM_DA*)da->data;
130
131 PetscFunctionBegin;
132 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
133 PetscValidLogicalCollectiveInt(da,dof,2);
134 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
135 dd->w = dof;
136 da->bs = dof;
137 PetscFunctionReturn(0);
138 }
139
140 /*@
141 DMDAGetDof - Gets the number of degrees of freedom per vertex
142
143 Not collective
144
145 Input Parameter:
146 . da - The DMDA
147
148 Output Parameter:
149 . dof - Number of degrees of freedom
150
151 Level: intermediate
152
153 .seealso: DMDASetDof(), DMDACreate(), DMDestroy(), DMDA
154 @*/
DMDAGetDof(DM da,PetscInt * dof)155 PetscErrorCode DMDAGetDof(DM da, PetscInt *dof)
156 {
157 DM_DA *dd = (DM_DA *) da->data;
158
159 PetscFunctionBegin;
160 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
161 PetscValidPointer(dof,2);
162 *dof = dd->w;
163 PetscFunctionReturn(0);
164 }
165
166 /*@
167 DMDAGetOverlap - Gets the size of the per-processor overlap.
168
169 Not collective
170
171 Input Parameters:
172 . da - The DMDA
173
174 Output Parameters:
175 + x - Overlap in the x direction
176 . y - Overlap in the y direction
177 - z - Overlap in the z direction
178
179 Level: intermediate
180
181 .seealso: DMDACreateDomainDecomposition(), DMDASetOverlap(), DMDA
182 @*/
DMDAGetOverlap(DM da,PetscInt * x,PetscInt * y,PetscInt * z)183 PetscErrorCode DMDAGetOverlap(DM da,PetscInt *x,PetscInt *y,PetscInt *z)
184 {
185 DM_DA *dd = (DM_DA*)da->data;
186
187 PetscFunctionBegin;
188 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
189 if (x) *x = dd->xol;
190 if (y) *y = dd->yol;
191 if (z) *z = dd->zol;
192 PetscFunctionReturn(0);
193 }
194
195 /*@
196 DMDASetOverlap - Sets the size of the per-processor overlap.
197
198 Not collective
199
200 Input Parameters:
201 + da - The DMDA
202 . x - Overlap in the x direction
203 . y - Overlap in the y direction
204 - z - Overlap in the z direction
205
206 Level: intermediate
207
208 .seealso: DMDACreateDomainDecomposition(), DMDAGetOverlap(), DMDA
209 @*/
DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)210 PetscErrorCode DMDASetOverlap(DM da,PetscInt x,PetscInt y,PetscInt z)
211 {
212 DM_DA *dd = (DM_DA*)da->data;
213
214 PetscFunctionBegin;
215 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
216 PetscValidLogicalCollectiveInt(da,x,2);
217 PetscValidLogicalCollectiveInt(da,y,3);
218 PetscValidLogicalCollectiveInt(da,z,4);
219 dd->xol = x;
220 dd->yol = y;
221 dd->zol = z;
222 PetscFunctionReturn(0);
223 }
224
225
226 /*@
227 DMDAGetNumLocalSubDomains - Gets the number of local subdomains created upon decomposition.
228
229 Not collective
230
231 Input Parameters:
232 . da - The DMDA
233
234 Output Parameters:
235 . Nsub - Number of local subdomains created upon decomposition
236
237 Level: intermediate
238
239 .seealso: DMDACreateDomainDecomposition(), DMDASetNumLocalSubDomains(), DMDA
240 @*/
DMDAGetNumLocalSubDomains(DM da,PetscInt * Nsub)241 PetscErrorCode DMDAGetNumLocalSubDomains(DM da,PetscInt *Nsub)
242 {
243 DM_DA *dd = (DM_DA*)da->data;
244
245 PetscFunctionBegin;
246 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
247 if (Nsub) *Nsub = dd->Nsub;
248 PetscFunctionReturn(0);
249 }
250
251 /*@
252 DMDASetNumLocalSubDomains - Sets the number of local subdomains created upon decomposition.
253
254 Not collective
255
256 Input Parameters:
257 + da - The DMDA
258 - Nsub - The number of local subdomains requested
259
260 Level: intermediate
261
262 .seealso: DMDACreateDomainDecomposition(), DMDAGetNumLocalSubDomains(), DMDA
263 @*/
DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)264 PetscErrorCode DMDASetNumLocalSubDomains(DM da,PetscInt Nsub)
265 {
266 DM_DA *dd = (DM_DA*)da->data;
267
268 PetscFunctionBegin;
269 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
270 PetscValidLogicalCollectiveInt(da,Nsub,2);
271 dd->Nsub = Nsub;
272 PetscFunctionReturn(0);
273 }
274
275 /*@
276 DMDASetOffset - Sets the index offset of the DA.
277
278 Collective on DA
279
280 Input Parameter:
281 + da - The DMDA
282 . xo - The offset in the x direction
283 . yo - The offset in the y direction
284 - zo - The offset in the z direction
285
286 Level: intermediate
287
288 Notes:
289 This is used primarily to overlap a computation on a local DA with that on a global DA without
290 changing boundary conditions or subdomain features that depend upon the global offsets.
291
292 .seealso: DMDAGetOffset(), DMDAVecGetArray()
293 @*/
DMDASetOffset(DM da,PetscInt xo,PetscInt yo,PetscInt zo,PetscInt Mo,PetscInt No,PetscInt Po)294 PetscErrorCode DMDASetOffset(DM da, PetscInt xo, PetscInt yo, PetscInt zo, PetscInt Mo, PetscInt No, PetscInt Po)
295 {
296 PetscErrorCode ierr;
297 DM_DA *dd = (DM_DA*)da->data;
298
299 PetscFunctionBegin;
300 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
301 PetscValidLogicalCollectiveInt(da,xo,2);
302 PetscValidLogicalCollectiveInt(da,yo,3);
303 PetscValidLogicalCollectiveInt(da,zo,4);
304 PetscValidLogicalCollectiveInt(da,Mo,5);
305 PetscValidLogicalCollectiveInt(da,No,6);
306 PetscValidLogicalCollectiveInt(da,Po,7);
307 dd->xo = xo;
308 dd->yo = yo;
309 dd->zo = zo;
310 dd->Mo = Mo;
311 dd->No = No;
312 dd->Po = Po;
313
314 if (da->coordinateDM) {
315 ierr = DMDASetOffset(da->coordinateDM,xo,yo,zo,Mo,No,Po);CHKERRQ(ierr);
316 }
317 PetscFunctionReturn(0);
318 }
319
320 /*@
321 DMDAGetOffset - Gets the index offset of the DA.
322
323 Not collective
324
325 Input Parameter:
326 . da - The DMDA
327
328 Output Parameters:
329 + xo - The offset in the x direction
330 . yo - The offset in the y direction
331 . zo - The offset in the z direction
332 . Mo - The global size in the x direction
333 . No - The global size in the y direction
334 - Po - The global size in the z direction
335
336 Level: intermediate
337
338 .seealso: DMDASetOffset(), DMDAVecGetArray()
339 @*/
DMDAGetOffset(DM da,PetscInt * xo,PetscInt * yo,PetscInt * zo,PetscInt * Mo,PetscInt * No,PetscInt * Po)340 PetscErrorCode DMDAGetOffset(DM da,PetscInt *xo,PetscInt *yo,PetscInt *zo,PetscInt *Mo,PetscInt *No,PetscInt *Po)
341 {
342 DM_DA *dd = (DM_DA*)da->data;
343
344 PetscFunctionBegin;
345 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
346 if (xo) *xo = dd->xo;
347 if (yo) *yo = dd->yo;
348 if (zo) *zo = dd->zo;
349 if (Mo) *Mo = dd->Mo;
350 if (No) *No = dd->No;
351 if (Po) *Po = dd->Po;
352 PetscFunctionReturn(0);
353 }
354
355 /*@
356 DMDAGetNonOverlappingRegion - Gets the indices of the nonoverlapping region of a subdomain DM.
357
358 Not collective
359
360 Input Parameter:
361 . da - The DMDA
362
363 Output Parameters:
364 + xs - The start of the region in x
365 . ys - The start of the region in y
366 . zs - The start of the region in z
367 . xs - The size of the region in x
368 . ys - The size of the region in y
369 - zs - The size of the region in z
370
371 Level: intermediate
372
373 .seealso: DMDAGetOffset(), DMDAVecGetArray()
374 @*/
DMDAGetNonOverlappingRegion(DM da,PetscInt * xs,PetscInt * ys,PetscInt * zs,PetscInt * xm,PetscInt * ym,PetscInt * zm)375 PetscErrorCode DMDAGetNonOverlappingRegion(DM da, PetscInt *xs, PetscInt *ys, PetscInt *zs, PetscInt *xm, PetscInt *ym, PetscInt *zm)
376 {
377 DM_DA *dd = (DM_DA*)da->data;
378
379 PetscFunctionBegin;
380 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
381 if (xs) *xs = dd->nonxs;
382 if (ys) *ys = dd->nonys;
383 if (zs) *zs = dd->nonzs;
384 if (xm) *xm = dd->nonxm;
385 if (ym) *ym = dd->nonym;
386 if (zm) *zm = dd->nonzm;
387 PetscFunctionReturn(0);
388 }
389
390
391 /*@
392 DMDASetNonOverlappingRegion - Sets the indices of the nonoverlapping region of a subdomain DM.
393
394 Collective on DA
395
396 Input Parameter:
397 + da - The DMDA
398 . xs - The start of the region in x
399 . ys - The start of the region in y
400 . zs - The start of the region in z
401 . xs - The size of the region in x
402 . ys - The size of the region in y
403 - zs - The size of the region in z
404
405 Level: intermediate
406
407 .seealso: DMDAGetOffset(), DMDAVecGetArray()
408 @*/
DMDASetNonOverlappingRegion(DM da,PetscInt xs,PetscInt ys,PetscInt zs,PetscInt xm,PetscInt ym,PetscInt zm)409 PetscErrorCode DMDASetNonOverlappingRegion(DM da, PetscInt xs, PetscInt ys, PetscInt zs, PetscInt xm, PetscInt ym, PetscInt zm)
410 {
411 DM_DA *dd = (DM_DA*)da->data;
412
413 PetscFunctionBegin;
414 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
415 PetscValidLogicalCollectiveInt(da,xs,2);
416 PetscValidLogicalCollectiveInt(da,ys,3);
417 PetscValidLogicalCollectiveInt(da,zs,4);
418 PetscValidLogicalCollectiveInt(da,xm,5);
419 PetscValidLogicalCollectiveInt(da,ym,6);
420 PetscValidLogicalCollectiveInt(da,zm,7);
421 dd->nonxs = xs;
422 dd->nonys = ys;
423 dd->nonzs = zs;
424 dd->nonxm = xm;
425 dd->nonym = ym;
426 dd->nonzm = zm;
427
428 PetscFunctionReturn(0);
429 }
430
431 /*@
432 DMDASetStencilType - Sets the type of the communication stencil
433
434 Logically Collective on da
435
436 Input Parameter:
437 + da - The DMDA
438 - stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
439
440 Level: intermediate
441
442 .seealso: DMDACreate(), DMDestroy(), DMDA
443 @*/
DMDASetStencilType(DM da,DMDAStencilType stype)444 PetscErrorCode DMDASetStencilType(DM da, DMDAStencilType stype)
445 {
446 DM_DA *dd = (DM_DA*)da->data;
447
448 PetscFunctionBegin;
449 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
450 PetscValidLogicalCollectiveEnum(da,stype,2);
451 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
452 dd->stencil_type = stype;
453 PetscFunctionReturn(0);
454 }
455
456 /*@
457 DMDAGetStencilType - Gets the type of the communication stencil
458
459 Not collective
460
461 Input Parameter:
462 . da - The DMDA
463
464 Output Parameter:
465 . stype - The stencil type, use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR.
466
467 Level: intermediate
468
469 .seealso: DMDACreate(), DMDestroy(), DMDA
470 @*/
DMDAGetStencilType(DM da,DMDAStencilType * stype)471 PetscErrorCode DMDAGetStencilType(DM da, DMDAStencilType *stype)
472 {
473 DM_DA *dd = (DM_DA*)da->data;
474
475 PetscFunctionBegin;
476 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
477 PetscValidPointer(stype,2);
478 *stype = dd->stencil_type;
479 PetscFunctionReturn(0);
480 }
481
482 /*@
483 DMDASetStencilWidth - Sets the width of the communication stencil
484
485 Logically Collective on da
486
487 Input Parameter:
488 + da - The DMDA
489 - width - The stencil width
490
491 Level: intermediate
492
493 .seealso: DMDACreate(), DMDestroy(), DMDA
494 @*/
DMDASetStencilWidth(DM da,PetscInt width)495 PetscErrorCode DMDASetStencilWidth(DM da, PetscInt width)
496 {
497 DM_DA *dd = (DM_DA*)da->data;
498
499 PetscFunctionBegin;
500 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
501 PetscValidLogicalCollectiveInt(da,width,2);
502 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
503 dd->s = width;
504 PetscFunctionReturn(0);
505 }
506
507 /*@
508 DMDAGetStencilWidth - Gets the width of the communication stencil
509
510 Not collective
511
512 Input Parameter:
513 . da - The DMDA
514
515 Output Parameter:
516 . width - The stencil width
517
518 Level: intermediate
519
520 .seealso: DMDACreate(), DMDestroy(), DMDA
521 @*/
DMDAGetStencilWidth(DM da,PetscInt * width)522 PetscErrorCode DMDAGetStencilWidth(DM da, PetscInt *width)
523 {
524 DM_DA *dd = (DM_DA *) da->data;
525
526 PetscFunctionBegin;
527 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
528 PetscValidPointer(width,2);
529 *width = dd->s;
530 PetscFunctionReturn(0);
531 }
532
DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])533 static PetscErrorCode DMDACheckOwnershipRanges_Private(DM da,PetscInt M,PetscInt m,const PetscInt lx[])
534 {
535 PetscInt i,sum;
536
537 PetscFunctionBegin;
538 if (M < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Global dimension not set");
539 for (i=sum=0; i<m; i++) sum += lx[i];
540 if (sum != M) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Ownership ranges sum to %D but global dimension is %D",sum,M);
541 PetscFunctionReturn(0);
542 }
543
544 /*@
545 DMDASetOwnershipRanges - Sets the number of nodes in each direction on each process
546
547 Logically Collective on da
548
549 Input Parameter:
550 + da - The DMDA
551 . lx - array containing number of nodes in the X direction on each process, or NULL. If non-null, must be of length da->m
552 . ly - array containing number of nodes in the Y direction on each process, or NULL. If non-null, must be of length da->n
553 - lz - array containing number of nodes in the Z direction on each process, or NULL. If non-null, must be of length da->p.
554
555 Level: intermediate
556
557 Note: these numbers are NOT multiplied by the number of dof per node.
558
559 .seealso: DMDACreate(), DMDestroy(), DMDA
560 @*/
DMDASetOwnershipRanges(DM da,const PetscInt lx[],const PetscInt ly[],const PetscInt lz[])561 PetscErrorCode DMDASetOwnershipRanges(DM da, const PetscInt lx[], const PetscInt ly[], const PetscInt lz[])
562 {
563 PetscErrorCode ierr;
564 DM_DA *dd = (DM_DA*)da->data;
565
566 PetscFunctionBegin;
567 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
568 if (da->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
569 if (lx) {
570 if (dd->m < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
571 ierr = DMDACheckOwnershipRanges_Private(da,dd->M,dd->m,lx);CHKERRQ(ierr);
572 if (!dd->lx) {
573 ierr = PetscMalloc1(dd->m, &dd->lx);CHKERRQ(ierr);
574 }
575 ierr = PetscArraycpy(dd->lx, lx, dd->m);CHKERRQ(ierr);
576 }
577 if (ly) {
578 if (dd->n < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
579 ierr = DMDACheckOwnershipRanges_Private(da,dd->N,dd->n,ly);CHKERRQ(ierr);
580 if (!dd->ly) {
581 ierr = PetscMalloc1(dd->n, &dd->ly);CHKERRQ(ierr);
582 }
583 ierr = PetscArraycpy(dd->ly, ly, dd->n);CHKERRQ(ierr);
584 }
585 if (lz) {
586 if (dd->p < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of procs");
587 ierr = DMDACheckOwnershipRanges_Private(da,dd->P,dd->p,lz);CHKERRQ(ierr);
588 if (!dd->lz) {
589 ierr = PetscMalloc1(dd->p, &dd->lz);CHKERRQ(ierr);
590 }
591 ierr = PetscArraycpy(dd->lz, lz, dd->p);CHKERRQ(ierr);
592 }
593 PetscFunctionReturn(0);
594 }
595
596 /*@
597 DMDASetInterpolationType - Sets the type of interpolation that will be
598 returned by DMCreateInterpolation()
599
600 Logically Collective on da
601
602 Input Parameter:
603 + da - initial distributed array
604 - ctype - DMDA_Q1 and DMDA_Q0 are currently the only supported forms
605
606 Level: intermediate
607
608 Notes:
609 you should call this on the coarser of the two DMDAs you pass to DMCreateInterpolation()
610
611 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType
612 @*/
DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)613 PetscErrorCode DMDASetInterpolationType(DM da,DMDAInterpolationType ctype)
614 {
615 DM_DA *dd = (DM_DA*)da->data;
616
617 PetscFunctionBegin;
618 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
619 PetscValidLogicalCollectiveEnum(da,ctype,2);
620 dd->interptype = ctype;
621 PetscFunctionReturn(0);
622 }
623
624 /*@
625 DMDAGetInterpolationType - Gets the type of interpolation that will be
626 used by DMCreateInterpolation()
627
628 Not Collective
629
630 Input Parameter:
631 . da - distributed array
632
633 Output Parameter:
634 . ctype - interpolation type (DMDA_Q1 and DMDA_Q0 are currently the only supported forms)
635
636 Level: intermediate
637
638 .seealso: DMDA, DMDAInterpolationType, DMDASetInterpolationType(), DMCreateInterpolation()
639 @*/
DMDAGetInterpolationType(DM da,DMDAInterpolationType * ctype)640 PetscErrorCode DMDAGetInterpolationType(DM da,DMDAInterpolationType *ctype)
641 {
642 DM_DA *dd = (DM_DA*)da->data;
643
644 PetscFunctionBegin;
645 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
646 PetscValidPointer(ctype,2);
647 *ctype = dd->interptype;
648 PetscFunctionReturn(0);
649 }
650
651 /*@C
652 DMDAGetNeighbors - Gets an array containing the MPI rank of all the current
653 processes neighbors.
654
655 Not Collective
656
657 Input Parameter:
658 . da - the DMDA object
659
660 Output Parameters:
661 . ranks - the neighbors ranks, stored with the x index increasing most rapidly.
662 this process itself is in the list
663
664 Notes:
665 In 2d the array is of length 9, in 3d of length 27
666 Not supported in 1d
667 Do not free the array, it is freed when the DMDA is destroyed.
668
669 Fortran Notes:
670 In fortran you must pass in an array of the appropriate length.
671
672 Level: intermediate
673
674 @*/
DMDAGetNeighbors(DM da,const PetscMPIInt * ranks[])675 PetscErrorCode DMDAGetNeighbors(DM da,const PetscMPIInt *ranks[])
676 {
677 DM_DA *dd = (DM_DA*)da->data;
678
679 PetscFunctionBegin;
680 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
681 *ranks = dd->neighbors;
682 PetscFunctionReturn(0);
683 }
684
685 /*@C
686 DMDAGetOwnershipRanges - Gets the ranges of indices in the x, y and z direction that are owned by each process
687
688 Not Collective
689
690 Input Parameter:
691 . da - the DMDA object
692
693 Output Parameter:
694 + lx - ownership along x direction (optional)
695 . ly - ownership along y direction (optional)
696 - lz - ownership along z direction (optional)
697
698 Level: intermediate
699
700 Note: these correspond to the optional final arguments passed to DMDACreate(), DMDACreate2d(), DMDACreate3d()
701
702 In Fortran one must pass in arrays lx, ly, and lz that are long enough to hold the values; the sixth, seventh and
703 eighth arguments from DMDAGetInfo()
704
705 In C you should not free these arrays, nor change the values in them. They will only have valid values while the
706 DMDA they came from still exists (has not been destroyed).
707
708 These numbers are NOT multiplied by the number of dof per node.
709
710 .seealso: DMDAGetCorners(), DMDAGetGhostCorners(), DMDACreate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), VecGetOwnershipRanges()
711 @*/
DMDAGetOwnershipRanges(DM da,const PetscInt * lx[],const PetscInt * ly[],const PetscInt * lz[])712 PetscErrorCode DMDAGetOwnershipRanges(DM da,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
713 {
714 DM_DA *dd = (DM_DA*)da->data;
715
716 PetscFunctionBegin;
717 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
718 if (lx) *lx = dd->lx;
719 if (ly) *ly = dd->ly;
720 if (lz) *lz = dd->lz;
721 PetscFunctionReturn(0);
722 }
723
724 /*@
725 DMDASetRefinementFactor - Set the ratios that the DMDA grid is refined
726
727 Logically Collective on da
728
729 Input Parameters:
730 + da - the DMDA object
731 . refine_x - ratio of fine grid to coarse in x direction (2 by default)
732 . refine_y - ratio of fine grid to coarse in y direction (2 by default)
733 - refine_z - ratio of fine grid to coarse in z direction (2 by default)
734
735 Options Database:
736 + -da_refine_x - refinement ratio in x direction
737 . -da_refine_y - refinement ratio in y direction
738 - -da_refine_z - refinement ratio in z direction
739
740 Level: intermediate
741
742 Notes:
743 Pass PETSC_IGNORE to leave a value unchanged
744
745 .seealso: DMRefine(), DMDAGetRefinementFactor()
746 @*/
DMDASetRefinementFactor(DM da,PetscInt refine_x,PetscInt refine_y,PetscInt refine_z)747 PetscErrorCode DMDASetRefinementFactor(DM da, PetscInt refine_x, PetscInt refine_y,PetscInt refine_z)
748 {
749 DM_DA *dd = (DM_DA*)da->data;
750
751 PetscFunctionBegin;
752 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
753 PetscValidLogicalCollectiveInt(da,refine_x,2);
754 PetscValidLogicalCollectiveInt(da,refine_y,3);
755 PetscValidLogicalCollectiveInt(da,refine_z,4);
756
757 if (refine_x > 0) dd->refine_x = refine_x;
758 if (refine_y > 0) dd->refine_y = refine_y;
759 if (refine_z > 0) dd->refine_z = refine_z;
760 PetscFunctionReturn(0);
761 }
762
763 /*@C
764 DMDAGetRefinementFactor - Gets the ratios that the DMDA grid is refined
765
766 Not Collective
767
768 Input Parameter:
769 . da - the DMDA object
770
771 Output Parameters:
772 + refine_x - ratio of fine grid to coarse in x direction (2 by default)
773 . refine_y - ratio of fine grid to coarse in y direction (2 by default)
774 - refine_z - ratio of fine grid to coarse in z direction (2 by default)
775
776 Level: intermediate
777
778 Notes:
779 Pass NULL for values you do not need
780
781 .seealso: DMRefine(), DMDASetRefinementFactor()
782 @*/
DMDAGetRefinementFactor(DM da,PetscInt * refine_x,PetscInt * refine_y,PetscInt * refine_z)783 PetscErrorCode DMDAGetRefinementFactor(DM da, PetscInt *refine_x, PetscInt *refine_y,PetscInt *refine_z)
784 {
785 DM_DA *dd = (DM_DA*)da->data;
786
787 PetscFunctionBegin;
788 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
789 if (refine_x) *refine_x = dd->refine_x;
790 if (refine_y) *refine_y = dd->refine_y;
791 if (refine_z) *refine_z = dd->refine_z;
792 PetscFunctionReturn(0);
793 }
794
795 /*@C
796 DMDASetGetMatrix - Sets the routine used by the DMDA to allocate a matrix.
797
798 Logically Collective on da
799
800 Input Parameters:
801 + da - the DMDA object
802 - f - the function that allocates the matrix for that specific DMDA
803
804 Level: developer
805
806 Notes:
807 See DMDASetBlockFills() that provides a simple way to provide the nonzero structure for
808 the diagonal and off-diagonal blocks of the matrix
809
810 Not supported from Fortran
811
812 .seealso: DMCreateMatrix(), DMDASetBlockFills()
813 @*/
DMDASetGetMatrix(DM da,PetscErrorCode (* f)(DM,Mat *))814 PetscErrorCode DMDASetGetMatrix(DM da,PetscErrorCode (*f)(DM, Mat*))
815 {
816 PetscFunctionBegin;
817 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
818 da->ops->creatematrix = f;
819 PetscFunctionReturn(0);
820 }
821
822 /*
823 Creates "balanced" ownership ranges after refinement, constrained by the need for the
824 fine grid boundaries to fall within one stencil width of the coarse partition.
825
826 Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
827 */
DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])828 static PetscErrorCode DMDARefineOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lc[],PetscInt lf[])
829 {
830 PetscInt i,totalc = 0,remaining,startc = 0,startf = 0;
831 PetscErrorCode ierr;
832
833 PetscFunctionBegin;
834 if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
835 if (ratio == 1) {
836 ierr = PetscArraycpy(lf,lc,m);CHKERRQ(ierr);
837 PetscFunctionReturn(0);
838 }
839 for (i=0; i<m; i++) totalc += lc[i];
840 remaining = (!periodic) + ratio * (totalc - (!periodic));
841 for (i=0; i<m; i++) {
842 PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
843 if (i == m-1) lf[i] = want;
844 else {
845 const PetscInt nextc = startc + lc[i];
846 /* Move the first fine node of the next subdomain to the right until the coarse node on its left is within one
847 * coarse stencil width of the first coarse node in the next subdomain. */
848 while ((startf+want)/ratio < nextc - stencil_width) want++;
849 /* Move the last fine node in the current subdomain to the left until the coarse node on its right is within one
850 * coarse stencil width of the last coarse node in the current subdomain. */
851 while ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width) want--;
852 /* Make sure all constraints are satisfied */
853 if (want < 0 || want > remaining || ((startf+want)/ratio < nextc - stencil_width)
854 || ((startf+want-1+ratio-1)/ratio > nextc-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible refined ownership range");
855 }
856 lf[i] = want;
857 startc += lc[i];
858 startf += lf[i];
859 remaining -= lf[i];
860 }
861 PetscFunctionReturn(0);
862 }
863
864 /*
865 Creates "balanced" ownership ranges after coarsening, constrained by the need for the
866 fine grid boundaries to fall within one stencil width of the coarse partition.
867
868 Uses a greedy algorithm to handle non-ideal layouts, could probably do something better.
869 */
DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])870 static PetscErrorCode DMDACoarsenOwnershipRanges(DM da,PetscBool periodic,PetscInt stencil_width,PetscInt ratio,PetscInt m,const PetscInt lf[],PetscInt lc[])
871 {
872 PetscInt i,totalf,remaining,startc,startf;
873 PetscErrorCode ierr;
874
875 PetscFunctionBegin;
876 if (ratio < 1) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Requested refinement ratio %D must be at least 1",ratio);
877 if (ratio == 1) {
878 ierr = PetscArraycpy(lc,lf,m);CHKERRQ(ierr);
879 PetscFunctionReturn(0);
880 }
881 for (i=0,totalf=0; i<m; i++) totalf += lf[i];
882 remaining = (!periodic) + (totalf - (!periodic)) / ratio;
883 for (i=0,startc=0,startf=0; i<m; i++) {
884 PetscInt want = remaining/(m-i) + !!(remaining%(m-i));
885 if (i == m-1) lc[i] = want;
886 else {
887 const PetscInt nextf = startf+lf[i];
888 /* Slide first coarse node of next subdomain to the left until the coarse node to the left of the first fine
889 * node is within one stencil width. */
890 while (nextf/ratio < startc+want-stencil_width) want--;
891 /* Slide the last coarse node of the current subdomain to the right until the coarse node to the right of the last
892 * fine node is within one stencil width. */
893 while ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width) want++;
894 if (want < 0 || want > remaining
895 || (nextf/ratio < startc+want-stencil_width) || ((nextf-1+ratio-1)/ratio > startc+want-1+stencil_width)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Could not find a compatible coarsened ownership range");
896 }
897 lc[i] = want;
898 startc += lc[i];
899 startf += lf[i];
900 remaining -= lc[i];
901 }
902 PetscFunctionReturn(0);
903 }
904
DMRefine_DA(DM da,MPI_Comm comm,DM * daref)905 PetscErrorCode DMRefine_DA(DM da,MPI_Comm comm,DM *daref)
906 {
907 PetscErrorCode ierr;
908 PetscInt M,N,P,i,dim;
909 DM da2;
910 DM_DA *dd = (DM_DA*)da->data,*dd2;
911
912 PetscFunctionBegin;
913 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA);
914 PetscValidPointer(daref,3);
915
916 ierr = DMGetDimension(da, &dim);CHKERRQ(ierr);
917 if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
918 M = dd->refine_x*dd->M;
919 } else {
920 M = 1 + dd->refine_x*(dd->M - 1);
921 }
922 if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
923 if (dim > 1) {
924 N = dd->refine_y*dd->N;
925 } else {
926 N = 1;
927 }
928 } else {
929 N = 1 + dd->refine_y*(dd->N - 1);
930 }
931 if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
932 if (dim > 2) {
933 P = dd->refine_z*dd->P;
934 } else {
935 P = 1;
936 }
937 } else {
938 P = 1 + dd->refine_z*(dd->P - 1);
939 }
940 ierr = DMDACreate(PetscObjectComm((PetscObject)da),&da2);CHKERRQ(ierr);
941 ierr = DMSetOptionsPrefix(da2,((PetscObject)da)->prefix);CHKERRQ(ierr);
942 ierr = DMSetDimension(da2,dim);CHKERRQ(ierr);
943 ierr = DMDASetSizes(da2,M,N,P);CHKERRQ(ierr);
944 ierr = DMDASetNumProcs(da2,dd->m,dd->n,dd->p);CHKERRQ(ierr);
945 ierr = DMDASetBoundaryType(da2,dd->bx,dd->by,dd->bz);CHKERRQ(ierr);
946 ierr = DMDASetDof(da2,dd->w);CHKERRQ(ierr);
947 ierr = DMDASetStencilType(da2,dd->stencil_type);CHKERRQ(ierr);
948 ierr = DMDASetStencilWidth(da2,dd->s);CHKERRQ(ierr);
949 if (dim == 3) {
950 PetscInt *lx,*ly,*lz;
951 ierr = PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);CHKERRQ(ierr);
952 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
953 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
954 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_z,dd->p,dd->lz,lz);CHKERRQ(ierr);
955 ierr = DMDASetOwnershipRanges(da2,lx,ly,lz);CHKERRQ(ierr);
956 ierr = PetscFree3(lx,ly,lz);CHKERRQ(ierr);
957 } else if (dim == 2) {
958 PetscInt *lx,*ly;
959 ierr = PetscMalloc2(dd->m,&lx,dd->n,&ly);CHKERRQ(ierr);
960 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
961 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
962 ierr = DMDASetOwnershipRanges(da2,lx,ly,NULL);CHKERRQ(ierr);
963 ierr = PetscFree2(lx,ly);CHKERRQ(ierr);
964 } else if (dim == 1) {
965 PetscInt *lx;
966 ierr = PetscMalloc1(dd->m,&lx);CHKERRQ(ierr);
967 ierr = DMDARefineOwnershipRanges(da,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->refine_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
968 ierr = DMDASetOwnershipRanges(da2,lx,NULL,NULL);CHKERRQ(ierr);
969 ierr = PetscFree(lx);CHKERRQ(ierr);
970 }
971 dd2 = (DM_DA*)da2->data;
972
973 /* allow overloaded (user replaced) operations to be inherited by refinement clones */
974 da2->ops->creatematrix = da->ops->creatematrix;
975 /* da2->ops->createinterpolation = da->ops->createinterpolation; this causes problem with SNESVI */
976 da2->ops->getcoloring = da->ops->getcoloring;
977 dd2->interptype = dd->interptype;
978
979 /* copy fill information if given */
980 if (dd->dfill) {
981 ierr = PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);CHKERRQ(ierr);
982 ierr = PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);CHKERRQ(ierr);
983 }
984 if (dd->ofill) {
985 ierr = PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);CHKERRQ(ierr);
986 ierr = PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);CHKERRQ(ierr);
987 }
988 /* copy the refine information */
989 dd2->coarsen_x = dd2->refine_x = dd->refine_x;
990 dd2->coarsen_y = dd2->refine_y = dd->refine_y;
991 dd2->coarsen_z = dd2->refine_z = dd->refine_z;
992
993 if (dd->refine_z_hier) {
994 if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_z_hier_n) {
995 dd2->refine_z = dd->refine_z_hier[da->levelup - da->leveldown + 1];
996 }
997 if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_z_hier_n) {
998 dd2->coarsen_z = dd->refine_z_hier[da->levelup - da->leveldown];
999 }
1000 dd2->refine_z_hier_n = dd->refine_z_hier_n;
1001 ierr = PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);CHKERRQ(ierr);
1002 ierr = PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);CHKERRQ(ierr);
1003 }
1004 if (dd->refine_y_hier) {
1005 if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_y_hier_n) {
1006 dd2->refine_y = dd->refine_y_hier[da->levelup - da->leveldown + 1];
1007 }
1008 if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_y_hier_n) {
1009 dd2->coarsen_y = dd->refine_y_hier[da->levelup - da->leveldown];
1010 }
1011 dd2->refine_y_hier_n = dd->refine_y_hier_n;
1012 ierr = PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);CHKERRQ(ierr);
1013 ierr = PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);CHKERRQ(ierr);
1014 }
1015 if (dd->refine_x_hier) {
1016 if (da->levelup - da->leveldown + 1 > -1 && da->levelup - da->leveldown + 1 < dd->refine_x_hier_n) {
1017 dd2->refine_x = dd->refine_x_hier[da->levelup - da->leveldown + 1];
1018 }
1019 if (da->levelup - da->leveldown > -1 && da->levelup - da->leveldown < dd->refine_x_hier_n) {
1020 dd2->coarsen_x = dd->refine_x_hier[da->levelup - da->leveldown];
1021 }
1022 dd2->refine_x_hier_n = dd->refine_x_hier_n;
1023 ierr = PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);CHKERRQ(ierr);
1024 ierr = PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);CHKERRQ(ierr);
1025 }
1026
1027
1028 /* copy vector type information */
1029 ierr = DMSetVecType(da2,da->vectype);CHKERRQ(ierr);
1030
1031 dd2->lf = dd->lf;
1032 dd2->lj = dd->lj;
1033
1034 da2->leveldown = da->leveldown;
1035 da2->levelup = da->levelup + 1;
1036
1037 ierr = DMSetUp(da2);CHKERRQ(ierr);
1038
1039 /* interpolate coordinates if they are set on the coarse grid */
1040 if (da->coordinates) {
1041 DM cdaf,cdac;
1042 Vec coordsc,coordsf;
1043 Mat II;
1044
1045 ierr = DMGetCoordinateDM(da,&cdac);CHKERRQ(ierr);
1046 ierr = DMGetCoordinates(da,&coordsc);CHKERRQ(ierr);
1047 ierr = DMGetCoordinateDM(da2,&cdaf);CHKERRQ(ierr);
1048 /* force creation of the coordinate vector */
1049 ierr = DMDASetUniformCoordinates(da2,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
1050 ierr = DMGetCoordinates(da2,&coordsf);CHKERRQ(ierr);
1051 ierr = DMCreateInterpolation(cdac,cdaf,&II,NULL);CHKERRQ(ierr);
1052 ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
1053 ierr = MatDestroy(&II);CHKERRQ(ierr);
1054 }
1055
1056 for (i=0; i<da->bs; i++) {
1057 const char *fieldname;
1058 ierr = DMDAGetFieldName(da,i,&fieldname);CHKERRQ(ierr);
1059 ierr = DMDASetFieldName(da2,i,fieldname);CHKERRQ(ierr);
1060 }
1061
1062 *daref = da2;
1063 PetscFunctionReturn(0);
1064 }
1065
1066
DMCoarsen_DA(DM dmf,MPI_Comm comm,DM * dmc)1067 PetscErrorCode DMCoarsen_DA(DM dmf, MPI_Comm comm,DM *dmc)
1068 {
1069 PetscErrorCode ierr;
1070 PetscInt M,N,P,i,dim;
1071 DM dmc2;
1072 DM_DA *dd = (DM_DA*)dmf->data,*dd2;
1073
1074 PetscFunctionBegin;
1075 PetscValidHeaderSpecificType(dmf,DM_CLASSID,1,DMDA);
1076 PetscValidPointer(dmc,3);
1077
1078 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr);
1079 if (dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1080 M = dd->M / dd->coarsen_x;
1081 } else {
1082 M = 1 + (dd->M - 1) / dd->coarsen_x;
1083 }
1084 if (dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1085 if (dim > 1) {
1086 N = dd->N / dd->coarsen_y;
1087 } else {
1088 N = 1;
1089 }
1090 } else {
1091 N = 1 + (dd->N - 1) / dd->coarsen_y;
1092 }
1093 if (dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) {
1094 if (dim > 2) {
1095 P = dd->P / dd->coarsen_z;
1096 } else {
1097 P = 1;
1098 }
1099 } else {
1100 P = 1 + (dd->P - 1) / dd->coarsen_z;
1101 }
1102 ierr = DMDACreate(PetscObjectComm((PetscObject)dmf),&dmc2);CHKERRQ(ierr);
1103 ierr = DMSetOptionsPrefix(dmc2,((PetscObject)dmf)->prefix);CHKERRQ(ierr);
1104 ierr = DMSetDimension(dmc2,dim);CHKERRQ(ierr);
1105 ierr = DMDASetSizes(dmc2,M,N,P);CHKERRQ(ierr);
1106 ierr = DMDASetNumProcs(dmc2,dd->m,dd->n,dd->p);CHKERRQ(ierr);
1107 ierr = DMDASetBoundaryType(dmc2,dd->bx,dd->by,dd->bz);CHKERRQ(ierr);
1108 ierr = DMDASetDof(dmc2,dd->w);CHKERRQ(ierr);
1109 ierr = DMDASetStencilType(dmc2,dd->stencil_type);CHKERRQ(ierr);
1110 ierr = DMDASetStencilWidth(dmc2,dd->s);CHKERRQ(ierr);
1111 if (dim == 3) {
1112 PetscInt *lx,*ly,*lz;
1113 ierr = PetscMalloc3(dd->m,&lx,dd->n,&ly,dd->p,&lz);CHKERRQ(ierr);
1114 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
1115 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
1116 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bz == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_z,dd->p,dd->lz,lz);CHKERRQ(ierr);
1117 ierr = DMDASetOwnershipRanges(dmc2,lx,ly,lz);CHKERRQ(ierr);
1118 ierr = PetscFree3(lx,ly,lz);CHKERRQ(ierr);
1119 } else if (dim == 2) {
1120 PetscInt *lx,*ly;
1121 ierr = PetscMalloc2(dd->m,&lx,dd->n,&ly);CHKERRQ(ierr);
1122 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
1123 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->by == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_y,dd->n,dd->ly,ly);CHKERRQ(ierr);
1124 ierr = DMDASetOwnershipRanges(dmc2,lx,ly,NULL);CHKERRQ(ierr);
1125 ierr = PetscFree2(lx,ly);CHKERRQ(ierr);
1126 } else if (dim == 1) {
1127 PetscInt *lx;
1128 ierr = PetscMalloc1(dd->m,&lx);CHKERRQ(ierr);
1129 ierr = DMDACoarsenOwnershipRanges(dmf,(PetscBool)(dd->bx == DM_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0),dd->s,dd->coarsen_x,dd->m,dd->lx,lx);CHKERRQ(ierr);
1130 ierr = DMDASetOwnershipRanges(dmc2,lx,NULL,NULL);CHKERRQ(ierr);
1131 ierr = PetscFree(lx);CHKERRQ(ierr);
1132 }
1133 dd2 = (DM_DA*)dmc2->data;
1134
1135 /* allow overloaded (user replaced) operations to be inherited by refinement clones; why are only some inherited and not all? */
1136 /* dmc2->ops->createinterpolation = dmf->ops->createinterpolation; copying this one causes trouble for DMSetVI */
1137 dmc2->ops->creatematrix = dmf->ops->creatematrix;
1138 dmc2->ops->getcoloring = dmf->ops->getcoloring;
1139 dd2->interptype = dd->interptype;
1140
1141 /* copy fill information if given */
1142 if (dd->dfill) {
1143 ierr = PetscMalloc1(dd->dfill[dd->w]+dd->w+1,&dd2->dfill);CHKERRQ(ierr);
1144 ierr = PetscArraycpy(dd2->dfill,dd->dfill,dd->dfill[dd->w]+dd->w+1);CHKERRQ(ierr);
1145 }
1146 if (dd->ofill) {
1147 ierr = PetscMalloc1(dd->ofill[dd->w]+dd->w+1,&dd2->ofill);CHKERRQ(ierr);
1148 ierr = PetscArraycpy(dd2->ofill,dd->ofill,dd->ofill[dd->w]+dd->w+1);CHKERRQ(ierr);
1149 }
1150 /* copy the refine information */
1151 dd2->coarsen_x = dd2->refine_x = dd->coarsen_x;
1152 dd2->coarsen_y = dd2->refine_y = dd->coarsen_y;
1153 dd2->coarsen_z = dd2->refine_z = dd->coarsen_z;
1154
1155 if (dd->refine_z_hier) {
1156 if (dmf->levelup - dmf->leveldown -1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_z_hier_n) {
1157 dd2->refine_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 1];
1158 }
1159 if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_z_hier_n) {
1160 dd2->coarsen_z = dd->refine_z_hier[dmf->levelup - dmf->leveldown - 2];
1161 }
1162 dd2->refine_z_hier_n = dd->refine_z_hier_n;
1163 ierr = PetscMalloc1(dd2->refine_z_hier_n,&dd2->refine_z_hier);CHKERRQ(ierr);
1164 ierr = PetscArraycpy(dd2->refine_z_hier,dd->refine_z_hier,dd2->refine_z_hier_n);CHKERRQ(ierr);
1165 }
1166 if (dd->refine_y_hier) {
1167 if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1< dd->refine_y_hier_n) {
1168 dd2->refine_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 1];
1169 }
1170 if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_y_hier_n) {
1171 dd2->coarsen_y = dd->refine_y_hier[dmf->levelup - dmf->leveldown - 2];
1172 }
1173 dd2->refine_y_hier_n = dd->refine_y_hier_n;
1174 ierr = PetscMalloc1(dd2->refine_y_hier_n,&dd2->refine_y_hier);CHKERRQ(ierr);
1175 ierr = PetscArraycpy(dd2->refine_y_hier,dd->refine_y_hier,dd2->refine_y_hier_n);CHKERRQ(ierr);
1176 }
1177 if (dd->refine_x_hier) {
1178 if (dmf->levelup - dmf->leveldown - 1 > -1 && dmf->levelup - dmf->leveldown - 1 < dd->refine_x_hier_n) {
1179 dd2->refine_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 1];
1180 }
1181 if (dmf->levelup - dmf->leveldown - 2 > -1 && dmf->levelup - dmf->leveldown - 2 < dd->refine_x_hier_n) {
1182 dd2->coarsen_x = dd->refine_x_hier[dmf->levelup - dmf->leveldown - 2];
1183 }
1184 dd2->refine_x_hier_n = dd->refine_x_hier_n;
1185 ierr = PetscMalloc1(dd2->refine_x_hier_n,&dd2->refine_x_hier);CHKERRQ(ierr);
1186 ierr = PetscArraycpy(dd2->refine_x_hier,dd->refine_x_hier,dd2->refine_x_hier_n);CHKERRQ(ierr);
1187 }
1188
1189 /* copy vector type information */
1190 ierr = DMSetVecType(dmc2,dmf->vectype);CHKERRQ(ierr);
1191
1192 dd2->lf = dd->lf;
1193 dd2->lj = dd->lj;
1194
1195 dmc2->leveldown = dmf->leveldown + 1;
1196 dmc2->levelup = dmf->levelup;
1197
1198 ierr = DMSetUp(dmc2);CHKERRQ(ierr);
1199
1200 /* inject coordinates if they are set on the fine grid */
1201 if (dmf->coordinates) {
1202 DM cdaf,cdac;
1203 Vec coordsc,coordsf;
1204 Mat inject;
1205 VecScatter vscat;
1206
1207 ierr = DMGetCoordinateDM(dmf,&cdaf);CHKERRQ(ierr);
1208 ierr = DMGetCoordinates(dmf,&coordsf);CHKERRQ(ierr);
1209 ierr = DMGetCoordinateDM(dmc2,&cdac);CHKERRQ(ierr);
1210 /* force creation of the coordinate vector */
1211 ierr = DMDASetUniformCoordinates(dmc2,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);
1212 ierr = DMGetCoordinates(dmc2,&coordsc);CHKERRQ(ierr);
1213
1214 ierr = DMCreateInjection(cdac,cdaf,&inject);CHKERRQ(ierr);
1215 ierr = MatScatterGetVecScatter(inject,&vscat);CHKERRQ(ierr);
1216 ierr = VecScatterBegin(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1217 ierr = VecScatterEnd(vscat,coordsf,coordsc,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
1218 ierr = MatDestroy(&inject);CHKERRQ(ierr);
1219 }
1220
1221 for (i=0; i<dmf->bs; i++) {
1222 const char *fieldname;
1223 ierr = DMDAGetFieldName(dmf,i,&fieldname);CHKERRQ(ierr);
1224 ierr = DMDASetFieldName(dmc2,i,fieldname);CHKERRQ(ierr);
1225 }
1226
1227 *dmc = dmc2;
1228 PetscFunctionReturn(0);
1229 }
1230
DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])1231 PetscErrorCode DMRefineHierarchy_DA(DM da,PetscInt nlevels,DM daf[])
1232 {
1233 PetscErrorCode ierr;
1234 PetscInt i,n,*refx,*refy,*refz;
1235
1236 PetscFunctionBegin;
1237 PetscValidHeaderSpecific(da,DM_CLASSID,1);
1238 if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1239 if (nlevels == 0) PetscFunctionReturn(0);
1240 PetscValidPointer(daf,3);
1241
1242 /* Get refinement factors, defaults taken from the coarse DMDA */
1243 ierr = PetscMalloc3(nlevels,&refx,nlevels,&refy,nlevels,&refz);CHKERRQ(ierr);
1244 for (i=0; i<nlevels; i++) {
1245 ierr = DMDAGetRefinementFactor(da,&refx[i],&refy[i],&refz[i]);CHKERRQ(ierr);
1246 }
1247 n = nlevels;
1248 ierr = PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);CHKERRQ(ierr);
1249 n = nlevels;
1250 ierr = PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);CHKERRQ(ierr);
1251 n = nlevels;
1252 ierr = PetscOptionsGetIntArray(((PetscObject)da)->options,((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);CHKERRQ(ierr);
1253
1254 ierr = DMDASetRefinementFactor(da,refx[0],refy[0],refz[0]);CHKERRQ(ierr);
1255 ierr = DMRefine(da,PetscObjectComm((PetscObject)da),&daf[0]);CHKERRQ(ierr);
1256 for (i=1; i<nlevels; i++) {
1257 ierr = DMDASetRefinementFactor(daf[i-1],refx[i],refy[i],refz[i]);CHKERRQ(ierr);
1258 ierr = DMRefine(daf[i-1],PetscObjectComm((PetscObject)da),&daf[i]);CHKERRQ(ierr);
1259 }
1260 ierr = PetscFree3(refx,refy,refz);CHKERRQ(ierr);
1261 PetscFunctionReturn(0);
1262 }
1263
DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])1264 PetscErrorCode DMCoarsenHierarchy_DA(DM da,PetscInt nlevels,DM dac[])
1265 {
1266 PetscErrorCode ierr;
1267 PetscInt i;
1268
1269 PetscFunctionBegin;
1270 PetscValidHeaderSpecific(da,DM_CLASSID,1);
1271 if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
1272 if (nlevels == 0) PetscFunctionReturn(0);
1273 PetscValidPointer(dac,3);
1274 ierr = DMCoarsen(da,PetscObjectComm((PetscObject)da),&dac[0]);CHKERRQ(ierr);
1275 for (i=1; i<nlevels; i++) {
1276 ierr = DMCoarsen(dac[i-1],PetscObjectComm((PetscObject)da),&dac[i]);CHKERRQ(ierr);
1277 }
1278 PetscFunctionReturn(0);
1279 }
1280
DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal * nodes)1281 PetscErrorCode DMDASetGLLCoordinates_1d(DM dm,PetscInt n,PetscReal *nodes)
1282 {
1283 PetscErrorCode ierr;
1284 PetscInt i,j,xs,xn,q;
1285 PetscScalar *xx;
1286 PetscReal h;
1287 Vec x;
1288 DM_DA *da = (DM_DA*)dm->data;
1289
1290 PetscFunctionBegin;
1291 if (da->bx != DM_BOUNDARY_PERIODIC) {
1292 ierr = DMDAGetInfo(dm,NULL,&q,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1293 q = (q-1)/(n-1); /* number of spectral elements */
1294 h = 2.0/q;
1295 ierr = DMDAGetCorners(dm,&xs,NULL,NULL,&xn,NULL,NULL);CHKERRQ(ierr);
1296 xs = xs/(n-1);
1297 xn = xn/(n-1);
1298 ierr = DMDASetUniformCoordinates(dm,-1.,1.,0.,0.,0.,0.);CHKERRQ(ierr);
1299 ierr = DMGetCoordinates(dm,&x);CHKERRQ(ierr);
1300 ierr = DMDAVecGetArray(dm,x,&xx);CHKERRQ(ierr);
1301
1302 /* loop over local spectral elements */
1303 for (j=xs; j<xs+xn; j++) {
1304 /*
1305 Except for the first process, each process starts on the second GLL point of the first element on that process
1306 */
1307 for (i= (j == xs && xs > 0)? 1 : 0; i<n; i++) {
1308 xx[j*(n-1) + i] = -1.0 + h*j + h*(nodes[i]+1.0)/2.;
1309 }
1310 }
1311 ierr = DMDAVecRestoreArray(dm,x,&xx);CHKERRQ(ierr);
1312 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for periodic");
1313 PetscFunctionReturn(0);
1314 }
1315
1316 /*@
1317
1318 DMDASetGLLCoordinates - Sets the global coordinates from -1 to 1 to the GLL points of as many GLL elements that fit the number of grid points
1319
1320 Collective on da
1321
1322 Input Parameters:
1323 + da - the DMDA object
1324 - n - the number of GLL nodes
1325 - nodes - the GLL nodes
1326
1327 Notes:
1328 the parallel decomposition of grid points must correspond to the degree of the GLL. That is, the number of grid points
1329 on each process much be divisible by the number of GLL elements needed per process. This depends on whether the DM is
1330 periodic or not.
1331
1332 Level: advanced
1333
1334 .seealso: DMDACreate(), PetscDTGaussLobattoLegendreQuadrature(), DMGetCoordinates()
1335 @*/
DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal * nodes)1336 PetscErrorCode DMDASetGLLCoordinates(DM da,PetscInt n,PetscReal *nodes)
1337 {
1338 PetscErrorCode ierr;
1339
1340 PetscFunctionBegin;
1341 if (da->dim == 1) {
1342 ierr = DMDASetGLLCoordinates_1d(da,n,nodes);CHKERRQ(ierr);
1343 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not yet implemented for 2 or 3d");
1344 PetscFunctionReturn(0);
1345 }
1346
DMGetCompatibility_DA(DM da1,DM dm2,PetscBool * compatible,PetscBool * set)1347 PETSC_INTERN PetscErrorCode DMGetCompatibility_DA(DM da1,DM dm2,PetscBool *compatible,PetscBool *set)
1348 {
1349 PetscErrorCode ierr;
1350 DM_DA *dd1 = (DM_DA*)da1->data,*dd2;
1351 DM da2;
1352 DMType dmtype2;
1353 PetscBool isda,compatibleLocal;
1354 PetscInt i;
1355
1356 PetscFunctionBegin;
1357 if (!da1->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da1),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on first DM before DMGetCompatibility()");
1358 ierr = DMGetType(dm2,&dmtype2);CHKERRQ(ierr);
1359 ierr = PetscStrcmp(dmtype2,DMDA,&isda);CHKERRQ(ierr);
1360 if (isda) {
1361 da2 = dm2;
1362 dd2 = (DM_DA*)da2->data;
1363 if (!da2->setupcalled) SETERRQ(PetscObjectComm((PetscObject)da2),PETSC_ERR_ARG_WRONGSTATE,"DMSetUp() must be called on second DM before DMGetCompatibility()");
1364 compatibleLocal = (PetscBool)(da1->dim == da2->dim);
1365 if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->s == dd2->s)); /* Stencil width */
1366 /* Global size ranks Boundary type */
1367 if (compatibleLocal) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->M == dd2->M) && (dd1->m == dd2->m) && (dd1->bx == dd2->bx));
1368 if (compatibleLocal && da1->dim > 1) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->N == dd2->N) && (dd1->n == dd2->n) && (dd1->by == dd2->by));
1369 if (compatibleLocal && da1->dim > 2) compatibleLocal = (PetscBool)(compatibleLocal && (dd1->P == dd2->P) && (dd1->p == dd2->p) && (dd1->bz == dd2->bz));
1370 if (compatibleLocal) {
1371 for (i=0; i<dd1->m; ++i) {
1372 compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lx[i] == dd2->lx[i])); /* Local size */
1373 }
1374 }
1375 if (compatibleLocal && da1->dim > 1) {
1376 for (i=0; i<dd1->n; ++i) {
1377 compatibleLocal = (PetscBool)(compatibleLocal && (dd1->ly[i] == dd2->ly[i]));
1378 }
1379 }
1380 if (compatibleLocal && da1->dim > 2) {
1381 for (i=0; i<dd1->p; ++i) {
1382 compatibleLocal = (PetscBool)(compatibleLocal && (dd1->lz[i] == dd2->lz[i]));
1383 }
1384 }
1385 *compatible = compatibleLocal;
1386 *set = PETSC_TRUE;
1387 } else {
1388 /* Decline to determine compatibility with other DM types */
1389 *set = PETSC_FALSE;
1390 }
1391 PetscFunctionReturn(0);
1392 }
1393