1
2 #include <../src/dm/impls/composite/packimpl.h> /*I "petscdmcomposite.h" I*/
3 #include <petsc/private/isimpl.h>
4 #include <petsc/private/glvisviewerimpl.h>
5 #include <petscds.h>
6
7 /*@C
8 DMCompositeSetCoupling - Sets user provided routines that compute the coupling between the
9 separate components (DMs) in a DMto build the correct matrix nonzero structure.
10
11
12 Logically Collective
13
14 Input Parameter:
15 + dm - the composite object
16 - formcouplelocations - routine to set the nonzero locations in the matrix
17
18 Level: advanced
19
20 Not available from Fortran
21
22 Notes:
23 See DMSetApplicationContext() and DMGetApplicationContext() for how to get user information into
24 this routine
25
26 @*/
DMCompositeSetCoupling(DM dm,PetscErrorCode (* FormCoupleLocations)(DM,Mat,PetscInt *,PetscInt *,PetscInt,PetscInt,PetscInt,PetscInt))27 PetscErrorCode DMCompositeSetCoupling(DM dm,PetscErrorCode (*FormCoupleLocations)(DM,Mat,PetscInt*,PetscInt*,PetscInt,PetscInt,PetscInt,PetscInt))
28 {
29 DM_Composite *com = (DM_Composite*)dm->data;
30 PetscBool flg;
31 PetscErrorCode ierr;
32
33 PetscFunctionBegin;
34 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
35 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
36 com->FormCoupleLocations = FormCoupleLocations;
37 PetscFunctionReturn(0);
38 }
39
DMDestroy_Composite(DM dm)40 PetscErrorCode DMDestroy_Composite(DM dm)
41 {
42 PetscErrorCode ierr;
43 struct DMCompositeLink *next, *prev;
44 DM_Composite *com = (DM_Composite*)dm->data;
45
46 PetscFunctionBegin;
47 next = com->next;
48 while (next) {
49 prev = next;
50 next = next->next;
51 ierr = DMDestroy(&prev->dm);CHKERRQ(ierr);
52 ierr = PetscFree(prev->grstarts);CHKERRQ(ierr);
53 ierr = PetscFree(prev);CHKERRQ(ierr);
54 }
55 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",NULL);CHKERRQ(ierr);
56 /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
57 ierr = PetscFree(com);CHKERRQ(ierr);
58 PetscFunctionReturn(0);
59 }
60
DMView_Composite(DM dm,PetscViewer v)61 PetscErrorCode DMView_Composite(DM dm,PetscViewer v)
62 {
63 PetscErrorCode ierr;
64 PetscBool iascii;
65 DM_Composite *com = (DM_Composite*)dm->data;
66
67 PetscFunctionBegin;
68 ierr = PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
69 if (iascii) {
70 struct DMCompositeLink *lnk = com->next;
71 PetscInt i;
72
73 ierr = PetscViewerASCIIPrintf(v,"DM (%s)\n",((PetscObject)dm)->prefix ? ((PetscObject)dm)->prefix : "no prefix");CHKERRQ(ierr);
74 ierr = PetscViewerASCIIPrintf(v," contains %D DMs\n",com->nDM);CHKERRQ(ierr);
75 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr);
76 for (i=0; lnk; lnk=lnk->next,i++) {
77 ierr = PetscViewerASCIIPrintf(v,"Link %D: DM of type %s\n",i,((PetscObject)lnk->dm)->type_name);CHKERRQ(ierr);
78 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr);
79 ierr = DMView(lnk->dm,v);CHKERRQ(ierr);
80 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr);
81 }
82 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr);
83 }
84 PetscFunctionReturn(0);
85 }
86
87 /* --------------------------------------------------------------------------------------*/
DMSetUp_Composite(DM dm)88 PetscErrorCode DMSetUp_Composite(DM dm)
89 {
90 PetscErrorCode ierr;
91 PetscInt nprev = 0;
92 PetscMPIInt rank,size;
93 DM_Composite *com = (DM_Composite*)dm->data;
94 struct DMCompositeLink *next = com->next;
95 PetscLayout map;
96
97 PetscFunctionBegin;
98 if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Packer has already been setup");
99 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject)dm),&map);CHKERRQ(ierr);
100 ierr = PetscLayoutSetLocalSize(map,com->n);CHKERRQ(ierr);
101 ierr = PetscLayoutSetSize(map,PETSC_DETERMINE);CHKERRQ(ierr);
102 ierr = PetscLayoutSetBlockSize(map,1);CHKERRQ(ierr);
103 ierr = PetscLayoutSetUp(map);CHKERRQ(ierr);
104 ierr = PetscLayoutGetSize(map,&com->N);CHKERRQ(ierr);
105 ierr = PetscLayoutGetRange(map,&com->rstart,NULL);CHKERRQ(ierr);
106 ierr = PetscLayoutDestroy(&map);CHKERRQ(ierr);
107
108 /* now set the rstart for each linked vector */
109 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
110 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);CHKERRQ(ierr);
111 while (next) {
112 next->rstart = nprev;
113 nprev += next->n;
114 next->grstart = com->rstart + next->rstart;
115 ierr = PetscMalloc1(size,&next->grstarts);CHKERRQ(ierr);
116 ierr = MPI_Allgather(&next->grstart,1,MPIU_INT,next->grstarts,1,MPIU_INT,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
117 next = next->next;
118 }
119 com->setup = PETSC_TRUE;
120 PetscFunctionReturn(0);
121 }
122
123 /* ----------------------------------------------------------------------------------*/
124
125 /*@
126 DMCompositeGetNumberDM - Get's the number of DM objects in the DMComposite
127 representation.
128
129 Not Collective
130
131 Input Parameter:
132 . dm - the packer object
133
134 Output Parameter:
135 . nDM - the number of DMs
136
137 Level: beginner
138
139 @*/
DMCompositeGetNumberDM(DM dm,PetscInt * nDM)140 PetscErrorCode DMCompositeGetNumberDM(DM dm,PetscInt *nDM)
141 {
142 DM_Composite *com = (DM_Composite*)dm->data;
143 PetscBool flg;
144 PetscErrorCode ierr;
145
146 PetscFunctionBegin;
147 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
148 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
149 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
150 *nDM = com->nDM;
151 PetscFunctionReturn(0);
152 }
153
154
155 /*@C
156 DMCompositeGetAccess - Allows one to access the individual packed vectors in their global
157 representation.
158
159 Collective on dm
160
161 Input Parameters:
162 + dm - the packer object
163 - gvec - the global vector
164
165 Output Parameters:
166 . Vec* ... - the packed parallel vectors, NULL for those that are not needed
167
168 Notes:
169 Use DMCompositeRestoreAccess() to return the vectors when you no longer need them
170
171 Fortran Notes:
172
173 Fortran callers must use numbered versions of this routine, e.g., DMCompositeGetAccess4(dm,gvec,vec1,vec2,vec3,vec4)
174 or use the alternative interface DMCompositeGetAccessArray().
175
176 Level: advanced
177
178 .seealso: DMCompositeGetEntries(), DMCompositeScatter()
179 @*/
DMCompositeGetAccess(DM dm,Vec gvec,...)180 PetscErrorCode DMCompositeGetAccess(DM dm,Vec gvec,...)
181 {
182 va_list Argp;
183 PetscErrorCode ierr;
184 struct DMCompositeLink *next;
185 DM_Composite *com = (DM_Composite*)dm->data;
186 PetscInt readonly;
187 PetscBool flg;
188
189 PetscFunctionBegin;
190 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
191 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
192 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
193 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
194 next = com->next;
195 if (!com->setup) {
196 ierr = DMSetUp(dm);CHKERRQ(ierr);
197 }
198
199 ierr = VecLockGet(gvec,&readonly);CHKERRQ(ierr);
200 /* loop over packed objects, handling one at at time */
201 va_start(Argp,gvec);
202 while (next) {
203 Vec *vec;
204 vec = va_arg(Argp, Vec*);
205 if (vec) {
206 ierr = DMGetGlobalVector(next->dm,vec);CHKERRQ(ierr);
207 if (readonly) {
208 const PetscScalar *array;
209 ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr);
210 ierr = VecPlaceArray(*vec,array+next->rstart);CHKERRQ(ierr);
211 ierr = VecLockReadPush(*vec);CHKERRQ(ierr);
212 ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr);
213 } else {
214 PetscScalar *array;
215 ierr = VecGetArray(gvec,&array);CHKERRQ(ierr);
216 ierr = VecPlaceArray(*vec,array+next->rstart);CHKERRQ(ierr);
217 ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr);
218 }
219 }
220 next = next->next;
221 }
222 va_end(Argp);
223 PetscFunctionReturn(0);
224 }
225
226 /*@C
227 DMCompositeGetAccessArray - Allows one to access the individual packed vectors in their global
228 representation.
229
230 Collective on dm
231
232 Input Parameters:
233 + dm - the packer object
234 . pvec - packed vector
235 . nwanted - number of vectors wanted
236 - wanted - sorted array of vectors wanted, or NULL to get all vectors
237
238 Output Parameters:
239 . vecs - array of requested global vectors (must be allocated)
240
241 Notes:
242 Use DMCompositeRestoreAccessArray() to return the vectors when you no longer need them
243
244 Level: advanced
245
246 .seealso: DMCompositeGetAccess(), DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
247 @*/
DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt * wanted,Vec * vecs)248 PetscErrorCode DMCompositeGetAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
249 {
250 PetscErrorCode ierr;
251 struct DMCompositeLink *link;
252 PetscInt i,wnum;
253 DM_Composite *com = (DM_Composite*)dm->data;
254 PetscInt readonly;
255 PetscBool flg;
256
257 PetscFunctionBegin;
258 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
259 PetscValidHeaderSpecific(pvec,VEC_CLASSID,2);
260 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
261 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
262 if (!com->setup) {
263 ierr = DMSetUp(dm);CHKERRQ(ierr);
264 }
265
266 ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr);
267 for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
268 if (!wanted || i == wanted[wnum]) {
269 Vec v;
270 ierr = DMGetGlobalVector(link->dm,&v);CHKERRQ(ierr);
271 if (readonly) {
272 const PetscScalar *array;
273 ierr = VecGetArrayRead(pvec,&array);CHKERRQ(ierr);
274 ierr = VecPlaceArray(v,array+link->rstart);CHKERRQ(ierr);
275 ierr = VecLockReadPush(v);CHKERRQ(ierr);
276 ierr = VecRestoreArrayRead(pvec,&array);CHKERRQ(ierr);
277 } else {
278 PetscScalar *array;
279 ierr = VecGetArray(pvec,&array);CHKERRQ(ierr);
280 ierr = VecPlaceArray(v,array+link->rstart);CHKERRQ(ierr);
281 ierr = VecRestoreArray(pvec,&array);CHKERRQ(ierr);
282 }
283 vecs[wnum++] = v;
284 }
285 }
286 PetscFunctionReturn(0);
287 }
288
289 /*@C
290 DMCompositeGetLocalAccessArray - Allows one to access the individual
291 packed vectors in their local representation.
292
293 Collective on dm.
294
295 Input Parameters:
296 + dm - the packer object
297 . pvec - packed vector
298 . nwanted - number of vectors wanted
299 - wanted - sorted array of vectors wanted, or NULL to get all vectors
300
301 Output Parameters:
302 . vecs - array of requested local vectors (must be allocated)
303
304 Notes:
305 Use DMCompositeRestoreLocalAccessArray() to return the vectors
306 when you no longer need them.
307
308 Level: advanced
309
310 .seealso: DMCompositeRestoreLocalAccessArray(), DMCompositeGetAccess(),
311 DMCompositeGetEntries(), DMCompositeScatter(), DMCompositeGather()
312 @*/
DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt * wanted,Vec * vecs)313 PetscErrorCode DMCompositeGetLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
314 {
315 PetscErrorCode ierr;
316 struct DMCompositeLink *link;
317 PetscInt i,wnum;
318 DM_Composite *com = (DM_Composite*)dm->data;
319 PetscInt readonly;
320 PetscInt nlocal = 0;
321 PetscBool flg;
322
323 PetscFunctionBegin;
324 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
325 PetscValidHeaderSpecific(pvec,VEC_CLASSID,2);
326 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
327 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
328 if (!com->setup) {
329 ierr = DMSetUp(dm);CHKERRQ(ierr);
330 }
331
332 ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr);
333 for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
334 if (!wanted || i == wanted[wnum]) {
335 Vec v;
336 ierr = DMGetLocalVector(link->dm,&v);CHKERRQ(ierr);
337 if (readonly) {
338 const PetscScalar *array;
339 ierr = VecGetArrayRead(pvec,&array);CHKERRQ(ierr);
340 ierr = VecPlaceArray(v,array+nlocal);CHKERRQ(ierr);
341 ierr = VecLockReadPush(v);CHKERRQ(ierr);
342 ierr = VecRestoreArrayRead(pvec,&array);CHKERRQ(ierr);
343 } else {
344 PetscScalar *array;
345 ierr = VecGetArray(pvec,&array);CHKERRQ(ierr);
346 ierr = VecPlaceArray(v,array+nlocal);CHKERRQ(ierr);
347 ierr = VecRestoreArray(pvec,&array);CHKERRQ(ierr);
348 }
349 vecs[wnum++] = v;
350 }
351
352 nlocal += link->nlocal;
353 }
354
355 PetscFunctionReturn(0);
356 }
357
358 /*@C
359 DMCompositeRestoreAccess - Returns the vectors obtained with DMCompositeGetAccess()
360 representation.
361
362 Collective on dm
363
364 Input Parameters:
365 + dm - the packer object
366 . gvec - the global vector
367 - Vec* ... - the individual parallel vectors, NULL for those that are not needed
368
369 Level: advanced
370
371 .seealso DMCompositeAddDM(), DMCreateGlobalVector(),
372 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeScatter(),
373 DMCompositeRestoreAccess(), DMCompositeGetAccess()
374
375 @*/
DMCompositeRestoreAccess(DM dm,Vec gvec,...)376 PetscErrorCode DMCompositeRestoreAccess(DM dm,Vec gvec,...)
377 {
378 va_list Argp;
379 PetscErrorCode ierr;
380 struct DMCompositeLink *next;
381 DM_Composite *com = (DM_Composite*)dm->data;
382 PetscInt readonly;
383 PetscBool flg;
384
385 PetscFunctionBegin;
386 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
387 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
388 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
389 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
390 next = com->next;
391 if (!com->setup) {
392 ierr = DMSetUp(dm);CHKERRQ(ierr);
393 }
394
395 ierr = VecLockGet(gvec,&readonly);CHKERRQ(ierr);
396 /* loop over packed objects, handling one at at time */
397 va_start(Argp,gvec);
398 while (next) {
399 Vec *vec;
400 vec = va_arg(Argp, Vec*);
401 if (vec) {
402 ierr = VecResetArray(*vec);CHKERRQ(ierr);
403 if (readonly) {
404 ierr = VecLockReadPop(*vec);CHKERRQ(ierr);
405 }
406 ierr = DMRestoreGlobalVector(next->dm,vec);CHKERRQ(ierr);
407 }
408 next = next->next;
409 }
410 va_end(Argp);
411 PetscFunctionReturn(0);
412 }
413
414 /*@C
415 DMCompositeRestoreAccessArray - Returns the vectors obtained with DMCompositeGetAccessArray()
416
417 Collective on dm
418
419 Input Parameters:
420 + dm - the packer object
421 . pvec - packed vector
422 . nwanted - number of vectors wanted
423 . wanted - sorted array of vectors wanted, or NULL to get all vectors
424 - vecs - array of global vectors to return
425
426 Level: advanced
427
428 .seealso: DMCompositeRestoreAccess(), DMCompositeRestoreEntries(), DMCompositeScatter(), DMCompositeGather()
429 @*/
DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt * wanted,Vec * vecs)430 PetscErrorCode DMCompositeRestoreAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
431 {
432 PetscErrorCode ierr;
433 struct DMCompositeLink *link;
434 PetscInt i,wnum;
435 DM_Composite *com = (DM_Composite*)dm->data;
436 PetscInt readonly;
437 PetscBool flg;
438
439 PetscFunctionBegin;
440 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
441 PetscValidHeaderSpecific(pvec,VEC_CLASSID,2);
442 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
443 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
444 if (!com->setup) {
445 ierr = DMSetUp(dm);CHKERRQ(ierr);
446 }
447
448 ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr);
449 for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
450 if (!wanted || i == wanted[wnum]) {
451 ierr = VecResetArray(vecs[wnum]);CHKERRQ(ierr);
452 if (readonly) {
453 ierr = VecLockReadPop(vecs[wnum]);CHKERRQ(ierr);
454 }
455 ierr = DMRestoreGlobalVector(link->dm,&vecs[wnum]);CHKERRQ(ierr);
456 wnum++;
457 }
458 }
459 PetscFunctionReturn(0);
460 }
461
462 /*@C
463 DMCompositeRestoreLocalAccessArray - Returns the vectors obtained with DMCompositeGetLocalAccessArray().
464
465 Collective on dm.
466
467 Input Parameters:
468 + dm - the packer object
469 . pvec - packed vector
470 . nwanted - number of vectors wanted
471 . wanted - sorted array of vectors wanted, or NULL to restore all vectors
472 - vecs - array of local vectors to return
473
474 Level: advanced
475
476 Notes:
477 nwanted and wanted must match the values given to DMCompositeGetLocalAccessArray()
478 otherwise the call will fail.
479
480 .seealso: DMCompositeGetLocalAccessArray(), DMCompositeRestoreAccessArray(),
481 DMCompositeRestoreAccess(), DMCompositeRestoreEntries(),
482 DMCompositeScatter(), DMCompositeGather()
483 @*/
DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt * wanted,Vec * vecs)484 PetscErrorCode DMCompositeRestoreLocalAccessArray(DM dm,Vec pvec,PetscInt nwanted,const PetscInt *wanted,Vec *vecs)
485 {
486 PetscErrorCode ierr;
487 struct DMCompositeLink *link;
488 PetscInt i,wnum;
489 DM_Composite *com = (DM_Composite*)dm->data;
490 PetscInt readonly;
491 PetscBool flg;
492
493 PetscFunctionBegin;
494 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
495 PetscValidHeaderSpecific(pvec,VEC_CLASSID,2);
496 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
497 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
498 if (!com->setup) {
499 ierr = DMSetUp(dm);CHKERRQ(ierr);
500 }
501
502 ierr = VecLockGet(pvec,&readonly);CHKERRQ(ierr);
503 for (i=0,wnum=0,link=com->next; link && wnum<nwanted; i++,link=link->next) {
504 if (!wanted || i == wanted[wnum]) {
505 ierr = VecResetArray(vecs[wnum]);CHKERRQ(ierr);
506 if (readonly) {
507 ierr = VecLockReadPop(vecs[wnum]);CHKERRQ(ierr);
508 }
509 ierr = DMRestoreLocalVector(link->dm,&vecs[wnum]);CHKERRQ(ierr);
510 wnum++;
511 }
512 }
513 PetscFunctionReturn(0);
514 }
515
516 /*@C
517 DMCompositeScatter - Scatters from a global packed vector into its individual local vectors
518
519 Collective on dm
520
521 Input Parameters:
522 + dm - the packer object
523 . gvec - the global vector
524 - Vec ... - the individual sequential vectors, NULL for those that are not needed
525
526 Level: advanced
527
528 Notes:
529 DMCompositeScatterArray() is a non-variadic alternative that is often more convenient for library callers and is
530 accessible from Fortran.
531
532 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
533 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
534 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
535 DMCompositeScatterArray()
536
537 @*/
DMCompositeScatter(DM dm,Vec gvec,...)538 PetscErrorCode DMCompositeScatter(DM dm,Vec gvec,...)
539 {
540 va_list Argp;
541 PetscErrorCode ierr;
542 struct DMCompositeLink *next;
543 PetscInt cnt;
544 DM_Composite *com = (DM_Composite*)dm->data;
545 PetscBool flg;
546
547 PetscFunctionBegin;
548 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
549 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
550 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
551 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
552 if (!com->setup) {
553 ierr = DMSetUp(dm);CHKERRQ(ierr);
554 }
555
556 /* loop over packed objects, handling one at at time */
557 va_start(Argp,gvec);
558 for (cnt=3,next=com->next; next; cnt++,next=next->next) {
559 Vec local;
560 local = va_arg(Argp, Vec);
561 if (local) {
562 Vec global;
563 const PetscScalar *array;
564 PetscValidHeaderSpecific(local,VEC_CLASSID,cnt);
565 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
566 ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr);
567 ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr);
568 ierr = DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,local);CHKERRQ(ierr);
569 ierr = DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,local);CHKERRQ(ierr);
570 ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr);
571 ierr = VecResetArray(global);CHKERRQ(ierr);
572 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
573 }
574 }
575 va_end(Argp);
576 PetscFunctionReturn(0);
577 }
578
579 /*@
580 DMCompositeScatterArray - Scatters from a global packed vector into its individual local vectors
581
582 Collective on dm
583
584 Input Parameters:
585 + dm - the packer object
586 . gvec - the global vector
587 - lvecs - array of local vectors, NULL for any that are not needed
588
589 Level: advanced
590
591 Note:
592 This is a non-variadic alternative to DMCompositeScatter()
593
594 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector()
595 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
596 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
597
598 @*/
DMCompositeScatterArray(DM dm,Vec gvec,Vec * lvecs)599 PetscErrorCode DMCompositeScatterArray(DM dm,Vec gvec,Vec *lvecs)
600 {
601 PetscErrorCode ierr;
602 struct DMCompositeLink *next;
603 PetscInt i;
604 DM_Composite *com = (DM_Composite*)dm->data;
605 PetscBool flg;
606
607 PetscFunctionBegin;
608 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
609 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
610 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
611 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
612 if (!com->setup) {
613 ierr = DMSetUp(dm);CHKERRQ(ierr);
614 }
615
616 /* loop over packed objects, handling one at at time */
617 for (i=0,next=com->next; next; next=next->next,i++) {
618 if (lvecs[i]) {
619 Vec global;
620 const PetscScalar *array;
621 PetscValidHeaderSpecific(lvecs[i],VEC_CLASSID,3);
622 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
623 ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr);
624 ierr = VecPlaceArray(global,(PetscScalar*)array+next->rstart);CHKERRQ(ierr);
625 ierr = DMGlobalToLocalBegin(next->dm,global,INSERT_VALUES,lvecs[i]);CHKERRQ(ierr);
626 ierr = DMGlobalToLocalEnd(next->dm,global,INSERT_VALUES,lvecs[i]);CHKERRQ(ierr);
627 ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr);
628 ierr = VecResetArray(global);CHKERRQ(ierr);
629 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
630 }
631 }
632 PetscFunctionReturn(0);
633 }
634
635 /*@C
636 DMCompositeGather - Gathers into a global packed vector from its individual local vectors
637
638 Collective on dm
639
640 Input Parameter:
641 + dm - the packer object
642 . gvec - the global vector
643 . imode - INSERT_VALUES or ADD_VALUES
644 - Vec ... - the individual sequential vectors, NULL for any that are not needed
645
646 Level: advanced
647
648 Not available from Fortran, Fortran users can use DMCompositeGatherArray()
649
650 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
651 DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
652 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
653
654 @*/
DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)655 PetscErrorCode DMCompositeGather(DM dm,InsertMode imode,Vec gvec,...)
656 {
657 va_list Argp;
658 PetscErrorCode ierr;
659 struct DMCompositeLink *next;
660 DM_Composite *com = (DM_Composite*)dm->data;
661 PetscInt cnt;
662 PetscBool flg;
663
664 PetscFunctionBegin;
665 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
666 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
667 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
668 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
669 if (!com->setup) {
670 ierr = DMSetUp(dm);CHKERRQ(ierr);
671 }
672
673 /* loop over packed objects, handling one at at time */
674 va_start(Argp,gvec);
675 for (cnt=3,next=com->next; next; cnt++,next=next->next) {
676 Vec local;
677 local = va_arg(Argp, Vec);
678 if (local) {
679 PetscScalar *array;
680 Vec global;
681 PetscValidHeaderSpecific(local,VEC_CLASSID,cnt);
682 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
683 ierr = VecGetArray(gvec,&array);CHKERRQ(ierr);
684 ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr);
685 ierr = DMLocalToGlobalBegin(next->dm,local,imode,global);CHKERRQ(ierr);
686 ierr = DMLocalToGlobalEnd(next->dm,local,imode,global);CHKERRQ(ierr);
687 ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr);
688 ierr = VecResetArray(global);CHKERRQ(ierr);
689 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
690 }
691 }
692 va_end(Argp);
693 PetscFunctionReturn(0);
694 }
695
696 /*@
697 DMCompositeGatherArray - Gathers into a global packed vector from its individual local vectors
698
699 Collective on dm
700
701 Input Parameter:
702 + dm - the packer object
703 . gvec - the global vector
704 . imode - INSERT_VALUES or ADD_VALUES
705 - lvecs - the individual sequential vectors, NULL for any that are not needed
706
707 Level: advanced
708
709 Notes:
710 This is a non-variadic alternative to DMCompositeGather().
711
712 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
713 DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
714 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries(),
715 @*/
DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec * lvecs)716 PetscErrorCode DMCompositeGatherArray(DM dm,InsertMode imode,Vec gvec,Vec *lvecs)
717 {
718 PetscErrorCode ierr;
719 struct DMCompositeLink *next;
720 DM_Composite *com = (DM_Composite*)dm->data;
721 PetscInt i;
722 PetscBool flg;
723
724 PetscFunctionBegin;
725 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
726 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
727 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
728 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
729 if (!com->setup) {
730 ierr = DMSetUp(dm);CHKERRQ(ierr);
731 }
732
733 /* loop over packed objects, handling one at at time */
734 for (next=com->next,i=0; next; next=next->next,i++) {
735 if (lvecs[i]) {
736 PetscScalar *array;
737 Vec global;
738 PetscValidHeaderSpecific(lvecs[i],VEC_CLASSID,3);
739 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
740 ierr = VecGetArray(gvec,&array);CHKERRQ(ierr);
741 ierr = VecPlaceArray(global,array+next->rstart);CHKERRQ(ierr);
742 ierr = DMLocalToGlobalBegin(next->dm,lvecs[i],imode,global);CHKERRQ(ierr);
743 ierr = DMLocalToGlobalEnd(next->dm,lvecs[i],imode,global);CHKERRQ(ierr);
744 ierr = VecRestoreArray(gvec,&array);CHKERRQ(ierr);
745 ierr = VecResetArray(global);CHKERRQ(ierr);
746 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
747 }
748 }
749 PetscFunctionReturn(0);
750 }
751
752 /*@
753 DMCompositeAddDM - adds a DM vector to a DMComposite
754
755 Collective on dm
756
757 Input Parameter:
758 + dmc - the DMComposite (packer) object
759 - dm - the DM object
760
761 Level: advanced
762
763 .seealso DMDestroy(), DMCompositeGather(), DMCompositeAddDM(), DMCreateGlobalVector(),
764 DMCompositeScatter(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
765 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
766
767 @*/
DMCompositeAddDM(DM dmc,DM dm)768 PetscErrorCode DMCompositeAddDM(DM dmc,DM dm)
769 {
770 PetscErrorCode ierr;
771 PetscInt n,nlocal;
772 struct DMCompositeLink *mine,*next;
773 Vec global,local;
774 DM_Composite *com = (DM_Composite*)dmc->data;
775 PetscBool flg;
776
777 PetscFunctionBegin;
778 PetscValidHeaderSpecific(dmc,DM_CLASSID,1);
779 PetscValidHeaderSpecific(dm,DM_CLASSID,2);
780 ierr = PetscObjectTypeCompare((PetscObject)dmc,DMCOMPOSITE,&flg);CHKERRQ(ierr);
781 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
782 next = com->next;
783 if (com->setup) SETERRQ(PetscObjectComm((PetscObject)dmc),PETSC_ERR_ARG_WRONGSTATE,"Cannot add a DM once you have used the DMComposite");
784
785 /* create new link */
786 ierr = PetscNew(&mine);CHKERRQ(ierr);
787 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);
788 ierr = DMGetGlobalVector(dm,&global);CHKERRQ(ierr);
789 ierr = VecGetLocalSize(global,&n);CHKERRQ(ierr);
790 ierr = DMRestoreGlobalVector(dm,&global);CHKERRQ(ierr);
791 ierr = DMGetLocalVector(dm,&local);CHKERRQ(ierr);
792 ierr = VecGetSize(local,&nlocal);CHKERRQ(ierr);
793 ierr = DMRestoreLocalVector(dm,&local);CHKERRQ(ierr);
794
795 mine->n = n;
796 mine->nlocal = nlocal;
797 mine->dm = dm;
798 mine->next = NULL;
799 com->n += n;
800 com->nghost += nlocal;
801
802 /* add to end of list */
803 if (!next) com->next = mine;
804 else {
805 while (next->next) next = next->next;
806 next->next = mine;
807 }
808 com->nDM++;
809 com->nmine++;
810 PetscFunctionReturn(0);
811 }
812
813 #include <petscdraw.h>
814 PETSC_EXTERN PetscErrorCode VecView_MPI(Vec,PetscViewer);
VecView_DMComposite(Vec gvec,PetscViewer viewer)815 PetscErrorCode VecView_DMComposite(Vec gvec,PetscViewer viewer)
816 {
817 DM dm;
818 PetscErrorCode ierr;
819 struct DMCompositeLink *next;
820 PetscBool isdraw;
821 DM_Composite *com;
822
823 PetscFunctionBegin;
824 ierr = VecGetDM(gvec, &dm);CHKERRQ(ierr);
825 if (!dm) SETERRQ(PetscObjectComm((PetscObject)gvec),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMComposite");
826 com = (DM_Composite*)dm->data;
827 next = com->next;
828
829 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
830 if (!isdraw) {
831 /* do I really want to call this? */
832 ierr = VecView_MPI(gvec,viewer);CHKERRQ(ierr);
833 } else {
834 PetscInt cnt = 0;
835
836 /* loop over packed objects, handling one at at time */
837 while (next) {
838 Vec vec;
839 const PetscScalar *array;
840 PetscInt bs;
841
842 /* Should use VecGetSubVector() eventually, but would need to forward the DM for that to work */
843 ierr = DMGetGlobalVector(next->dm,&vec);CHKERRQ(ierr);
844 ierr = VecGetArrayRead(gvec,&array);CHKERRQ(ierr);
845 ierr = VecPlaceArray(vec,(PetscScalar*)array+next->rstart);CHKERRQ(ierr);
846 ierr = VecRestoreArrayRead(gvec,&array);CHKERRQ(ierr);
847 ierr = VecView(vec,viewer);CHKERRQ(ierr);
848 ierr = VecResetArray(vec);CHKERRQ(ierr);
849 ierr = VecGetBlockSize(vec,&bs);CHKERRQ(ierr);
850 ierr = DMRestoreGlobalVector(next->dm,&vec);CHKERRQ(ierr);
851 ierr = PetscViewerDrawBaseAdd(viewer,bs);CHKERRQ(ierr);
852 cnt += bs;
853 next = next->next;
854 }
855 ierr = PetscViewerDrawBaseAdd(viewer,-cnt);CHKERRQ(ierr);
856 }
857 PetscFunctionReturn(0);
858 }
859
DMCreateGlobalVector_Composite(DM dm,Vec * gvec)860 PetscErrorCode DMCreateGlobalVector_Composite(DM dm,Vec *gvec)
861 {
862 PetscErrorCode ierr;
863 DM_Composite *com = (DM_Composite*)dm->data;
864
865 PetscFunctionBegin;
866 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
867 ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
868 ierr = DMSetUp(dm);CHKERRQ(ierr);
869 ierr = VecCreate(PetscObjectComm((PetscObject)dm),gvec);CHKERRQ(ierr);
870 ierr = VecSetType(*gvec,dm->vectype);CHKERRQ(ierr);
871 ierr = VecSetSizes(*gvec,com->n,com->N);CHKERRQ(ierr);
872 ierr = VecSetDM(*gvec, dm);CHKERRQ(ierr);
873 ierr = VecSetOperation(*gvec,VECOP_VIEW,(void (*)(void))VecView_DMComposite);CHKERRQ(ierr);
874 PetscFunctionReturn(0);
875 }
876
DMCreateLocalVector_Composite(DM dm,Vec * lvec)877 PetscErrorCode DMCreateLocalVector_Composite(DM dm,Vec *lvec)
878 {
879 PetscErrorCode ierr;
880 DM_Composite *com = (DM_Composite*)dm->data;
881
882 PetscFunctionBegin;
883 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
884 if (!com->setup) {
885 ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
886 ierr = DMSetUp(dm);CHKERRQ(ierr);
887 }
888 ierr = VecCreate(PETSC_COMM_SELF,lvec);CHKERRQ(ierr);
889 ierr = VecSetType(*lvec,dm->vectype);CHKERRQ(ierr);
890 ierr = VecSetSizes(*lvec,com->nghost,PETSC_DECIDE);CHKERRQ(ierr);
891 ierr = VecSetDM(*lvec, dm);CHKERRQ(ierr);
892 PetscFunctionReturn(0);
893 }
894
895 /*@C
896 DMCompositeGetISLocalToGlobalMappings - gets an ISLocalToGlobalMapping for each DM in the DMComposite, maps to the composite global space
897
898 Collective on DM
899
900 Input Parameter:
901 . dm - the packer object
902
903 Output Parameters:
904 . ltogs - the individual mappings for each packed vector. Note that this includes
905 all the ghost points that individual ghosted DMDA's may have.
906
907 Level: advanced
908
909 Notes:
910 Each entry of ltogs should be destroyed with ISLocalToGlobalMappingDestroy(), the ltogs array should be freed with PetscFree().
911
912 Not available from Fortran
913
914 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
915 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
916 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
917
918 @*/
DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping ** ltogs)919 PetscErrorCode DMCompositeGetISLocalToGlobalMappings(DM dm,ISLocalToGlobalMapping **ltogs)
920 {
921 PetscErrorCode ierr;
922 PetscInt i,*idx,n,cnt;
923 struct DMCompositeLink *next;
924 PetscMPIInt rank;
925 DM_Composite *com = (DM_Composite*)dm->data;
926 PetscBool flg;
927
928 PetscFunctionBegin;
929 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
930 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
931 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
932 ierr = DMSetUp(dm);CHKERRQ(ierr);
933 ierr = PetscMalloc1(com->nDM,ltogs);CHKERRQ(ierr);
934 next = com->next;
935 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
936
937 /* loop over packed objects, handling one at at time */
938 cnt = 0;
939 while (next) {
940 ISLocalToGlobalMapping ltog;
941 PetscMPIInt size;
942 const PetscInt *suboff,*indices;
943 Vec global;
944
945 /* Get sub-DM global indices for each local dof */
946 ierr = DMGetLocalToGlobalMapping(next->dm,<og);CHKERRQ(ierr);
947 ierr = ISLocalToGlobalMappingGetSize(ltog,&n);CHKERRQ(ierr);
948 ierr = ISLocalToGlobalMappingGetIndices(ltog,&indices);CHKERRQ(ierr);
949 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr);
950
951 /* Get the offsets for the sub-DM global vector */
952 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
953 ierr = VecGetOwnershipRanges(global,&suboff);CHKERRQ(ierr);
954 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)global),&size);CHKERRQ(ierr);
955
956 /* Shift the sub-DM definition of the global space to the composite global space */
957 for (i=0; i<n; i++) {
958 PetscInt subi = indices[i],lo = 0,hi = size,t;
959 /* There's no consensus on what a negative index means,
960 except for skipping when setting the values in vectors and matrices */
961 if (subi < 0) { idx[i] = subi - next->grstarts[rank]; continue; }
962 /* Binary search to find which rank owns subi */
963 while (hi-lo > 1) {
964 t = lo + (hi-lo)/2;
965 if (suboff[t] > subi) hi = t;
966 else lo = t;
967 }
968 idx[i] = subi - suboff[lo] + next->grstarts[lo];
969 }
970 ierr = ISLocalToGlobalMappingRestoreIndices(ltog,&indices);CHKERRQ(ierr);
971 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm),1,n,idx,PETSC_OWN_POINTER,&(*ltogs)[cnt]);CHKERRQ(ierr);
972 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
973 next = next->next;
974 cnt++;
975 }
976 PetscFunctionReturn(0);
977 }
978
979 /*@C
980 DMCompositeGetLocalISs - Gets index sets for each component of a composite local vector
981
982 Not Collective
983
984 Input Arguments:
985 . dm - composite DM
986
987 Output Arguments:
988 . is - array of serial index sets for each each component of the DMComposite
989
990 Level: intermediate
991
992 Notes:
993 At present, a composite local vector does not normally exist. This function is used to provide index sets for
994 MatGetLocalSubMatrix(). In the future, the scatters for each entry in the DMComposite may be be merged into a single
995 scatter to a composite local vector. The user should not typically need to know which is being done.
996
997 To get the composite global indices at all local points (including ghosts), use DMCompositeGetISLocalToGlobalMappings().
998
999 To get index sets for pieces of the composite global vector, use DMCompositeGetGlobalISs().
1000
1001 Each returned IS should be destroyed with ISDestroy(), the array should be freed with PetscFree().
1002
1003 Not available from Fortran
1004
1005 .seealso: DMCompositeGetGlobalISs(), DMCompositeGetISLocalToGlobalMappings(), MatGetLocalSubMatrix(), MatCreateLocalRef()
1006 @*/
DMCompositeGetLocalISs(DM dm,IS ** is)1007 PetscErrorCode DMCompositeGetLocalISs(DM dm,IS **is)
1008 {
1009 PetscErrorCode ierr;
1010 DM_Composite *com = (DM_Composite*)dm->data;
1011 struct DMCompositeLink *link;
1012 PetscInt cnt,start;
1013 PetscBool flg;
1014
1015 PetscFunctionBegin;
1016 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1017 PetscValidPointer(is,2);
1018 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1019 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1020 ierr = PetscMalloc1(com->nmine,is);CHKERRQ(ierr);
1021 for (cnt=0,start=0,link=com->next; link; start+=link->nlocal,cnt++,link=link->next) {
1022 PetscInt bs;
1023 ierr = ISCreateStride(PETSC_COMM_SELF,link->nlocal,start,1,&(*is)[cnt]);CHKERRQ(ierr);
1024 ierr = DMGetBlockSize(link->dm,&bs);CHKERRQ(ierr);
1025 ierr = ISSetBlockSize((*is)[cnt],bs);CHKERRQ(ierr);
1026 }
1027 PetscFunctionReturn(0);
1028 }
1029
1030 /*@C
1031 DMCompositeGetGlobalISs - Gets the index sets for each composed object
1032
1033 Collective on dm
1034
1035 Input Parameter:
1036 . dm - the packer object
1037
1038 Output Parameters:
1039 . is - the array of index sets
1040
1041 Level: advanced
1042
1043 Notes:
1044 The is entries should be destroyed with ISDestroy(), the is array should be freed with PetscFree()
1045
1046 These could be used to extract a subset of vector entries for a "multi-physics" preconditioner
1047
1048 Use DMCompositeGetLocalISs() for index sets in the packed local numbering, and
1049 DMCompositeGetISLocalToGlobalMappings() for to map local sub-DM (including ghost) indices to packed global
1050 indices.
1051
1052 Fortran Notes:
1053
1054 The output argument 'is' must be an allocated array of sufficient length, which can be learned using DMCompositeGetNumberDM().
1055
1056 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1057 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetAccess(), DMCompositeScatter(),
1058 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(),DMCompositeGetEntries()
1059
1060 @*/
DMCompositeGetGlobalISs(DM dm,IS * is[])1061 PetscErrorCode DMCompositeGetGlobalISs(DM dm,IS *is[])
1062 {
1063 PetscErrorCode ierr;
1064 PetscInt cnt = 0;
1065 struct DMCompositeLink *next;
1066 PetscMPIInt rank;
1067 DM_Composite *com = (DM_Composite*)dm->data;
1068 PetscBool flg;
1069
1070 PetscFunctionBegin;
1071 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1072 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1073 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1074 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Must call DMSetUp() before");
1075 ierr = PetscMalloc1(com->nDM,is);CHKERRQ(ierr);
1076 next = com->next;
1077 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
1078
1079 /* loop over packed objects, handling one at at time */
1080 while (next) {
1081 PetscDS prob;
1082
1083 ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),next->n,next->grstart,1,&(*is)[cnt]);CHKERRQ(ierr);
1084 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
1085 if (prob) {
1086 MatNullSpace space;
1087 Mat pmat;
1088 PetscObject disc;
1089 PetscInt Nf;
1090
1091 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr);
1092 if (cnt < Nf) {
1093 ierr = PetscDSGetDiscretization(prob, cnt, &disc);CHKERRQ(ierr);
1094 ierr = PetscObjectQuery(disc, "nullspace", (PetscObject*) &space);CHKERRQ(ierr);
1095 if (space) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "nullspace", (PetscObject) space);CHKERRQ(ierr);}
1096 ierr = PetscObjectQuery(disc, "nearnullspace", (PetscObject*) &space);CHKERRQ(ierr);
1097 if (space) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "nearnullspace", (PetscObject) space);CHKERRQ(ierr);}
1098 ierr = PetscObjectQuery(disc, "pmat", (PetscObject*) &pmat);CHKERRQ(ierr);
1099 if (pmat) {ierr = PetscObjectCompose((PetscObject) (*is)[cnt], "pmat", (PetscObject) pmat);CHKERRQ(ierr);}
1100 }
1101 }
1102 cnt++;
1103 next = next->next;
1104 }
1105 PetscFunctionReturn(0);
1106 }
1107
DMCreateFieldIS_Composite(DM dm,PetscInt * numFields,char *** fieldNames,IS ** fields)1108 PetscErrorCode DMCreateFieldIS_Composite(DM dm, PetscInt *numFields,char ***fieldNames, IS **fields)
1109 {
1110 PetscInt nDM;
1111 DM *dms;
1112 PetscInt i;
1113 PetscErrorCode ierr;
1114
1115 PetscFunctionBegin;
1116 ierr = DMCompositeGetNumberDM(dm, &nDM);CHKERRQ(ierr);
1117 if (numFields) *numFields = nDM;
1118 ierr = DMCompositeGetGlobalISs(dm, fields);CHKERRQ(ierr);
1119 if (fieldNames) {
1120 ierr = PetscMalloc1(nDM, &dms);CHKERRQ(ierr);
1121 ierr = PetscMalloc1(nDM, fieldNames);CHKERRQ(ierr);
1122 ierr = DMCompositeGetEntriesArray(dm, dms);CHKERRQ(ierr);
1123 for (i=0; i<nDM; i++) {
1124 char buf[256];
1125 const char *splitname;
1126
1127 /* Split naming precedence: object name, prefix, number */
1128 splitname = ((PetscObject) dm)->name;
1129 if (!splitname) {
1130 ierr = PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);CHKERRQ(ierr);
1131 if (splitname) {
1132 size_t len;
1133 ierr = PetscStrncpy(buf,splitname,sizeof(buf));CHKERRQ(ierr);
1134 buf[sizeof(buf) - 1] = 0;
1135 ierr = PetscStrlen(buf,&len);CHKERRQ(ierr);
1136 if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */
1137 splitname = buf;
1138 }
1139 }
1140 if (!splitname) {
1141 ierr = PetscSNPrintf(buf,sizeof(buf),"%D",i);CHKERRQ(ierr);
1142 splitname = buf;
1143 }
1144 ierr = PetscStrallocpy(splitname,&(*fieldNames)[i]);CHKERRQ(ierr);
1145 }
1146 ierr = PetscFree(dms);CHKERRQ(ierr);
1147 }
1148 PetscFunctionReturn(0);
1149 }
1150
1151 /*
1152 This could take over from DMCreateFieldIS(), as it is more general,
1153 making DMCreateFieldIS() a special case -- calling with dmlist == NULL;
1154 At this point it's probably best to be less intrusive, however.
1155 */
DMCreateFieldDecomposition_Composite(DM dm,PetscInt * len,char *** namelist,IS ** islist,DM ** dmlist)1156 PetscErrorCode DMCreateFieldDecomposition_Composite(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist)
1157 {
1158 PetscInt nDM;
1159 PetscInt i;
1160 PetscErrorCode ierr;
1161
1162 PetscFunctionBegin;
1163 ierr = DMCreateFieldIS_Composite(dm, len, namelist, islist);CHKERRQ(ierr);
1164 if (dmlist) {
1165 ierr = DMCompositeGetNumberDM(dm, &nDM);CHKERRQ(ierr);
1166 ierr = PetscMalloc1(nDM, dmlist);CHKERRQ(ierr);
1167 ierr = DMCompositeGetEntriesArray(dm, *dmlist);CHKERRQ(ierr);
1168 for (i=0; i<nDM; i++) {
1169 ierr = PetscObjectReference((PetscObject)((*dmlist)[i]));CHKERRQ(ierr);
1170 }
1171 }
1172 PetscFunctionReturn(0);
1173 }
1174
1175
1176
1177 /* -------------------------------------------------------------------------------------*/
1178 /*@C
1179 DMCompositeGetLocalVectors - Gets local vectors for each part of a DMComposite.
1180 Use DMCompositeRestoreLocalVectors() to return them.
1181
1182 Not Collective
1183
1184 Input Parameter:
1185 . dm - the packer object
1186
1187 Output Parameter:
1188 . Vec ... - the individual sequential Vecs
1189
1190 Level: advanced
1191
1192 Not available from Fortran
1193
1194 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1195 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1196 DMCompositeRestoreLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1197
1198 @*/
DMCompositeGetLocalVectors(DM dm,...)1199 PetscErrorCode DMCompositeGetLocalVectors(DM dm,...)
1200 {
1201 va_list Argp;
1202 PetscErrorCode ierr;
1203 struct DMCompositeLink *next;
1204 DM_Composite *com = (DM_Composite*)dm->data;
1205 PetscBool flg;
1206
1207 PetscFunctionBegin;
1208 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1209 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1210 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1211 next = com->next;
1212 /* loop over packed objects, handling one at at time */
1213 va_start(Argp,dm);
1214 while (next) {
1215 Vec *vec;
1216 vec = va_arg(Argp, Vec*);
1217 if (vec) {ierr = DMGetLocalVector(next->dm,vec);CHKERRQ(ierr);}
1218 next = next->next;
1219 }
1220 va_end(Argp);
1221 PetscFunctionReturn(0);
1222 }
1223
1224 /*@C
1225 DMCompositeRestoreLocalVectors - Restores local vectors for each part of a DMComposite.
1226
1227 Not Collective
1228
1229 Input Parameter:
1230 . dm - the packer object
1231
1232 Output Parameter:
1233 . Vec ... - the individual sequential Vecs
1234
1235 Level: advanced
1236
1237 Not available from Fortran
1238
1239 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(),
1240 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1241 DMCompositeGetLocalVectors(), DMCompositeScatter(), DMCompositeGetEntries()
1242
1243 @*/
DMCompositeRestoreLocalVectors(DM dm,...)1244 PetscErrorCode DMCompositeRestoreLocalVectors(DM dm,...)
1245 {
1246 va_list Argp;
1247 PetscErrorCode ierr;
1248 struct DMCompositeLink *next;
1249 DM_Composite *com = (DM_Composite*)dm->data;
1250 PetscBool flg;
1251
1252 PetscFunctionBegin;
1253 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1254 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1255 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1256 next = com->next;
1257 /* loop over packed objects, handling one at at time */
1258 va_start(Argp,dm);
1259 while (next) {
1260 Vec *vec;
1261 vec = va_arg(Argp, Vec*);
1262 if (vec) {ierr = DMRestoreLocalVector(next->dm,vec);CHKERRQ(ierr);}
1263 next = next->next;
1264 }
1265 va_end(Argp);
1266 PetscFunctionReturn(0);
1267 }
1268
1269 /* -------------------------------------------------------------------------------------*/
1270 /*@C
1271 DMCompositeGetEntries - Gets the DM for each entry in a DMComposite.
1272
1273 Not Collective
1274
1275 Input Parameter:
1276 . dm - the packer object
1277
1278 Output Parameter:
1279 . DM ... - the individual entries (DMs)
1280
1281 Level: advanced
1282
1283 Fortran Notes:
1284 Available as DMCompositeGetEntries() for one output DM, DMCompositeGetEntries2() for 2, etc
1285
1286 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntriesArray()
1287 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1288 DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1289 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1290
1291 @*/
DMCompositeGetEntries(DM dm,...)1292 PetscErrorCode DMCompositeGetEntries(DM dm,...)
1293 {
1294 va_list Argp;
1295 struct DMCompositeLink *next;
1296 DM_Composite *com = (DM_Composite*)dm->data;
1297 PetscBool flg;
1298 PetscErrorCode ierr;
1299
1300 PetscFunctionBegin;
1301 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1302 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1303 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1304 next = com->next;
1305 /* loop over packed objects, handling one at at time */
1306 va_start(Argp,dm);
1307 while (next) {
1308 DM *dmn;
1309 dmn = va_arg(Argp, DM*);
1310 if (dmn) *dmn = next->dm;
1311 next = next->next;
1312 }
1313 va_end(Argp);
1314 PetscFunctionReturn(0);
1315 }
1316
1317 /*@C
1318 DMCompositeGetEntriesArray - Gets the DM for each entry in a DMComposite.
1319
1320 Not Collective
1321
1322 Input Parameter:
1323 . dm - the packer object
1324
1325 Output Parameter:
1326 . dms - array of sufficient length (see DMCompositeGetNumberDM()) to hold the individual DMs
1327
1328 Level: advanced
1329
1330 .seealso DMDestroy(), DMCompositeAddDM(), DMCreateGlobalVector(), DMCompositeGetEntries()
1331 DMCompositeGather(), DMCompositeCreate(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess(),
1332 DMCompositeRestoreLocalVectors(), DMCompositeGetLocalVectors(), DMCompositeScatter(),
1333 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors()
1334
1335 @*/
DMCompositeGetEntriesArray(DM dm,DM dms[])1336 PetscErrorCode DMCompositeGetEntriesArray(DM dm,DM dms[])
1337 {
1338 struct DMCompositeLink *next;
1339 DM_Composite *com = (DM_Composite*)dm->data;
1340 PetscInt i;
1341 PetscBool flg;
1342 PetscErrorCode ierr;
1343
1344 PetscFunctionBegin;
1345 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1346 ierr = PetscObjectTypeCompare((PetscObject)dm,DMCOMPOSITE,&flg);CHKERRQ(ierr);
1347 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for type %s",((PetscObject)dm)->type_name);
1348 /* loop over packed objects, handling one at at time */
1349 for (next=com->next,i=0; next; next=next->next,i++) dms[i] = next->dm;
1350 PetscFunctionReturn(0);
1351 }
1352
1353 typedef struct {
1354 DM dm;
1355 PetscViewer *subv;
1356 Vec *vecs;
1357 } GLVisViewerCtx;
1358
DestroyGLVisViewerCtx_Private(void * vctx)1359 static PetscErrorCode DestroyGLVisViewerCtx_Private(void *vctx)
1360 {
1361 GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1362 PetscInt i,n;
1363 PetscErrorCode ierr;
1364
1365 PetscFunctionBegin;
1366 ierr = DMCompositeGetNumberDM(ctx->dm,&n);CHKERRQ(ierr);
1367 for (i = 0; i < n; i++) {
1368 ierr = PetscViewerDestroy(&ctx->subv[i]);CHKERRQ(ierr);
1369 }
1370 ierr = PetscFree2(ctx->subv,ctx->vecs);CHKERRQ(ierr);
1371 ierr = DMDestroy(&ctx->dm);CHKERRQ(ierr);
1372 ierr = PetscFree(ctx);CHKERRQ(ierr);
1373 PetscFunctionReturn(0);
1374 }
1375
DMCompositeSampleGLVisFields_Private(PetscObject oX,PetscInt nf,PetscObject oXfield[],void * vctx)1376 static PetscErrorCode DMCompositeSampleGLVisFields_Private(PetscObject oX, PetscInt nf, PetscObject oXfield[], void *vctx)
1377 {
1378 Vec X = (Vec)oX;
1379 GLVisViewerCtx *ctx = (GLVisViewerCtx*)vctx;
1380 PetscInt i,n,cumf;
1381 PetscErrorCode ierr;
1382
1383 PetscFunctionBegin;
1384 ierr = DMCompositeGetNumberDM(ctx->dm,&n);CHKERRQ(ierr);
1385 ierr = DMCompositeGetAccessArray(ctx->dm,X,n,NULL,ctx->vecs);CHKERRQ(ierr);
1386 for (i = 0, cumf = 0; i < n; i++) {
1387 PetscErrorCode (*g2l)(PetscObject,PetscInt,PetscObject[],void*);
1388 void *fctx;
1389 PetscInt nfi;
1390
1391 ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nfi,NULL,NULL,&g2l,NULL,&fctx);CHKERRQ(ierr);
1392 if (!nfi) continue;
1393 if (g2l) {
1394 ierr = (*g2l)((PetscObject)ctx->vecs[i],nfi,oXfield+cumf,fctx);CHKERRQ(ierr);
1395 } else {
1396 ierr = VecCopy(ctx->vecs[i],(Vec)(oXfield[cumf]));CHKERRQ(ierr);
1397 }
1398 cumf += nfi;
1399 }
1400 ierr = DMCompositeRestoreAccessArray(ctx->dm,X,n,NULL,ctx->vecs);CHKERRQ(ierr);
1401 PetscFunctionReturn(0);
1402 }
1403
DMSetUpGLVisViewer_Composite(PetscObject odm,PetscViewer viewer)1404 static PetscErrorCode DMSetUpGLVisViewer_Composite(PetscObject odm, PetscViewer viewer)
1405 {
1406 DM dm = (DM)odm, *dms;
1407 Vec *Ufds;
1408 GLVisViewerCtx *ctx;
1409 PetscInt i,n,tnf,*sdim;
1410 char **fecs;
1411 PetscErrorCode ierr;
1412
1413 PetscFunctionBegin;
1414 ierr = PetscNew(&ctx);CHKERRQ(ierr);
1415 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);
1416 ctx->dm = dm;
1417 ierr = DMCompositeGetNumberDM(dm,&n);CHKERRQ(ierr);
1418 ierr = PetscMalloc1(n,&dms);CHKERRQ(ierr);
1419 ierr = DMCompositeGetEntriesArray(dm,dms);CHKERRQ(ierr);
1420 ierr = PetscMalloc2(n,&ctx->subv,n,&ctx->vecs);CHKERRQ(ierr);
1421 for (i = 0, tnf = 0; i < n; i++) {
1422 PetscInt nf;
1423
1424 ierr = PetscViewerCreate(PetscObjectComm(odm),&ctx->subv[i]);CHKERRQ(ierr);
1425 ierr = PetscViewerSetType(ctx->subv[i],PETSCVIEWERGLVIS);CHKERRQ(ierr);
1426 ierr = PetscViewerGLVisSetDM_Private(ctx->subv[i],(PetscObject)dms[i]);CHKERRQ(ierr);
1427 ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
1428 tnf += nf;
1429 }
1430 ierr = PetscFree(dms);CHKERRQ(ierr);
1431 ierr = PetscMalloc3(tnf,&fecs,tnf,&sdim,tnf,&Ufds);CHKERRQ(ierr);
1432 for (i = 0, tnf = 0; i < n; i++) {
1433 PetscInt *sd,nf,f;
1434 const char **fec;
1435 Vec *Uf;
1436
1437 ierr = PetscViewerGLVisGetFields_Private(ctx->subv[i],&nf,&fec,&sd,NULL,(PetscObject**)&Uf,NULL);CHKERRQ(ierr);
1438 for (f = 0; f < nf; f++) {
1439 ierr = PetscStrallocpy(fec[f],&fecs[tnf+f]);CHKERRQ(ierr);
1440 Ufds[tnf+f] = Uf[f];
1441 sdim[tnf+f] = sd[f];
1442 }
1443 tnf += nf;
1444 }
1445 ierr = PetscViewerGLVisSetFields(viewer,tnf,(const char**)fecs,sdim,DMCompositeSampleGLVisFields_Private,(PetscObject*)Ufds,ctx,DestroyGLVisViewerCtx_Private);CHKERRQ(ierr);
1446 for (i = 0; i < tnf; i++) {
1447 ierr = PetscFree(fecs[i]);CHKERRQ(ierr);
1448 }
1449 ierr = PetscFree3(fecs,sdim,Ufds);CHKERRQ(ierr);
1450 PetscFunctionReturn(0);
1451 }
1452
DMRefine_Composite(DM dmi,MPI_Comm comm,DM * fine)1453 PetscErrorCode DMRefine_Composite(DM dmi,MPI_Comm comm,DM *fine)
1454 {
1455 PetscErrorCode ierr;
1456 struct DMCompositeLink *next;
1457 DM_Composite *com = (DM_Composite*)dmi->data;
1458 DM dm;
1459
1460 PetscFunctionBegin;
1461 PetscValidHeaderSpecific(dmi,DM_CLASSID,1);
1462 if (comm == MPI_COMM_NULL) {
1463 ierr = PetscObjectGetComm((PetscObject)dmi,&comm);CHKERRQ(ierr);
1464 }
1465 ierr = DMSetUp(dmi);CHKERRQ(ierr);
1466 next = com->next;
1467 ierr = DMCompositeCreate(comm,fine);CHKERRQ(ierr);
1468
1469 /* loop over packed objects, handling one at at time */
1470 while (next) {
1471 ierr = DMRefine(next->dm,comm,&dm);CHKERRQ(ierr);
1472 ierr = DMCompositeAddDM(*fine,dm);CHKERRQ(ierr);
1473 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr);
1474 next = next->next;
1475 }
1476 PetscFunctionReturn(0);
1477 }
1478
DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM * fine)1479 PetscErrorCode DMCoarsen_Composite(DM dmi,MPI_Comm comm,DM *fine)
1480 {
1481 PetscErrorCode ierr;
1482 struct DMCompositeLink *next;
1483 DM_Composite *com = (DM_Composite*)dmi->data;
1484 DM dm;
1485
1486 PetscFunctionBegin;
1487 PetscValidHeaderSpecific(dmi,DM_CLASSID,1);
1488 ierr = DMSetUp(dmi);CHKERRQ(ierr);
1489 if (comm == MPI_COMM_NULL) {
1490 ierr = PetscObjectGetComm((PetscObject)dmi,&comm);CHKERRQ(ierr);
1491 }
1492 next = com->next;
1493 ierr = DMCompositeCreate(comm,fine);CHKERRQ(ierr);
1494
1495 /* loop over packed objects, handling one at at time */
1496 while (next) {
1497 ierr = DMCoarsen(next->dm,comm,&dm);CHKERRQ(ierr);
1498 ierr = DMCompositeAddDM(*fine,dm);CHKERRQ(ierr);
1499 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr);
1500 next = next->next;
1501 }
1502 PetscFunctionReturn(0);
1503 }
1504
DMCreateInterpolation_Composite(DM coarse,DM fine,Mat * A,Vec * v)1505 PetscErrorCode DMCreateInterpolation_Composite(DM coarse,DM fine,Mat *A,Vec *v)
1506 {
1507 PetscErrorCode ierr;
1508 PetscInt m,n,M,N,nDM,i;
1509 struct DMCompositeLink *nextc;
1510 struct DMCompositeLink *nextf;
1511 Vec gcoarse,gfine,*vecs;
1512 DM_Composite *comcoarse = (DM_Composite*)coarse->data;
1513 DM_Composite *comfine = (DM_Composite*)fine->data;
1514 Mat *mats;
1515
1516 PetscFunctionBegin;
1517 PetscValidHeaderSpecific(coarse,DM_CLASSID,1);
1518 PetscValidHeaderSpecific(fine,DM_CLASSID,2);
1519 ierr = DMSetUp(coarse);CHKERRQ(ierr);
1520 ierr = DMSetUp(fine);CHKERRQ(ierr);
1521 /* use global vectors only for determining matrix layout */
1522 ierr = DMGetGlobalVector(coarse,&gcoarse);CHKERRQ(ierr);
1523 ierr = DMGetGlobalVector(fine,&gfine);CHKERRQ(ierr);
1524 ierr = VecGetLocalSize(gcoarse,&n);CHKERRQ(ierr);
1525 ierr = VecGetLocalSize(gfine,&m);CHKERRQ(ierr);
1526 ierr = VecGetSize(gcoarse,&N);CHKERRQ(ierr);
1527 ierr = VecGetSize(gfine,&M);CHKERRQ(ierr);
1528 ierr = DMRestoreGlobalVector(coarse,&gcoarse);CHKERRQ(ierr);
1529 ierr = DMRestoreGlobalVector(fine,&gfine);CHKERRQ(ierr);
1530
1531 nDM = comfine->nDM;
1532 if (nDM != comcoarse->nDM) SETERRQ2(PetscObjectComm((PetscObject)fine),PETSC_ERR_ARG_INCOMP,"Fine DMComposite has %D entries, but coarse has %D",nDM,comcoarse->nDM);
1533 ierr = PetscCalloc1(nDM*nDM,&mats);CHKERRQ(ierr);
1534 if (v) {
1535 ierr = PetscCalloc1(nDM,&vecs);CHKERRQ(ierr);
1536 }
1537
1538 /* loop over packed objects, handling one at at time */
1539 for (nextc=comcoarse->next,nextf=comfine->next,i=0; nextc; nextc=nextc->next,nextf=nextf->next,i++) {
1540 if (!v) {
1541 ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],NULL);CHKERRQ(ierr);
1542 } else {
1543 ierr = DMCreateInterpolation(nextc->dm,nextf->dm,&mats[i*nDM+i],&vecs[i]);CHKERRQ(ierr);
1544 }
1545 }
1546 ierr = MatCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,nDM,NULL,mats,A);CHKERRQ(ierr);
1547 if (v) {
1548 ierr = VecCreateNest(PetscObjectComm((PetscObject)fine),nDM,NULL,vecs,v);CHKERRQ(ierr);
1549 }
1550 for (i=0; i<nDM*nDM; i++) {ierr = MatDestroy(&mats[i]);CHKERRQ(ierr);}
1551 ierr = PetscFree(mats);CHKERRQ(ierr);
1552 if (v) {
1553 for (i=0; i<nDM; i++) {ierr = VecDestroy(&vecs[i]);CHKERRQ(ierr);}
1554 ierr = PetscFree(vecs);CHKERRQ(ierr);
1555 }
1556 PetscFunctionReturn(0);
1557 }
1558
DMGetLocalToGlobalMapping_Composite(DM dm)1559 static PetscErrorCode DMGetLocalToGlobalMapping_Composite(DM dm)
1560 {
1561 DM_Composite *com = (DM_Composite*)dm->data;
1562 ISLocalToGlobalMapping *ltogs;
1563 PetscInt i;
1564 PetscErrorCode ierr;
1565
1566 PetscFunctionBegin;
1567 /* Set the ISLocalToGlobalMapping on the new matrix */
1568 ierr = DMCompositeGetISLocalToGlobalMappings(dm,<ogs);CHKERRQ(ierr);
1569 ierr = ISLocalToGlobalMappingConcatenate(PetscObjectComm((PetscObject)dm),com->nDM,ltogs,&dm->ltogmap);CHKERRQ(ierr);
1570 for (i=0; i<com->nDM; i++) {ierr = ISLocalToGlobalMappingDestroy(<ogs[i]);CHKERRQ(ierr);}
1571 ierr = PetscFree(ltogs);CHKERRQ(ierr);
1572 PetscFunctionReturn(0);
1573 }
1574
1575
DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring * coloring)1576 PetscErrorCode DMCreateColoring_Composite(DM dm,ISColoringType ctype,ISColoring *coloring)
1577 {
1578 PetscErrorCode ierr;
1579 PetscInt n,i,cnt;
1580 ISColoringValue *colors;
1581 PetscBool dense = PETSC_FALSE;
1582 ISColoringValue maxcol = 0;
1583 DM_Composite *com = (DM_Composite*)dm->data;
1584
1585 PetscFunctionBegin;
1586 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1587 if (ctype == IS_COLORING_LOCAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only global coloring supported");
1588 else if (ctype == IS_COLORING_GLOBAL) {
1589 n = com->n;
1590 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unknown ISColoringType");
1591 ierr = PetscMalloc1(n,&colors);CHKERRQ(ierr); /* freed in ISColoringDestroy() */
1592
1593 ierr = PetscOptionsGetBool(((PetscObject)dm)->options,((PetscObject)dm)->prefix,"-dmcomposite_dense_jacobian",&dense,NULL);CHKERRQ(ierr);
1594 if (dense) {
1595 for (i=0; i<n; i++) {
1596 colors[i] = (ISColoringValue)(com->rstart + i);
1597 }
1598 maxcol = com->N;
1599 } else {
1600 struct DMCompositeLink *next = com->next;
1601 PetscMPIInt rank;
1602
1603 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
1604 cnt = 0;
1605 while (next) {
1606 ISColoring lcoloring;
1607
1608 ierr = DMCreateColoring(next->dm,IS_COLORING_GLOBAL,&lcoloring);CHKERRQ(ierr);
1609 for (i=0; i<lcoloring->N; i++) {
1610 colors[cnt++] = maxcol + lcoloring->colors[i];
1611 }
1612 maxcol += lcoloring->n;
1613 ierr = ISColoringDestroy(&lcoloring);CHKERRQ(ierr);
1614 next = next->next;
1615 }
1616 }
1617 ierr = ISColoringCreate(PetscObjectComm((PetscObject)dm),maxcol,n,colors,PETSC_OWN_POINTER,coloring);CHKERRQ(ierr);
1618 PetscFunctionReturn(0);
1619 }
1620
DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)1621 PetscErrorCode DMGlobalToLocalBegin_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1622 {
1623 PetscErrorCode ierr;
1624 struct DMCompositeLink *next;
1625 PetscScalar *garray,*larray;
1626 DM_Composite *com = (DM_Composite*)dm->data;
1627
1628 PetscFunctionBegin;
1629 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1630 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
1631
1632 if (!com->setup) {
1633 ierr = DMSetUp(dm);CHKERRQ(ierr);
1634 }
1635
1636 ierr = VecGetArray(gvec,&garray);CHKERRQ(ierr);
1637 ierr = VecGetArray(lvec,&larray);CHKERRQ(ierr);
1638
1639 /* loop over packed objects, handling one at at time */
1640 next = com->next;
1641 while (next) {
1642 Vec local,global;
1643 PetscInt N;
1644
1645 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
1646 ierr = VecGetLocalSize(global,&N);CHKERRQ(ierr);
1647 ierr = VecPlaceArray(global,garray);CHKERRQ(ierr);
1648 ierr = DMGetLocalVector(next->dm,&local);CHKERRQ(ierr);
1649 ierr = VecPlaceArray(local,larray);CHKERRQ(ierr);
1650 ierr = DMGlobalToLocalBegin(next->dm,global,mode,local);CHKERRQ(ierr);
1651 ierr = DMGlobalToLocalEnd(next->dm,global,mode,local);CHKERRQ(ierr);
1652 ierr = VecResetArray(global);CHKERRQ(ierr);
1653 ierr = VecResetArray(local);CHKERRQ(ierr);
1654 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
1655 ierr = DMRestoreLocalVector(next->dm,&local);CHKERRQ(ierr);
1656
1657 larray += next->nlocal;
1658 garray += next->n;
1659 next = next->next;
1660 }
1661
1662 ierr = VecRestoreArray(gvec,NULL);CHKERRQ(ierr);
1663 ierr = VecRestoreArray(lvec,NULL);CHKERRQ(ierr);
1664 PetscFunctionReturn(0);
1665 }
1666
DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)1667 PetscErrorCode DMGlobalToLocalEnd_Composite(DM dm,Vec gvec,InsertMode mode,Vec lvec)
1668 {
1669 PetscFunctionBegin;
1670 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1671 PetscValidHeaderSpecific(gvec,VEC_CLASSID,2);
1672 PetscValidHeaderSpecific(lvec,VEC_CLASSID,4);
1673 PetscFunctionReturn(0);
1674 }
1675
DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)1676 PetscErrorCode DMLocalToGlobalBegin_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1677 {
1678 PetscErrorCode ierr;
1679 struct DMCompositeLink *next;
1680 PetscScalar *larray,*garray;
1681 DM_Composite *com = (DM_Composite*)dm->data;
1682
1683 PetscFunctionBegin;
1684 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1685 PetscValidHeaderSpecific(lvec,VEC_CLASSID,2);
1686 PetscValidHeaderSpecific(gvec,VEC_CLASSID,4);
1687
1688 if (!com->setup) {
1689 ierr = DMSetUp(dm);CHKERRQ(ierr);
1690 }
1691
1692 ierr = VecGetArray(lvec,&larray);CHKERRQ(ierr);
1693 ierr = VecGetArray(gvec,&garray);CHKERRQ(ierr);
1694
1695 /* loop over packed objects, handling one at at time */
1696 next = com->next;
1697 while (next) {
1698 Vec global,local;
1699
1700 ierr = DMGetLocalVector(next->dm,&local);CHKERRQ(ierr);
1701 ierr = VecPlaceArray(local,larray);CHKERRQ(ierr);
1702 ierr = DMGetGlobalVector(next->dm,&global);CHKERRQ(ierr);
1703 ierr = VecPlaceArray(global,garray);CHKERRQ(ierr);
1704 ierr = DMLocalToGlobalBegin(next->dm,local,mode,global);CHKERRQ(ierr);
1705 ierr = DMLocalToGlobalEnd(next->dm,local,mode,global);CHKERRQ(ierr);
1706 ierr = VecResetArray(local);CHKERRQ(ierr);
1707 ierr = VecResetArray(global);CHKERRQ(ierr);
1708 ierr = DMRestoreGlobalVector(next->dm,&global);CHKERRQ(ierr);
1709 ierr = DMRestoreLocalVector(next->dm,&local);CHKERRQ(ierr);
1710
1711 garray += next->n;
1712 larray += next->nlocal;
1713 next = next->next;
1714 }
1715
1716 ierr = VecRestoreArray(gvec,NULL);CHKERRQ(ierr);
1717 ierr = VecRestoreArray(lvec,NULL);CHKERRQ(ierr);
1718
1719 PetscFunctionReturn(0);
1720 }
1721
DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)1722 PetscErrorCode DMLocalToGlobalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1723 {
1724 PetscFunctionBegin;
1725 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1726 PetscValidHeaderSpecific(lvec,VEC_CLASSID,2);
1727 PetscValidHeaderSpecific(gvec,VEC_CLASSID,4);
1728 PetscFunctionReturn(0);
1729 }
1730
DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)1731 PetscErrorCode DMLocalToLocalBegin_Composite(DM dm,Vec vec1,InsertMode mode,Vec vec2)
1732 {
1733 PetscErrorCode ierr;
1734 struct DMCompositeLink *next;
1735 PetscScalar *array1,*array2;
1736 DM_Composite *com = (DM_Composite*)dm->data;
1737
1738 PetscFunctionBegin;
1739 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1740 PetscValidHeaderSpecific(vec1,VEC_CLASSID,2);
1741 PetscValidHeaderSpecific(vec2,VEC_CLASSID,4);
1742
1743 if (!com->setup) {
1744 ierr = DMSetUp(dm);CHKERRQ(ierr);
1745 }
1746
1747 ierr = VecGetArray(vec1,&array1);CHKERRQ(ierr);
1748 ierr = VecGetArray(vec2,&array2);CHKERRQ(ierr);
1749
1750 /* loop over packed objects, handling one at at time */
1751 next = com->next;
1752 while (next) {
1753 Vec local1,local2;
1754
1755 ierr = DMGetLocalVector(next->dm,&local1);CHKERRQ(ierr);
1756 ierr = VecPlaceArray(local1,array1);CHKERRQ(ierr);
1757 ierr = DMGetLocalVector(next->dm,&local2);CHKERRQ(ierr);
1758 ierr = VecPlaceArray(local2,array2);CHKERRQ(ierr);
1759 ierr = DMLocalToLocalBegin(next->dm,local1,mode,local2);CHKERRQ(ierr);
1760 ierr = DMLocalToLocalEnd(next->dm,local1,mode,local2);CHKERRQ(ierr);
1761 ierr = VecResetArray(local2);CHKERRQ(ierr);
1762 ierr = DMRestoreLocalVector(next->dm,&local2);CHKERRQ(ierr);
1763 ierr = VecResetArray(local1);CHKERRQ(ierr);
1764 ierr = DMRestoreLocalVector(next->dm,&local1);CHKERRQ(ierr);
1765
1766 array1 += next->nlocal;
1767 array2 += next->nlocal;
1768 next = next->next;
1769 }
1770
1771 ierr = VecRestoreArray(vec1,NULL);CHKERRQ(ierr);
1772 ierr = VecRestoreArray(vec2,NULL);CHKERRQ(ierr);
1773
1774 PetscFunctionReturn(0);
1775 }
1776
DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)1777 PetscErrorCode DMLocalToLocalEnd_Composite(DM dm,Vec lvec,InsertMode mode,Vec gvec)
1778 {
1779 PetscFunctionBegin;
1780 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1781 PetscValidHeaderSpecific(lvec,VEC_CLASSID,2);
1782 PetscValidHeaderSpecific(gvec,VEC_CLASSID,4);
1783 PetscFunctionReturn(0);
1784 }
1785
1786 /*MC
1787 DMCOMPOSITE = "composite" - A DM object that is used to manage data for a collection of DMs
1788
1789 Level: intermediate
1790
1791 .seealso: DMType, DM, DMDACreate(), DMCreate(), DMSetType(), DMCompositeCreate()
1792 M*/
1793
1794
DMCreate_Composite(DM p)1795 PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM p)
1796 {
1797 PetscErrorCode ierr;
1798 DM_Composite *com;
1799
1800 PetscFunctionBegin;
1801 ierr = PetscNewLog(p,&com);CHKERRQ(ierr);
1802 p->data = com;
1803 com->n = 0;
1804 com->nghost = 0;
1805 com->next = NULL;
1806 com->nDM = 0;
1807
1808 p->ops->createglobalvector = DMCreateGlobalVector_Composite;
1809 p->ops->createlocalvector = DMCreateLocalVector_Composite;
1810 p->ops->getlocaltoglobalmapping = DMGetLocalToGlobalMapping_Composite;
1811 p->ops->createfieldis = DMCreateFieldIS_Composite;
1812 p->ops->createfielddecomposition = DMCreateFieldDecomposition_Composite;
1813 p->ops->refine = DMRefine_Composite;
1814 p->ops->coarsen = DMCoarsen_Composite;
1815 p->ops->createinterpolation = DMCreateInterpolation_Composite;
1816 p->ops->creatematrix = DMCreateMatrix_Composite;
1817 p->ops->getcoloring = DMCreateColoring_Composite;
1818 p->ops->globaltolocalbegin = DMGlobalToLocalBegin_Composite;
1819 p->ops->globaltolocalend = DMGlobalToLocalEnd_Composite;
1820 p->ops->localtoglobalbegin = DMLocalToGlobalBegin_Composite;
1821 p->ops->localtoglobalend = DMLocalToGlobalEnd_Composite;
1822 p->ops->localtolocalbegin = DMLocalToLocalBegin_Composite;
1823 p->ops->localtolocalend = DMLocalToLocalEnd_Composite;
1824 p->ops->destroy = DMDestroy_Composite;
1825 p->ops->view = DMView_Composite;
1826 p->ops->setup = DMSetUp_Composite;
1827
1828 ierr = PetscObjectComposeFunction((PetscObject)p,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Composite);CHKERRQ(ierr);
1829 PetscFunctionReturn(0);
1830 }
1831
1832 /*@
1833 DMCompositeCreate - Creates a vector packer, used to generate "composite"
1834 vectors made up of several subvectors.
1835
1836 Collective
1837
1838 Input Parameter:
1839 . comm - the processors that will share the global vector
1840
1841 Output Parameters:
1842 . packer - the packer object
1843
1844 Level: advanced
1845
1846 .seealso DMDestroy(), DMCompositeAddDM(), DMCompositeScatter(), DMCOMPOSITE,DMCreate()
1847 DMCompositeGather(), DMCreateGlobalVector(), DMCompositeGetISLocalToGlobalMappings(), DMCompositeGetAccess()
1848 DMCompositeGetLocalVectors(), DMCompositeRestoreLocalVectors(), DMCompositeGetEntries()
1849
1850 @*/
DMCompositeCreate(MPI_Comm comm,DM * packer)1851 PetscErrorCode DMCompositeCreate(MPI_Comm comm,DM *packer)
1852 {
1853 PetscErrorCode ierr;
1854
1855 PetscFunctionBegin;
1856 PetscValidPointer(packer,2);
1857 ierr = DMCreate(comm,packer);CHKERRQ(ierr);
1858 ierr = DMSetType(*packer,DMCOMPOSITE);CHKERRQ(ierr);
1859 PetscFunctionReturn(0);
1860 }
1861