1 static char help[] = "Simple MOAB example\n\n"; 2 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 12 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; 21 22 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 23 { 24 PetscBool flg; 25 26 PetscFunctionBegin; 27 PetscCall(PetscStrncpy(options->filename, "", sizeof(options->filename))); 28 PetscCall(PetscStrncpy(options->tagname, "petsc_tag", sizeof(options->tagname))); 29 options->dim = -1; 30 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(); 36 37 PetscCall(PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent)); 38 PetscFunctionReturn(PETSC_SUCCESS); 39 } 40 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; 49 50 PetscFunctionBegin; 51 PetscCall(PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0)); 52 PetscCall(DMMoabCreateMoab(comm, iface, <og_tag, &range, dm)); 53 std::cout << "Created DMMoab using DMMoabCreateMoab." << std::endl; 54 PetscCall(DMMoabGetInterface(*dm, &iface)); 55 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); 65 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)); 84 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); 102 103 PetscCall(DMSetUp(*dm)); 104 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 } 111 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; 119 120 PetscFunctionBeginUser; 121 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 122 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 123 124 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &user.dm)); /* create the MOAB dm and the mesh */ 125 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 */ 130 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 } 139 140 /*TEST 141 142 build: 143 requires: moab !complex 144 145 test: 146 147 TEST*/ 148