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

  3: #include <petscdmmoab.h>
  4: #if defined(__GNUC__) || defined(__GNUG__)
  5:   #pragma GCC diagnostic push
  6:   #pragma GCC diagnostic ignored "-Wunused-result"
  7: #endif
  8: #include "moab/ScdInterface.hpp"
  9: #if defined(__GNUC__) || defined(__GNUG__)
 10:   #pragma GCC diagnostic pop
 11: #endif

 13: typedef struct {
 14:   DM            dm; /* DM implementation using the MOAB interface */
 15:   PetscLogEvent createMeshEvent;
 16:   /* Domain and mesh definition */
 17:   PetscInt dim;
 18:   char     filename[PETSC_MAX_PATH_LEN];
 19:   char     tagname[PETSC_MAX_PATH_LEN];
 20: } AppCtx;

 22: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 23: {
 24:   PetscBool flg;

 26:   PetscFunctionBegin;
 27:   PetscCall(PetscStrncpy(options->filename, "", sizeof(options->filename)));
 28:   PetscCall(PetscStrncpy(options->tagname, "petsc_tag", sizeof(options->tagname)));
 29:   options->dim = -1;

 31:   PetscOptionsBegin(comm, "", "MOAB example options", "DMMOAB");
 32:   PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex1.cxx", options->dim, &options->dim, NULL, PETSC_DECIDE, 3));
 33:   PetscCall(PetscOptionsString("-filename", "The file containing the mesh", "ex1.cxx", options->filename, options->filename, sizeof(options->filename), NULL));
 34:   PetscCall(PetscOptionsString("-tagname", "The tag name from which to create a vector", "ex1.cxx", options->tagname, options->tagname, sizeof(options->tagname), &flg));
 35:   PetscOptionsEnd();

 37:   PetscCall(PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent));
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 42: {
 43:   moab::Interface *iface    = NULL;
 44:   moab::Tag        tag      = NULL;
 45:   moab::Tag        ltog_tag = NULL;
 46:   moab::Range      range;
 47:   PetscInt         tagsize;
 48:   moab::ErrorCode  merr;

 50:   PetscFunctionBegin;
 51:   PetscCall(PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0));
 52:   PetscCall(DMMoabCreateMoab(comm, iface, &ltog_tag, &range, dm));
 53:   std::cout << "Created DMMoab using DMMoabCreateMoab." << std::endl;
 54:   PetscCall(DMMoabGetInterface(*dm, &iface));

 56:   // load file and get entities of requested or max dimension
 57:   if (strlen(user->filename) > 0) {
 58:     merr = iface->load_file(user->filename);
 59:     MBERRNM(merr);
 60:     std::cout << "Read mesh from file " << user->filename << std::endl;
 61:   } else {
 62:     // make a simple structured mesh
 63:     moab::ScdInterface *scdi;
 64:     merr = iface->query_interface(scdi);

 66:     moab::ScdBox *box;
 67:     merr = scdi->construct_box(moab::HomCoord(0, 0, 0), moab::HomCoord(5, 5, 5), NULL, 0, box);
 68:     MBERRNM(merr);
 69:     user->dim = 3;
 70:     merr      = iface->set_dimension(user->dim);
 71:     MBERRNM(merr);
 72:     std::cout << "Created structured 5x5x5 mesh." << std::endl;
 73:   }
 74:   if (-1 == user->dim) {
 75:     moab::Range tmp_range;
 76:     merr = iface->get_entities_by_handle(0, tmp_range);
 77:     MBERRNM(merr);
 78:     if (tmp_range.empty()) MBERRNM(moab::MB_FAILURE);
 79:     user->dim = iface->dimension_from_handle(*tmp_range.rbegin());
 80:   }
 81:   merr = iface->get_entities_by_dimension(0, user->dim, range);
 82:   MBERRNM(merr);
 83:   PetscCall(DMMoabSetLocalVertices(*dm, &range));

 85:   // get the requested tag and create if necessary
 86:   std::cout << "Creating tag with name: " << user->tagname << ";\n";
 87:   merr = iface->tag_get_handle(user->tagname, 1, moab::MB_TYPE_DOUBLE, tag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE);
 88:   MBERRNM(merr);
 89:   {
 90:     // initialize new tag with gids
 91:     std::vector<double> tag_vals(range.size());
 92:     merr = iface->tag_get_data(ltog_tag, range, tag_vals.data());
 93:     MBERRNM(merr); // read them into ints
 94:     double *dval = tag_vals.data();
 95:     int    *ival = reinterpret_cast<int *>(dval); // "stretch" them into doubles, from the end
 96:     for (int i = tag_vals.size() - 1; i >= 0; i--) dval[i] = ival[i];
 97:     merr = iface->tag_set_data(tag, range, tag_vals.data());
 98:     MBERRNM(merr); // write them into doubles
 99:   }
100:   merr = iface->tag_get_length(tag, tagsize);
101:   MBERRNM(merr);

103:   PetscCall(DMSetUp(*dm));

105:   // create the dmmoab and initialize its data
106:   PetscCall(PetscObjectSetName((PetscObject)*dm, "MOAB mesh"));
107:   PetscCall(PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0));
108:   user->dm = *dm;
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }

112: int main(int argc, char **argv)
113: {
114:   AppCtx           user; /* user-defined work context */
115:   moab::ErrorCode  merr;
116:   Vec              vec;
117:   moab::Interface *mbImpl  = NULL;
118:   moab::Tag        datatag = NULL;

120:   PetscFunctionBeginUser;
121:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
122:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));

124:   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &user.dm)); /* create the MOAB dm and the mesh */

126:   PetscCall(DMMoabGetInterface(user.dm, &mbImpl));
127:   merr = mbImpl->tag_get_handle(user.tagname, datatag);
128:   MBERRNM(merr);
129:   PetscCall(DMMoabCreateVector(user.dm, datatag, NULL, PETSC_TRUE, PETSC_FALSE, &vec)); /* create a vec from user-input tag */

131:   std::cout << "Created VecMoab from existing tag." << std::endl;
132:   PetscCall(VecDestroy(&vec));
133:   std::cout << "Destroyed VecMoab." << std::endl;
134:   PetscCall(DMDestroy(&user.dm));
135:   std::cout << "Destroyed DMMoab." << std::endl;
136:   PetscCall(PetscFinalize());
137:   return 0;
138: }

140: /*TEST

142:      build:
143:        requires: moab !complex

145:      test:

147: TEST*/