Actual source code: dmmoab.cxx
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/
2: #include <petsc-private/vecimpl.h> /*I "petscdm.h" I*/
4: #include <petscdmmoab.h>
5: #include <MBTagConventions.hpp>
6: #include <sstream>
8: typedef struct {
9: PetscInt bs; /* Number of degrees of freedom on each entity, aka tag size in moab */
10: PetscBool icreatedinstance; /* true if DM created moab instance internally, will destroy instance in DMDestroy */
11: moab::ParallelComm *pcomm;
12: moab::Interface *mbiface;
13: moab::Tag ltog_tag; /* moab supports "global id" tags, which are usually local to global numbering */
14: moab::Range range;
15: } DM_Moab;
17: typedef struct {
18: moab::Interface *mbiface;
19: moab::ParallelComm *pcomm;
20: moab::Range tag_range; /* entities to which this tag applies */
21: moab::Tag tag;
22: moab::Tag ltog_tag;
23: PetscInt tag_size;
24: PetscBool new_tag;
25: PetscBool serial;
27: } Vec_MOAB;
31: PetscErrorCode DMCreateGlobalVector_Moab(DM dm,Vec *gvec)
32: {
33: PetscErrorCode ierr;
34: DM_Moab *dmmoab = (DM_Moab*)dm->data;
39: PetscInt block_size = ((DM_Moab*)dm->data)->bs;
40: moab::Tag tag = 0;
41: DMMoabCreateVector(dm,tag,block_size,dmmoab->range,PETSC_FALSE,PETSC_TRUE,gvec);
42: return(0);
43: }
48: PetscErrorCode DMCreateLocalVector_Moab(DM dm,Vec *gvec)
49: {
50: PetscErrorCode ierr;
51: DM_Moab *dmmoab = (DM_Moab*)dm->data;
54: PetscInt bs = 1;
57: moab::Tag tag = 0;
58: DMMoabCreateVector(dm,tag,bs,dmmoab->range,PETSC_TRUE,PETSC_TRUE,gvec);
59: return(0);
60: }
64: PetscErrorCode DMDestroy_Moab(DM dm)
65: {
70: if (((DM_Moab*)dm->data)->icreatedinstance) {
71: delete ((DM_Moab*)dm->data)->mbiface;
72: ((DM_Moab*)dm->data)->mbiface = NULL;
73: ((DM_Moab*)dm->data)->pcomm = NULL;
74: ((DM_Moab*)dm->data)->range.~Range();
75: }
76: PetscFree(dm->data);
77: return(0);
78: }
82: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
83: {
84: DM_Moab *moab;
89: PetscNewLog(dm,DM_Moab,&moab);
90: dm->data = moab;
91: new (moab) DM_Moab();
93: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
94: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
95: dm->ops->destroy = DMDestroy_Moab;
96: return(0);
97: }
101: /*@
102: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
104: Collective on MPI_Comm
106: Input Parameter:
107: . comm - The communicator for the DMMoab object
109: Output Parameter:
110: . moab - The DMMoab object
112: Level: beginner
114: .keywords: DMMoab, create
115: @*/
116: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *moab)
117: {
122: DMCreate(comm, moab);
123: DMSetType(*moab, DMMOAB);
124: return(0);
125: }
129: /*@
130: DMMoabCreate - Creates a DMMoab object, optionally from an instance and other data
132: Collective on MPI_Comm
134: Input Parameter:
135: . comm - The communicator for the DMMoab object
136: . moab - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
137: along with the DMMoab
138: . pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
139: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
140: . range - If non-NULL, contains range of entities to which DOFs will be assigned
142: Output Parameter:
143: . moab - The DMMoab object
145: Level: beginner
147: .keywords: DMMoab, create
148: @*/
149: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::ParallelComm *pcomm, moab::Tag ltog_tag, moab::Range *range, DM *moab)
150: {
152: DM_Moab *dmmoab;
156: DMMoabCreate(comm, moab);
157: dmmoab = (DM_Moab*)(*moab)->data;
159: if (!mbiface) {
160: mbiface = new moab::Core();
161: dmmoab->icreatedinstance = PETSC_TRUE;
162: }
163: else
164: dmmoab->icreatedinstance = PETSC_FALSE;
166: if (!pcomm) {
167: PetscInt rank, nprocs;
168: MPI_Comm_rank(comm, &rank);
169: MPI_Comm_size(comm, &nprocs);
170: pcomm = new moab::ParallelComm(mbiface, comm);
171: }
173: // do the initialization of the DM
174: dmmoab->bs = 0;
175: dmmoab->pcomm = pcomm;
176: dmmoab->mbiface = mbiface;
177: dmmoab->ltog_tag = ltog_tag;
179: DMMoabSetInterface(*moab, mbiface);
180: if (!pcomm) pcomm = new moab::ParallelComm(mbiface, comm);
181: DMMoabSetParallelComm(*moab, pcomm);
182: if (!ltog_tag) {
183: moab::ErrorCode merr = mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, ltog_tag);MBERRNM(merr);
184: }
185: if (ltog_tag) {
186: DMMoabSetLocalToGlobalTag(*moab, ltog_tag);
187: }
188: if (range) {
189: DMMoabSetRange(*moab, *range);
190: }
191: return(0);
192: }
196: /*@
197: DMMoabSetParallelComm - Set the ParallelComm used with this DMMoab
199: Collective on MPI_Comm
201: Input Parameter:
202: . dm - The DMMoab object being set
203: . pcomm - The ParallelComm being set on the DMMoab
205: Level: beginner
207: .keywords: DMMoab, create
208: @*/
209: PetscErrorCode DMMoabSetParallelComm(DM dm,moab::ParallelComm *pcomm)
210: {
213: ((DM_Moab*)dm->data)->pcomm = pcomm;
214: ((DM_Moab*)dm->data)->mbiface = pcomm->get_moab();
215: return(0);
216: }
221: /*@
222: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
224: Collective on MPI_Comm
226: Input Parameter:
227: . dm - The DMMoab object being set
229: Output Parameter:
230: . pcomm - The ParallelComm for the DMMoab
232: Level: beginner
234: .keywords: DMMoab, create
235: @*/
236: PetscErrorCode DMMoabGetParallelComm(DM dm,moab::ParallelComm **pcomm)
237: {
240: *pcomm = ((DM_Moab*)dm->data)->pcomm;
241: return(0);
242: }
247: /*@
248: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
250: Collective on MPI_Comm
252: Input Parameter:
253: . dm - The DMMoab object being set
254: . mbiface - The MOAB instance being set on this DMMoab
256: Level: beginner
258: .keywords: DMMoab, create
259: @*/
260: PetscErrorCode DMMoabSetInterface(DM dm,moab::Interface *mbiface)
261: {
264: ((DM_Moab*)dm->data)->pcomm = NULL;
265: ((DM_Moab*)dm->data)->mbiface = mbiface;
266: return(0);
267: }
272: /*@
273: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
275: Collective on MPI_Comm
277: Input Parameter:
278: . dm - The DMMoab object being set
280: Output Parameter:
281: . mbiface - The MOAB instance set on this DMMoab
283: Level: beginner
285: .keywords: DMMoab, create
286: @*/
287: PetscErrorCode DMMoabGetInterface(DM dm,moab::Interface **mbiface)
288: {
291: *mbiface = ((DM_Moab*)dm->data)->mbiface;
292: return(0);
293: }
298: /*@
299: DMMoabSetRange - Set the entities having DOFs on this DMMoab
301: Collective on MPI_Comm
303: Input Parameter:
304: . dm - The DMMoab object being set
305: . range - The entities treated by this DMMoab
307: Level: beginner
309: .keywords: DMMoab, create
310: @*/
311: PetscErrorCode DMMoabSetRange(DM dm,moab::Range range)
312: {
315: ((DM_Moab*)dm->data)->range = range;
316: return(0);
317: }
322: /*@
323: DMMoabGetRange - Get the entities having DOFs on this DMMoab
325: Collective on MPI_Comm
327: Input Parameter:
328: . dm - The DMMoab object being set
330: Output Parameter:
331: . range - The entities treated by this DMMoab
333: Level: beginner
335: .keywords: DMMoab, create
336: @*/
337: PetscErrorCode DMMoabGetRange(DM dm,moab::Range *range)
338: {
341: *range = ((DM_Moab*)dm->data)->range;
342: return(0);
343: }
347: /*@
348: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
350: Collective on MPI_Comm
352: Input Parameter:
353: . dm - The DMMoab object being set
354: . ltogtag - The MOAB tag used for local to global ids
356: Level: beginner
358: .keywords: DMMoab, create
359: @*/
360: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm,moab::Tag ltogtag)
361: {
364: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
365: return(0);
366: }
371: /*@
372: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
374: Collective on MPI_Comm
376: Input Parameter:
377: . dm - The DMMoab object being set
379: Output Parameter:
380: . ltogtag - The MOAB tag used for local to global ids
382: Level: beginner
384: .keywords: DMMoab, create
385: @*/
386: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm,moab::Tag *ltog_tag)
387: {
390: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
391: return(0);
392: }
397: /*@
398: DMMoabSetBlockSize - Set the block size used with this DMMoab
400: Collective on MPI_Comm
402: Input Parameter:
403: . dm - The DMMoab object being set
404: . bs - The block size used with this DMMoab
406: Level: beginner
408: .keywords: DMMoab, create
409: @*/
410: PetscErrorCode DMMoabSetBlockSize(DM dm,PetscInt bs)
411: {
414: ((DM_Moab*)dm->data)->bs = bs;
415: return(0);
416: }
421: /*@
422: DMMoabGetBlockSize - Get the block size used with this DMMoab
424: Collective on MPI_Comm
426: Input Parameter:
427: . dm - The DMMoab object being set
429: Output Parameter:
430: . bs - The block size used with this DMMoab
432: Level: beginner
434: .keywords: DMMoab, create
435: @*/
436: PetscErrorCode DMMoabGetBlockSize(DM dm,PetscInt *bs)
437: {
440: *bs = ((DM_Moab*)dm->data)->bs;
441: return(0);
442: }
445: // declare for use later but before they're defined
446: PetscErrorCode DMMoab_VecUserDestroy(void *user);
447: PetscErrorCode DMMoab_VecDuplicate(Vec x,Vec *y);
448: PetscErrorCode DMMoab_CreateTagName(const moab::ParallelComm *pcomm,std::string& tag_name);
449: PetscErrorCode DMMoab_CreateVector(moab::Interface *iface,moab::ParallelComm *pcomm,moab::Tag tag,PetscInt tag_size,moab::Tag ltog_tag,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec);
453: /*@
454: DMMoabCreateVector - Create a Vec from either an existing tag, or a specified tag size, and a range of entities
456: Collective on MPI_Comm
458: Input Parameter:
459: . dm - The DMMoab object being set
460: . tag - If non-zero, block size will be taken from the tag size
461: . tag_size - If tag was zero, this parameter specifies the block size; unique tag name will be generated automatically
462: . range - If non-empty, Vec corresponds to these entities, otherwise to the entities set on the DMMoab
463: . serial - If true, this is a serial Vec, otherwise a parallel one
464: . destroy_tag - If true, MOAB tag is destroyed with Vec, otherwise it is left on MOAB
466: Output Parameter:
467: . vec - The created vector
469: Level: beginner
471: .keywords: DMMoab, create
472: @*/
473: PetscErrorCode DMMoabCreateVector(DM dm,moab::Tag tag,PetscInt tag_size,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec)
474: {
475: PetscErrorCode ierr;
479: DM_Moab *dmmoab = (DM_Moab*)dm->data;
480: moab::ParallelComm *pcomm = dmmoab->pcomm;
481: moab::Interface *mbiface = dmmoab->mbiface;
482: moab::Tag ltog_tag = dmmoab->ltog_tag;
484: if (!tag && !tag_size) {
485: PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
486: }
487: else {
488: DMMoab_CreateVector(mbiface,pcomm,tag,tag_size,ltog_tag,range,serial,destroy_tag,vec);
489: }
490: return(0);
491: }
496: PetscErrorCode DMMoab_CreateVector(moab::Interface *mbiface,moab::ParallelComm *pcomm,moab::Tag tag,PetscInt tag_size,moab::Tag ltog_tag,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec)
497: {
498: PetscErrorCode ierr;
499: moab::ErrorCode merr;
503: if (!tag) {
504: std::string tag_name;
505: DMMoab_CreateTagName(pcomm,tag_name);
507: // Create the default value for the tag (all zeros):
508: std::vector<PetscScalar> default_value(tag_size, 0.0);
510: // Create the tag:
511: merr = mbiface->tag_get_handle(tag_name.c_str(),tag_size,moab::MB_TYPE_DOUBLE,tag,
512: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,default_value.data());MBERRNM(merr);
513: }
514: else {
516: // Make sure the tag data is of type "double":
517: moab::DataType tag_type;
518: merr = mbiface->tag_get_data_type(tag, tag_type);MBERRNM(merr);
519: if(tag_type != moab::MB_TYPE_DOUBLE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
520: }
522: // Create the MOAB internal data object
523: Vec_MOAB *vmoab;
524: PetscMalloc(sizeof(Vec_MOAB),&vmoab);
525: new (vmoab) Vec_MOAB();
526: vmoab->tag = tag;
527: vmoab->ltog_tag = ltog_tag;
528: vmoab->mbiface = mbiface;
529: vmoab->pcomm = pcomm;
530: vmoab->tag_range = range;
531: vmoab->new_tag = destroy_tag;
532: vmoab->serial = serial;
533: merr = mbiface->tag_get_length(tag,vmoab->tag_size);MBERR("tag_get_size", merr);
535: // Call tag_iterate. This will cause MOAB to allocate memory for the
536: // tag data if it hasn't already happened:
537: int count;
538: void *void_ptr;
539: merr = mbiface->tag_iterate(tag,range.begin(),range.end(),count,void_ptr);MBERRNM(merr);
541: // Check to make sure the tag data is in a single sequence:
542: if ((unsigned)count != range.size()) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only create MOAB Vector for single sequence");
543: PetscScalar *data_ptr = (PetscScalar*)void_ptr;
545: // Create the PETSc Vector:
546: if(!serial) {
547: // This is an MPI Vector:
548: VecCreateMPIWithArray(vmoab->pcomm->comm(),vmoab->tag_size,vmoab->tag_size*range.size(),
549: PETSC_DECIDE,data_ptr,vec);
551: // Vector created, manually set local to global mapping:
552: ISLocalToGlobalMapping ltog;
553: PetscInt *gindices = new PetscInt[range.size()];
554: PetscInt count = 0;
555: moab::Range::iterator iter;
556: for(iter = range.begin(); iter != range.end(); iter++) {
557: int dof;
558: merr = mbiface->tag_get_data(ltog_tag,&(*iter),1,&dof);MBERRNM(merr);
559: gindices[count] = dof;
560: count++;
561: }
563: ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,range.size(),gindices,
564: PETSC_COPY_VALUES,<og);
565: VecSetLocalToGlobalMappingBlock(*vec,ltog);
567: // Clean up:
568: ISLocalToGlobalMappingDestroy(<og);
569: delete [] gindices;
570: } else {
571: // This is a serial vector:
572: VecCreateSeqWithArray(PETSC_COMM_SELF,vmoab->tag_size,vmoab->tag_size*range.size(),data_ptr,vec);
573: }
576: PetscContainer moabdata;
577: PetscContainerCreate(PETSC_COMM_SELF,&moabdata);
578: PetscContainerSetPointer(moabdata,vmoab);
579: PetscContainerSetUserDestroy(moabdata,DMMoab_VecUserDestroy);
580: PetscObjectCompose((PetscObject)*vec,"MOABData",(PetscObject)moabdata);
581: (*vec)->ops->duplicate = DMMoab_VecDuplicate;
583: PetscContainerDestroy(&moabdata);
584: return(0);
585: }
589: /*@
590: DMMoabGetVecTag - Get the MOAB tag associated with this Vec
592: Collective on MPI_Comm
594: Input Parameter:
595: . vec - Vec being queried
597: Output Parameter:
598: . tag - Tag associated with this Vec
600: Level: beginner
602: .keywords: DMMoab, create
603: @*/
604: PetscErrorCode DMMoabGetVecTag(Vec vec,moab::Tag *tag)
605: {
606: PetscContainer moabdata;
607: Vec_MOAB *vmoab;
608: PetscErrorCode ierr;
612: // Get the MOAB private data:
613: PetscObjectQuery((PetscObject)vec,"MOABData", (PetscObject*) &moabdata);
614: PetscContainerGetPointer(moabdata, (void**) &vmoab);
616: *tag = vmoab->tag;
618: return(0);
619: }
624: /*@
625: DMMoabGetVecRange - Get the MOAB entities associated with this Vec
627: Collective on MPI_Comm
629: Input Parameter:
630: . vec - Vec being queried
632: Output Parameter:
633: . range - Entities associated with this Vec
635: Level: beginner
637: .keywords: DMMoab, create
638: @*/
639: PetscErrorCode DMMoabGetVecRange(Vec vec,moab::Range *range)
640: {
641: PetscContainer moabdata;
642: Vec_MOAB *vmoab;
643: PetscErrorCode ierr;
647: // Get the MOAB private data:
648: PetscObjectQuery((PetscObject)vec,"MOABData", (PetscObject*) &moabdata);
649: PetscContainerGetPointer(moabdata, (void**) &vmoab);
651: *range = vmoab->tag_range;
653: return(0);
654: }
659: PetscErrorCode DMMoab_VecDuplicate(Vec x,Vec *y)
660: {
666: // Get the Vec_MOAB struct for the original vector:
667: PetscContainer moabdata;
668: Vec_MOAB *vmoab;
669: PetscObjectQuery((PetscObject)x,"MOABData", (PetscObject*) &moabdata);
670: PetscContainerGetPointer(moabdata, (void**)&vmoab);
672: DMMoab_CreateVector(vmoab->mbiface,vmoab->pcomm,0,vmoab->tag_size,vmoab->ltog_tag,vmoab->tag_range,vmoab->serial,PETSC_TRUE,y);
673: return(0);
674: }
679: /* DMMoab_CreateTagName
680: *
681: * Creates a unique tag name that will be shared across processes. If
682: * pcomm is NULL, then this is a serial vector. A unique tag name
683: * will be returned in tag_name in either case.
684: *
685: * The tag names have the format _PETSC_VEC_N where N is some integer.
686: */
687: PetscErrorCode DMMoab_CreateTagName(const moab::ParallelComm *pcomm,std::string& tag_name)
688: {
689: moab::ErrorCode mberr;
690: PetscErrorCode ierr;
693: const std::string PVEC_PREFIX = "_PETSC_VEC_";
694: const PetscInt PVEC_PREFIX_SIZE = PVEC_PREFIX.size();
696: // Check to see if there are any PETSc vectors defined:
697: const moab::Interface *mbiface = pcomm->get_moab();
698: std::vector<moab::Tag> tags;
699: PetscInt n = 0;
700: mberr = mbiface->tag_get_tags(tags);MBERRNM(mberr);
701: for(unsigned i = 0; i < tags.size(); i++) {
702: std::string s;
703: mberr = mbiface->tag_get_name(tags[i],s);MBERRNM(mberr);
704: if(s.find(PVEC_PREFIX) != std::string::npos){
705: // This tag represents a PETSc vector. Find the vector number:
706: PetscInt m;
707: std::istringstream(s.substr(PVEC_PREFIX_SIZE)) >> m;
708: if(m >= n) n = m+1;
709: }
710: }
712: // Make sure that n is consistent across all processes:
713: PetscInt global_n;
714: MPI_Comm comm = PETSC_COMM_SELF;
715: if(pcomm) comm = pcomm->comm();
716: MPI_Allreduce(&n,&global_n,1,MPI_INT,MPI_MAX,comm);
718: // Set the answer and return:
719: std::ostringstream ss;
720: ss << PVEC_PREFIX << global_n;
721: tag_name = ss.str();
722: return(0);
723: }
728: PetscErrorCode DMMoab_VecUserDestroy(void *user)
729: {
730: Vec_MOAB *vmoab;
731: PetscErrorCode ierr;
732: moab::ErrorCode merr;
735: vmoab = (Vec_MOAB*)user;
736: vmoab->tag_range.~Range();
737: if(vmoab->new_tag) {
738: // Tag created via a call to VecDuplicate, delete the underlying tag in MOAB...
739: merr = vmoab->mbiface->tag_delete(vmoab->tag);MBERRNM(merr);
740: }
742: PetscFree(vmoab);
743: return(0);
744: }