Actual source code: ex1.cxx

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

  3: #include <petscdmmoab.h>
  4: #include "moab/ScdInterface.hpp"

  6: typedef struct {
  7:   DM            dm; /* DM implementation using the MOAB interface */
  8:   PetscLogEvent createMeshEvent;
  9:   /* Domain and mesh definition */
 10:   PetscInt dim;
 11:   char     filename[PETSC_MAX_PATH_LEN];
 12:   char     tagname[PETSC_MAX_PATH_LEN];
 13: } AppCtx;

 15: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 16: {
 17:   PetscBool flg;

 19:   PetscStrcpy(options->filename, "");
 20:   PetscStrcpy(options->tagname, "petsc_tag");
 21:   options->dim = -1;

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

 29:   PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);
 30:   return 0;
 31: }

 33: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 34: {
 35:   moab::Interface *iface    = NULL;
 36:   moab::Tag        tag      = NULL;
 37:   moab::Tag        ltog_tag = NULL;
 38:   moab::Range      range;
 39:   PetscInt         tagsize;
 40:   moab::ErrorCode  merr;

 42:   PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0);
 43:   DMMoabCreateMoab(comm, iface, &ltog_tag, &range, dm);
 44:   std::cout << "Created DMMoab using DMMoabCreateMoab." << std::endl;
 45:   DMMoabGetInterface(*dm, &iface);

 47:   // load file and get entities of requested or max dimension
 48:   if (strlen(user->filename) > 0) {
 49:     merr = iface->load_file(user->filename);
 50:     MBERRNM(merr);
 51:     std::cout << "Read mesh from file " << user->filename << std::endl;
 52:   } else {
 53:     // make a simple structured mesh
 54:     moab::ScdInterface *scdi;
 55:     merr = iface->query_interface(scdi);

 57:     moab::ScdBox *box;
 58:     merr = scdi->construct_box(moab::HomCoord(0, 0, 0), moab::HomCoord(5, 5, 5), NULL, 0, box);
 59:     MBERRNM(merr);
 60:     user->dim = 3;
 61:     merr      = iface->set_dimension(user->dim);
 62:     MBERRNM(merr);
 63:     std::cout << "Created structured 5x5x5 mesh." << std::endl;
 64:   }
 65:   if (-1 == user->dim) {
 66:     moab::Range tmp_range;
 67:     merr = iface->get_entities_by_handle(0, tmp_range);
 68:     MBERRNM(merr);
 69:     if (tmp_range.empty()) MBERRNM(moab::MB_FAILURE);
 70:     user->dim = iface->dimension_from_handle(*tmp_range.rbegin());
 71:   }
 72:   merr = iface->get_entities_by_dimension(0, user->dim, range);
 73:   MBERRNM(merr);
 74:   DMMoabSetLocalVertices(*dm, &range);

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

 94:   DMSetUp(*dm);

 96:   // create the dmmoab and initialize its data
 97:   PetscObjectSetName((PetscObject)*dm, "MOAB mesh");
 98:   PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0);
 99:   user->dm = *dm;
100:   return 0;
101: }

103: int main(int argc, char **argv)
104: {
105:   AppCtx           user; /* user-defined work context */
106:   moab::ErrorCode  merr;
107:   Vec              vec;
108:   moab::Interface *mbImpl  = NULL;
109:   moab::Tag        datatag = NULL;

112:   PetscInitialize(&argc, &argv, NULL, help);
113:   ProcessOptions(PETSC_COMM_WORLD, &user);

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

117:   DMMoabGetInterface(user.dm, &mbImpl);
118:   merr = mbImpl->tag_get_handle(user.tagname, datatag);
119:   MBERRNM(merr);
120:   DMMoabCreateVector(user.dm, datatag, NULL, PETSC_TRUE, PETSC_FALSE, &vec); /* create a vec from user-input tag */

122:   std::cout << "Created VecMoab from existing tag." << std::endl;
123:   VecDestroy(&vec);
124:   std::cout << "Destroyed VecMoab." << std::endl;
125:   DMDestroy(&user.dm);
126:   std::cout << "Destroyed DMMoab." << std::endl;
127:   PetscFinalize();
128:   return 0;
129: }

131: /*TEST

133:      build:
134:        requires: moab

136:      test:

138: TEST*/