1 static char help[] = "Simple MOAB example\n\n"; 2 3 #include <petscdmmoab.h> 4 #include "moab/ScdInterface.hpp" 5 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; 14 15 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) { 16 PetscBool flg; 17 18 PetscFunctionBegin; 19 PetscCall(PetscStrcpy(options->filename, "")); 20 PetscCall(PetscStrcpy(options->tagname, "petsc_tag")); 21 options->dim = -1; 22 23 PetscOptionsBegin(comm, "", "MOAB example options", "DMMOAB"); 24 PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex1.cxx", options->dim, &options->dim, NULL, PETSC_DECIDE, 3)); 25 PetscCall(PetscOptionsString("-filename", "The file containing the mesh", "ex1.cxx", options->filename, options->filename, sizeof(options->filename), NULL)); 26 PetscCall(PetscOptionsString("-tagname", "The tag name from which to create a vector", "ex1.cxx", options->tagname, options->tagname, sizeof(options->tagname), &flg)); 27 PetscOptionsEnd(); 28 29 PetscCall(PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent)); 30 PetscFunctionReturn(0); 31 } 32 33 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) { 34 moab::Interface *iface = NULL; 35 moab::Tag tag = NULL; 36 moab::Tag ltog_tag = NULL; 37 moab::Range range; 38 PetscInt tagsize; 39 moab::ErrorCode merr; 40 41 PetscFunctionBegin; 42 PetscCall(PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0)); 43 PetscCall(DMMoabCreateMoab(comm, iface, <og_tag, &range, dm)); 44 std::cout << "Created DMMoab using DMMoabCreateMoab." << std::endl; 45 PetscCall(DMMoabGetInterface(*dm, &iface)); 46 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); 56 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 PetscCall(DMMoabSetLocalVertices(*dm, &range)); 75 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); 93 94 PetscCall(DMSetUp(*dm)); 95 96 // create the dmmoab and initialize its data 97 PetscCall(PetscObjectSetName((PetscObject)*dm, "MOAB mesh")); 98 PetscCall(PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0)); 99 user->dm = *dm; 100 PetscFunctionReturn(0); 101 } 102 103 int main(int argc, char **argv) { 104 AppCtx user; /* user-defined work context */ 105 moab::ErrorCode merr; 106 Vec vec; 107 moab::Interface *mbImpl = NULL; 108 moab::Tag datatag = NULL; 109 110 PetscFunctionBeginUser; 111 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 112 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 113 114 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &user.dm)); /* create the MOAB dm and the mesh */ 115 116 PetscCall(DMMoabGetInterface(user.dm, &mbImpl)); 117 merr = mbImpl->tag_get_handle(user.tagname, datatag); 118 MBERRNM(merr); 119 PetscCall(DMMoabCreateVector(user.dm, datatag, NULL, PETSC_TRUE, PETSC_FALSE, &vec)); /* create a vec from user-input tag */ 120 121 std::cout << "Created VecMoab from existing tag." << std::endl; 122 PetscCall(VecDestroy(&vec)); 123 std::cout << "Destroyed VecMoab." << std::endl; 124 PetscCall(DMDestroy(&user.dm)); 125 std::cout << "Destroyed DMMoab." << std::endl; 126 PetscCall(PetscFinalize()); 127 return 0; 128 } 129 130 /*TEST 131 132 build: 133 requires: moab 134 135 test: 136 137 TEST*/ 138