1: static char help[] = "Simple MOAB example\n\n";

  3: #include <petscdmmoab.h>
  4: #include <iostream>
  5: #include "moab/Interface.hpp"
  6: #include "moab/ScdInterface.hpp"
  7: #include "MBTagConventions.hpp"

  9: class AppCtx {
 10: public:
 11:   DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
 12:   PetscLogEvent createMeshEvent;
 13:   /* Domain and mesh definition */
 14:   PetscInt dim;
 15:   char filename[PETSC_MAX_PATH_LEN];
 16:   char tagname[PETSC_MAX_PATH_LEN];
 17:   moab::Interface *iface;
 18:   moab::ParallelComm *pcomm;
 19:   moab::Tag tag;
 20:   moab::Tag ltog_tag;
 21:   moab::Range range;
 22:   PetscInt tagsize;

 24:   AppCtx()
 25:           : dm(NULL), dim(-1), iface(NULL), pcomm(NULL), tag(0), ltog_tag(0), tagsize(0)
 26:       {strcpy(filename, ""); strcpy(tagname, "");}

 28: };

 32: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 33: {

 37:   strcpy(options->filename, "");
 38:   strcpy(options->tagname, "");
 39:   options->dim = -1;

 41:   PetscOptionsBegin(comm, "", "MOAB example options", "DMMOAB");
 42:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);
 43:   PetscOptionsString("-filename", "The file containing the mesh", "ex1.c", options->filename, options->filename, sizeof(options->filename), NULL);
 44:   PetscOptionsString("-tagname", "The tag name from which to create a vector", "ex1.c", options->tagname, options->tagname, sizeof(options->tagname), NULL);
 45:   PetscOptionsEnd();

 47:   PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);
 48:   return(0);
 49: };

 53: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 54: {

 58:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
 59:   DMMoabCreateMoab(comm, user->iface, user->pcomm, user->ltog_tag, &user->range, dm);
 60:   std::cout << "Created DMMoab using DMMoabCreateDMAndInstance." << std::endl;
 61:   DMMoabGetInterface(*dm, &user->iface);

 63:     // load file and get entities of requested or max dimension
 64:   moab::ErrorCode merr;
 65:   if (strlen(user->filename) > 0) {
 66:     merr = user->iface->load_file(user->filename);MBERRNM(merr);
 67:     std::cout << "Read mesh from file " << user->filename << std::endl;
 68:   }
 69:   else {
 70:       // make a simple structured mesh
 71:     moab::ScdInterface *scdi;
 72:     merr = user->iface->query_interface(scdi);
 73:     moab::ScdBox *box;
 74:     merr = scdi->construct_box (moab::HomCoord(0,0,0), moab::HomCoord(2,2,2), NULL, 0, box);MBERRNM(merr);
 75:     user->dim = 3;
 76:     std::cout << "Created structured 2x2x2 mesh." << std::endl;
 77:   }
 78:   if (-1 == user->dim) {
 79:     moab::Range tmp_range;
 80:     merr = user->iface->get_entities_by_handle(0, tmp_range);MBERRNM(merr);
 81:     if (tmp_range.empty()) {
 82:       MBERRNM(moab::MB_FAILURE);
 83:     }
 84:     user->dim = user->iface->dimension_from_handle(*tmp_range.rbegin());
 85:   }
 86:   merr = user->iface->get_entities_by_dimension(0, user->dim, user->range);MBERRNM(merr);
 87:   DMMoabSetRange(*dm, user->range);

 89:     // get the requested tag if a name was input
 90:   if (strlen(user->tagname)) {
 91:     merr = user->iface->tag_get_handle(user->tagname, user->tag);MBERRNM(merr);
 92:     moab::DataType ttype;
 93:     merr = user->iface->tag_get_data_type(user->tag, ttype);MBERRNM(merr);
 94:     if (ttype != moab::MB_TYPE_DOUBLE) {
 95:       printf("Tag type must be double!.\n");
 96:       return 1;
 97:     }
 98:   }
 99:   else {
100:       // make a new tag
101:     merr = user->iface->tag_get_handle("petsc_tag", 1, moab::MB_TYPE_DOUBLE, user->tag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE); MBERRNM(merr);
102:       // initialize new tag with gids
103:     std::vector<double> tag_vals(user->range.size());
104:     moab::Tag gid_tag;
105:     merr = user->iface->tag_get_handle("GLOBAL_ID", gid_tag);MBERRNM(merr);
106:     merr = user->iface->tag_get_data(gid_tag, user->range, tag_vals.data());MBERRNM(merr); // read them into ints
107:     double *dval = tag_vals.data(); int *ival = reinterpret_cast<int*>(dval); // "stretch" them into doubles, from the end
108:     for (int i = tag_vals.size()-1; i >= 0; i--) dval[i] = ival[i];
109:     merr = user->iface->tag_set_data(user->tag, user->range, tag_vals.data());MBERRNM(merr); // write them into doubles
110:   }
111:   merr = user->iface->tag_get_length(user->tag, user->tagsize);MBERRNM(merr);

113:     // create the dmmoab and initialize its data
114:   PetscObjectSetName((PetscObject) *dm, "MOAB mesh");
115:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
116:   user->dm = *dm;
117:   return(0);
118: }

122: int main(int argc, char **argv)
123: {
124:   AppCtx         user;                 /* user-defined work context */
126:   Vec            vec;

128:   PetscInitialize(&argc, &argv, NULL, help);
129:   ProcessOptions(PETSC_COMM_WORLD, &user);

131:   CreateMesh(PETSC_COMM_WORLD, &user, &user.dm); /* create the MOAB dm and the mesh */
132:   DMMoabCreateVector(user.dm, user.tag, 1, user.range, PETSC_TRUE, PETSC_FALSE,
133:                               &vec); /* create a vec from user-input tag */
134:   std::cout << "Created VecMoab from existing tag." << std::endl;
135:   VecDestroy(&vec);
136:   std::cout << "Destroyed VecMoab." << std::endl;
137:   DMDestroy(&user.dm);
138:   std::cout << "Destroyed DMMoab." << std::endl;
139:   PetscFinalize();
140:   return 0;
141: }