Actual source code: tagger.c
1: #include <petsc/private/vecimpl.h>
3: /*@C
4: VecTaggerCreate - create a Vec tagger context. This object is used to control the tagging/selection of index sets
5: based on the values in a vector. This is used, for example, in adaptive simulations when aspects are selected for
6: refinement or coarsening. The primary intent is that the selected index sets are based purely on the values in the
7: vector, though implementations that do not follow this intent are possible.
9: Once a VecTagger is created (VecTaggerCreate()), optionally modified by options (VecTaggerSetFromOptions()), and
10: set up (VecTaggerSetUp()), it is applied to vectors with VecTaggerComputeIS() to comute the selected index sets.
12: In many cases, the selection criteria for an index is whether the corresponding value falls within a collection of
13: boxes: for this common case, VecTaggerCreateBoxes() can also be used to determine those boxes.
15: Provided implementations support tagging based on a box/interval of values (VECTAGGERABSOLUTE), based on a box of
16: values of relative to the range of values present in the vector (VECTAGGERRELATIVE), based on where values fall in
17: the cumulative distribution of values in the vector (VECTAGGERCDF), and based on unions (VECTAGGEROR) or
18: intersections (VECTAGGERAND) of other criteria.
20: Collective
22: Input Parameter:
23: . comm - communicator on which the vec tagger will operate
25: Output Parameter:
26: . tagger - new Vec tagger context
28: Level: advanced
30: .seealso: VecTaggerSetBlockSize(), VecTaggerSetFromOptions(), VecTaggerSetUp(), VecTaggerComputeIS(), VecTaggerComputeBoxes(), VecTaggerDestroy()
31: @*/
32: PetscErrorCode VecTaggerCreate(MPI_Comm comm,VecTagger *tagger)
33: {
34: VecTagger b;
37: VecTaggerInitializePackage();
39: PetscHeaderCreate(b,VEC_TAGGER_CLASSID,"VecTagger","Vec Tagger","Vec",comm,VecTaggerDestroy,VecTaggerView);
41: b->blocksize = 1;
42: b->invert = PETSC_FALSE;
43: b->setupcalled = PETSC_FALSE;
45: *tagger = b;
46: return 0;
47: }
49: /*@C
50: VecTaggerSetType - set the Vec tagger implementation
52: Collective on VecTagger
54: Input Parameters:
55: + tagger - the VecTagger context
56: - type - a known method
58: Options Database Key:
59: . -vec_tagger_type <type> - Sets the method; use -help for a list
60: of available methods (for instance, absolute, relative, cdf, or, and)
62: Notes:
63: See "include/petscvec.h" for available methods (for instance)
64: + VECTAGGERABSOLUTE - tag based on a box of values
65: . VECTAGGERRELATIVE - tag based on a box relative to the range of values present in the vector
66: . VECTAGGERCDF - tag based on a box in the cumulative distribution of values present in the vector
67: . VECTAGGEROR - tag based on the union of a set of VecTagger contexts
68: - VECTAGGERAND - tag based on the intersection of a set of other VecTagger contexts
70: Level: advanced
72: .seealso: VecTaggerType, VecTaggerCreate(), VecTagger
73: @*/
74: PetscErrorCode VecTaggerSetType(VecTagger tagger,VecTaggerType type)
75: {
76: PetscBool match;
77: PetscErrorCode (*r)(VecTagger);
82: PetscObjectTypeCompare((PetscObject)tagger,type,&match);
83: if (match) return 0;
85: PetscFunctionListFind(VecTaggerList,type,&r);
87: /* Destroy the previous private VecTagger context */
88: if (tagger->ops->destroy) (*(tagger)->ops->destroy)(tagger);
89: PetscMemzero(tagger->ops,sizeof(*tagger->ops));
90: PetscObjectChangeTypeName((PetscObject)tagger,type);
91: tagger->ops->create = r;
92: (*r)(tagger);
93: return 0;
94: }
96: /*@C
97: VecTaggerGetType - Gets the VecTagger type name (as a string) from the VecTagger.
99: Not Collective
101: Input Parameter:
102: . tagger - The Vec tagger context
104: Output Parameter:
105: . type - The VecTagger type name
107: Level: advanced
109: .seealso: VecTaggerSetType(), VecTaggerCreate(), VecTaggerSetFromOptions(), VecTagger, VecTaggerType
110: @*/
111: PetscErrorCode VecTaggerGetType(VecTagger tagger, VecTaggerType *type)
112: {
115: VecTaggerRegisterAll();
116: *type = ((PetscObject)tagger)->type_name;
117: return 0;
118: }
120: /*@
121: VecTaggerDestroy - destroy a VecTagger context
123: Collective
125: Input Parameter:
126: . tagger - address of tagger
128: Level: advanced
130: .seealso: VecTaggerCreate(), VecTaggerSetType(), VecTagger
131: @*/
132: PetscErrorCode VecTaggerDestroy(VecTagger *tagger)
133: {
134: if (!*tagger) return 0;
136: if (--((PetscObject)(*tagger))->refct > 0) {*tagger = NULL; return 0;}
137: if ((*tagger)->ops->destroy) (*(*tagger)->ops->destroy)(*tagger);
138: PetscHeaderDestroy(tagger);
139: return 0;
140: }
142: /*@
143: VecTaggerSetUp - set up a VecTagger context
145: Collective
147: Input Parameter:
148: . tagger - Vec tagger object
150: Level: advanced
152: .seealso: VecTaggerSetFromOptions(), VecTaggerSetType(), VecTagger, VecTaggerCreate(), VecTaggerSetUp()
153: @*/
154: PetscErrorCode VecTaggerSetUp(VecTagger tagger)
155: {
156: if (tagger->setupcalled) return 0;
157: if (!((PetscObject)tagger)->type_name) VecTaggerSetType(tagger,VECTAGGERABSOLUTE);
158: if (tagger->ops->setup) (*tagger->ops->setup)(tagger);
159: tagger->setupcalled = PETSC_TRUE;
160: return 0;
161: }
163: /*@C
164: VecTaggerSetFromOptions - set VecTagger options using the options database
166: Logically Collective
168: Input Parameter:
169: . tagger - vec tagger
171: Options Database Keys:
172: + -vec_tagger_type - implementation type, see VecTaggerSetType()
173: . -vec_tagger_block_size - set the block size, see VecTaggerSetBlockSize()
174: - -vec_tagger_invert - invert the index set returned by VecTaggerComputeIS()
176: Level: advanced
178: .seealso: VecTagger, VecTaggerCreate(), VecTaggerSetUp()
180: @*/
181: PetscErrorCode VecTaggerSetFromOptions(VecTagger tagger)
182: {
183: VecTaggerType deft;
184: char type[256];
186: PetscBool flg;
189: PetscObjectOptionsBegin((PetscObject)tagger);
190: deft = ((PetscObject)tagger)->type_name ? ((PetscObject)tagger)->type_name : VECTAGGERABSOLUTE;
191: PetscOptionsFList("-vec_tagger_type","VecTagger implementation type","VecTaggerSetType",VecTaggerList,deft,type,256,&flg);
192: VecTaggerSetType(tagger,flg ? type : deft);
193: PetscOptionsInt("-vec_tagger_block_size","block size of the vectors the tagger operates on","VecTaggerSetBlockSize",tagger->blocksize,&tagger->blocksize,NULL);
194: PetscOptionsBool("-vec_tagger_invert","invert the set of indices returned by VecTaggerComputeIS()","VecTaggerSetInvert",tagger->invert,&tagger->invert,NULL);
195: if (tagger->ops->setfromoptions) (*tagger->ops->setfromoptions)(PetscOptionsObject,tagger);
196: PetscOptionsEnd();
197: return 0;
198: }
200: /*@C
201: VecTaggerSetBlockSize - block size of the set of indices returned by VecTaggerComputeIS(). Values greater than one
202: are useful when there are multiple criteria for determining which indices to include in the set. For example,
203: consider adaptive mesh refinement in a multiphysics problem, with metrics of solution quality for multiple fields
204: measure on each cell. The size of the vector will be [numCells * numFields]; the VecTagger block size should be
205: numFields; VecTaggerComputeIS() will return indices in the range [0,numCells), i.e., one index is given for each
206: block of values.
208: Note that the block size of the vector does not have to match.
210: Note also that the index set created in VecTaggerComputeIS() has block size: it is an index set over the list of
211: items that the vector refers to, not to the vector itself.
213: Logically Collective
215: Input Parameters:
216: + tagger - vec tagger
217: - blocksize - block size of the criteria used to tagger vectors
219: Level: advanced
221: .seealso: VecTaggerComputeIS(), VecTaggerGetBlockSize(), VecSetBlockSize(), VecGetBlockSize(), VecTagger, VecTaggerCreate()
222: @*/
223: PetscErrorCode VecTaggerSetBlockSize(VecTagger tagger, PetscInt blocksize)
224: {
227: tagger->blocksize = blocksize;
228: return 0;
229: }
231: /*@C
232: VecTaggerGetBlockSize - get the block size of the indices created by VecTaggerComputeIS().
234: Logically Collective
236: Input Parameter:
237: . tagger - vec tagger
239: Output Parameter:
240: . blocksize - block size of the vectors the tagger operates on
242: Level: advanced
244: .seealso: VecTaggerComputeIS(), VecTaggerSetBlockSize(), VecTagger, VecTaggerCreate()
245: @*/
246: PetscErrorCode VecTaggerGetBlockSize(VecTagger tagger, PetscInt *blocksize)
247: {
250: *blocksize = tagger->blocksize;
251: return 0;
252: }
254: /*@C
255: VecTaggerSetInvert - If the tagged index sets are based on boxes that can be returned by VecTaggerComputeBoxes(),
256: then this option inverts values used to compute the IS, i.e., from being in the union of the boxes to being in the
257: intersection of their exteriors.
259: Logically Collective
261: Input Parameters:
262: + tagger - vec tagger
263: - invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
265: Level: advanced
267: .seealso: VecTaggerComputeIS(), VecTaggerGetInvert(), VecTagger, VecTaggerCreate()
268: @*/
269: PetscErrorCode VecTaggerSetInvert(VecTagger tagger, PetscBool invert)
270: {
273: tagger->invert = invert;
274: return 0;
275: }
277: /*@C
278: VecTaggerGetInvert - get whether the set of indices returned by VecTaggerComputeIS() are inverted
280: Logically Collective
282: Input Parameter:
283: . tagger - vec tagger
285: Output Parameter:
286: . invert - PETSC_TRUE to invert, PETSC_FALSE to use the indices as is
288: Level: advanced
290: .seealso: VecTaggerComputeIS(), VecTaggerSetInvert(), VecTagger, VecTaggerCreate()
291: @*/
292: PetscErrorCode VecTaggerGetInvert(VecTagger tagger, PetscBool *invert)
293: {
296: *invert = tagger->invert;
297: return 0;
298: }
300: /*@C
301: VecTaggerView - view a VecTagger context
303: Collective
305: Input Parameters:
306: + tagger - vec tagger
307: - viewer - viewer to display tagger, for example PETSC_VIEWER_STDOUT_WORLD
309: Level: advanced
311: .seealso: VecTaggerCreate(), VecTagger, VecTaggerCreate()
312: @*/
313: PetscErrorCode VecTaggerView(VecTagger tagger,PetscViewer viewer)
314: {
315: PetscBool iascii;
318: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)tagger),&viewer);
321: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
322: if (iascii) {
323: PetscObjectPrintClassNamePrefixType((PetscObject)tagger,viewer);
324: PetscViewerASCIIPushTab(viewer);
325: PetscViewerASCIIPrintf(viewer,"Block size: %" PetscInt_FMT "\n",tagger->blocksize);
326: if (tagger->ops->view) (*tagger->ops->view)(tagger,viewer);
327: if (tagger->invert) PetscViewerASCIIPrintf(viewer,"Inverting ISs.\n");
328: PetscViewerASCIIPopTab(viewer);
329: }
330: return 0;
331: }
333: /*@C
334: VecTaggerComputeBoxes - If the tagged index set can be summarized as a list of boxes of values, returns that list, otherwise returns
335: in listed PETSC_FALSE
337: Collective on VecTagger
339: Input Parameters:
340: + tagger - the VecTagger context
341: - vec - the vec to tag
343: Output Parameters:
344: + numBoxes - the number of boxes in the tag definition
345: . boxes - a newly allocated list of boxes. This is a flat array of (BlockSize * numBoxes) pairs that the user can free with PetscFree().
346: - listed - PETSC_TRUE if a list was created, pass in NULL if not needed
348: Notes:
349: A value is tagged if it is in any of the boxes, unless the tagger has been inverted (see VecTaggerSetInvert()/VecTaggerGetInvert()), in which case a value is tagged if it is in none of the boxes.
351: Level: advanced
353: .seealso: VecTaggerComputeIS(), VecTagger, VecTaggerCreate()
354: @*/
355: PetscErrorCode VecTaggerComputeBoxes(VecTagger tagger,Vec vec,PetscInt *numBoxes,VecTaggerBox **boxes,PetscBool *listed)
356: {
357: PetscInt vls, tbs;
363: VecGetLocalSize(vec,&vls);
364: VecTaggerGetBlockSize(tagger,&tbs);
366: if (tagger->ops->computeboxes) {
367: *listed = PETSC_TRUE;
368: (*tagger->ops->computeboxes) (tagger,vec,numBoxes,boxes,listed);
369: } else *listed = PETSC_FALSE;
370: return 0;
371: }
373: /*@C
374: VecTaggerComputeIS - Use a VecTagger context to tag a set of indices based on a vector's values
376: Collective on VecTagger
378: Input Parameters:
379: + tagger - the VecTagger context
380: - vec - the vec to tag
382: Output Parameters:
383: + IS - a list of the indices tagged by the tagger, i.e., if the number of local indices will be n / bs, where n is VecGetLocalSize() and bs is VecTaggerGetBlockSize().
384: - listed - routine was able to compute the IS, pass in NULL if not needed
386: Level: advanced
388: .seealso: VecTaggerComputeBoxes(), VecTagger, VecTaggerCreate()
389: @*/
390: PetscErrorCode VecTaggerComputeIS(VecTagger tagger,Vec vec,IS *is,PetscBool *listed)
391: {
392: PetscInt vls, tbs;
397: VecGetLocalSize(vec,&vls);
398: VecTaggerGetBlockSize(tagger,&tbs);
400: if (tagger->ops->computeis) {
401: (*tagger->ops->computeis) (tagger,vec,is,listed);
402: } else if (listed) *listed = PETSC_FALSE;
403: return 0;
404: }
406: PetscErrorCode VecTaggerComputeIS_FromBoxes(VecTagger tagger, Vec vec, IS *is,PetscBool *listed)
407: {
408: PetscInt numBoxes;
409: VecTaggerBox *boxes;
410: PetscInt numTagged, offset;
411: PetscInt *tagged;
412: PetscInt bs, b, i, j, k, n;
413: PetscBool invert;
414: const PetscScalar *vecArray;
415: PetscBool boxlisted;
417: VecTaggerGetBlockSize(tagger,&bs);
418: VecTaggerComputeBoxes(tagger,vec,&numBoxes,&boxes,&boxlisted);
419: if (!boxlisted) {
420: if (listed) *listed = PETSC_FALSE;
421: return 0;
422: }
423: VecGetArrayRead(vec, &vecArray);
424: VecGetLocalSize(vec, &n);
425: invert = tagger->invert;
426: numTagged = 0;
427: offset = 0;
428: tagged = NULL;
430: n /= bs;
431: for (i = 0; i < 2; i++) {
432: if (i) {
433: PetscMalloc1(numTagged,&tagged);
434: }
435: for (j = 0; j < n; j++) {
437: for (k = 0; k < numBoxes; k++) {
438: for (b = 0; b < bs; b++) {
439: PetscScalar val = vecArray[j * bs + b];
440: PetscInt l = k * bs + b;
441: VecTaggerBox box;
442: PetscBool in;
444: box = boxes[l];
445: #if !defined(PETSC_USE_COMPLEX)
446: in = (PetscBool) ((box.min <= val) && (val <= box.max));
447: #else
448: in = (PetscBool) ((PetscRealPart (box.min) <= PetscRealPart (val)) &&
449: (PetscImaginaryPart(box.min) <= PetscImaginaryPart(val)) &&
450: (PetscRealPart (val) <= PetscRealPart (box.max)) &&
451: (PetscImaginaryPart(val) <= PetscImaginaryPart(box.max)));
452: #endif
453: if (!in) break;
454: }
455: if (b == bs) break;
456: }
457: if ((PetscBool)(k < numBoxes) ^ invert) {
458: if (!i) numTagged++;
459: else tagged[offset++] = j;
460: }
461: }
462: }
463: VecRestoreArrayRead(vec, &vecArray);
464: PetscFree(boxes);
465: ISCreateGeneral(PetscObjectComm((PetscObject)vec),numTagged,tagged,PETSC_OWN_POINTER,is);
466: ISSort(*is);
467: if (listed) *listed = PETSC_TRUE;
468: return 0;
469: }