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: {
 18:   PetscBool      flg;

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

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

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

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

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

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

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

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

 88:   DMSetUp(*dm);

 90:     // create the dmmoab and initialize its data
 91:   PetscObjectSetName((PetscObject) *dm, "MOAB mesh");
 92:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
 93:   user->dm = *dm;
 94:   return 0;
 95: }

 97: int main(int argc, char **argv)
 98: {
 99:   AppCtx            user;                 /* user-defined work context */
100:   moab::ErrorCode   merr;
101:   Vec               vec;
102:   moab::Interface*  mbImpl=NULL;
103:   moab::Tag         datatag=NULL;

105:   PetscInitialize(&argc, &argv, NULL,help);
106:   ProcessOptions(PETSC_COMM_WORLD, &user);

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

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

114:   std::cout << "Created VecMoab from existing tag." << std::endl;
115:   VecDestroy(&vec);
116:   std::cout << "Destroyed VecMoab." << std::endl;
117:   DMDestroy(&user.dm);
118:   std::cout << "Destroyed DMMoab." << std::endl;
119:   PetscFinalize();
120:   return 0;
121: }

123: /*TEST

125:      build:
126:        requires: moab

128:      test:

130: TEST*/