1 static char help[] = "Create a box mesh with DMMoab and test defining a tag on the mesh\n\n";
2
3 #include <petscdmmoab.h>
4
5 typedef struct {
6 DM dm; /* DM implementation using the MOAB interface */
7 PetscBool debug; /* The debugging level */
8 PetscLogEvent createMeshEvent;
9 /* Domain and mesh definition */
10 PetscInt dim; /* The topological mesh dimension */
11 PetscInt nele; /* Elements in each dimension */
12 PetscInt degree; /* Degree of refinement */
13 PetscBool simplex; /* Use simplex elements */
14 PetscInt nlevels; /* Number of levels in mesh hierarchy */
15 PetscInt nghost; /* Number of ghost layers in the mesh */
16 char input_file[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
17 char output_file[PETSC_MAX_PATH_LEN]; /* Output mesh file name */
18 PetscBool write_output; /* Write output mesh and data to file */
19 } AppCtx;
20
ProcessOptions(MPI_Comm comm,AppCtx * options)21 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
22 {
23 PetscFunctionBegin;
24 options->debug = PETSC_FALSE;
25 options->nlevels = 1;
26 options->nghost = 1;
27 options->dim = 2;
28 options->nele = 5;
29 options->degree = 2;
30 options->simplex = PETSC_FALSE;
31 options->write_output = PETSC_FALSE;
32 options->input_file[0] = '\0';
33 PetscCall(PetscStrncpy(options->output_file, "ex3.h5m", sizeof(options->output_file)));
34
35 PetscOptionsBegin(comm, "", "Uniform Mesh Refinement Options", "DMMOAB");
36 PetscCall(PetscOptionsBool("-debug", "Enable debug messages", "ex2.cxx", options->debug, &options->debug, NULL));
37 PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex3.cxx", options->dim, &options->dim, NULL, 0, 3));
38 PetscCall(PetscOptionsBoundedInt("-n", "The number of elements in each dimension", "ex3.cxx", options->nele, &options->nele, NULL, 1));
39 PetscCall(PetscOptionsBoundedInt("-levels", "Number of levels in the hierarchy", "ex3.cxx", options->nlevels, &options->nlevels, NULL, 0));
40 PetscCall(PetscOptionsBoundedInt("-degree", "Number of degrees at each level of refinement", "ex3.cxx", options->degree, &options->degree, NULL, 0));
41 PetscCall(PetscOptionsBoundedInt("-ghost", "Number of ghost layers in the mesh", "ex3.cxx", options->nghost, &options->nghost, NULL, 0));
42 PetscCall(PetscOptionsBool("-simplex", "Create simplices instead of tensor product elements", "ex3.cxx", options->simplex, &options->simplex, NULL));
43 PetscCall(PetscOptionsString("-input", "The input mesh file", "ex3.cxx", options->input_file, options->input_file, sizeof(options->input_file), NULL));
44 PetscCall(PetscOptionsString("-io", "Write out the mesh and solution that is defined on it (Default H5M format)", "ex3.cxx", options->output_file, options->output_file, sizeof(options->output_file), &options->write_output));
45 PetscOptionsEnd();
46
47 PetscCall(PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent));
48 PetscFunctionReturn(PETSC_SUCCESS);
49 }
50
CreateMesh(MPI_Comm comm,AppCtx * user)51 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user)
52 {
53 size_t len;
54 PetscMPIInt rank;
55
56 PetscFunctionBegin;
57 PetscCall(PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0));
58 PetscCallMPI(MPI_Comm_rank(comm, &rank));
59 PetscCall(PetscStrlen(user->input_file, &len));
60 if (len) {
61 if (user->debug) PetscCall(PetscPrintf(comm, "Loading mesh from file: %s and creating the coarse level DM object.\n", user->input_file));
62 PetscCall(DMMoabLoadFromFile(comm, user->dim, user->nghost, user->input_file, "", &user->dm));
63 } else {
64 if (user->debug)
65 PetscCall(PetscPrintf(comm, "Creating a %" PetscInt_FMT "-dimensional structured %s mesh of %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " in memory and creating a DM object.\n", user->dim, user->simplex ? "simplex" : "regular", user->nele,
66 user->nele, user->nele));
67 PetscCall(DMMoabCreateBoxMesh(comm, user->dim, user->simplex, NULL, user->nele, user->nghost, &user->dm));
68 }
69 PetscCall(PetscObjectSetName((PetscObject)user->dm, "Coarse Mesh"));
70 PetscCall(PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0));
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
main(int argc,char ** argv)74 int main(int argc, char **argv)
75 {
76 AppCtx user; /* user-defined work context */
77 MPI_Comm comm;
78 PetscInt i;
79 Mat R;
80 DM *dmhierarchy;
81 PetscInt *degrees;
82 PetscBool createR;
83
84 PetscFunctionBeginUser;
85 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
86 comm = PETSC_COMM_WORLD;
87 createR = PETSC_FALSE;
88
89 PetscCall(ProcessOptions(comm, &user));
90 PetscCall(CreateMesh(comm, &user));
91 PetscCall(DMSetFromOptions(user.dm));
92
93 /* SetUp the data structures for DMMOAB */
94 PetscCall(DMSetUp(user.dm));
95
96 PetscCall(PetscMalloc(sizeof(DM) * (user.nlevels + 1), &dmhierarchy));
97 for (i = 0; i <= user.nlevels; i++) dmhierarchy[i] = NULL;
98
99 // coarsest grid = 0
100 // finest grid = nlevels
101 dmhierarchy[0] = user.dm;
102 PetscCall(PetscObjectReference((PetscObject)user.dm));
103
104 if (user.nlevels) {
105 PetscCall(PetscMalloc1(user.nlevels, °rees));
106 for (i = 0; i < user.nlevels; i++) degrees[i] = user.degree;
107 if (user.debug) PetscCall(PetscPrintf(comm, "Generate the MOAB mesh hierarchy with %" PetscInt_FMT " levels.\n", user.nlevels));
108 PetscCall(DMMoabGenerateHierarchy(user.dm, user.nlevels, degrees));
109
110 PetscBool usehierarchy = PETSC_FALSE;
111 if (usehierarchy) PetscCall(DMRefineHierarchy(user.dm, user.nlevels, &dmhierarchy[1]));
112 else {
113 if (user.debug) {
114 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Level %" PetscInt_FMT "\n", 0));
115 PetscCall(DMView(user.dm, 0));
116 }
117 for (i = 1; i <= user.nlevels; i++) {
118 if (user.debug) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Level %" PetscInt_FMT "\n", i));
119 PetscCall(DMRefine(dmhierarchy[i - 1], MPI_COMM_NULL, &dmhierarchy[i]));
120 if (createR) PetscCall(DMCreateInterpolation(dmhierarchy[i - 1], dmhierarchy[i], &R, NULL));
121 if (user.debug) {
122 PetscCall(DMView(dmhierarchy[i], 0));
123 if (createR) PetscCall(MatView(R, 0));
124 }
125 /* Solvers could now set operator "R" to the multigrid PC object for level i
126 PCMGSetInterpolation(pc,i,R)
127 */
128 if (createR) PetscCall(MatDestroy(&R));
129 }
130 }
131 }
132
133 if (user.write_output) {
134 if (user.debug) PetscCall(PetscPrintf(comm, "Output mesh hierarchy to file: %s.\n", user.output_file));
135 PetscCall(DMMoabOutput(dmhierarchy[user.nlevels], (const char *)user.output_file, ""));
136 }
137
138 for (i = 0; i <= user.nlevels; i++) PetscCall(DMDestroy(&dmhierarchy[i]));
139 PetscCall(PetscFree(degrees));
140 PetscCall(PetscFree(dmhierarchy));
141 PetscCall(DMDestroy(&user.dm));
142 PetscCall(PetscFinalize());
143 return 0;
144 }
145
146 /*TEST
147
148 build:
149 requires: moab !complex
150
151 test:
152 args: -debug -n 2 -dim 2 -levels 2 -simplex
153 filter: grep -v "DM_0x"
154
155 test:
156 args: -debug -n 2 -dim 3 -levels 2
157 filter: grep -v "DM_0x"
158 suffix: 1_2
159
160 test:
161 args: -debug -n 2 -dim 3 -ghost 1 -levels 2
162 filter: grep -v "DM_0x"
163 nsize: 2
164 suffix: 2_1
165
166 TEST*/
167