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
ProcessOptions(MPI_Comm comm,AppCtx * options)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
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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
main(int argc,char ** argv)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