xref: /petsc/src/dm/impls/moab/tests/ex3.cxx (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
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 
21 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
22 {
23   PetscErrorCode ierr;
24 
25   PetscFunctionBegin;
26   options->debug             = PETSC_FALSE;
27   options->nlevels           = 1;
28   options->nghost            = 1;
29   options->dim               = 2;
30   options->nele              = 5;
31   options->degree            = 2;
32   options->simplex           = PETSC_FALSE;
33   options->write_output      = PETSC_FALSE;
34   options->input_file[0]     = '\0';
35   CHKERRQ(PetscStrcpy(options->output_file,"ex3.h5m"));
36 
37   ierr = PetscOptionsBegin(comm, "", "Uniform Mesh Refinement Options", "DMMOAB");CHKERRQ(ierr);
38   CHKERRQ(PetscOptionsBool("-debug", "Enable debug messages", "ex2.cxx", options->debug, &options->debug, NULL));
39   CHKERRQ(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex3.cxx", options->dim, &options->dim, NULL,0,3));
40   CHKERRQ(PetscOptionsBoundedInt("-n", "The number of elements in each dimension", "ex3.cxx", options->nele, &options->nele, NULL,1));
41   CHKERRQ(PetscOptionsBoundedInt("-levels", "Number of levels in the hierarchy", "ex3.cxx", options->nlevels, &options->nlevels, NULL,0));
42   CHKERRQ(PetscOptionsBoundedInt("-degree", "Number of degrees at each level of refinement", "ex3.cxx", options->degree, &options->degree, NULL,0));
43   CHKERRQ(PetscOptionsBoundedInt("-ghost", "Number of ghost layers in the mesh", "ex3.cxx", options->nghost, &options->nghost, NULL,0));
44   CHKERRQ(PetscOptionsBool("-simplex", "Create simplices instead of tensor product elements", "ex3.cxx", options->simplex, &options->simplex, NULL));
45   CHKERRQ(PetscOptionsString("-input", "The input mesh file", "ex3.cxx", options->input_file, options->input_file, sizeof(options->input_file), NULL));
46   CHKERRQ(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));
47   ierr = PetscOptionsEnd();CHKERRQ(ierr);
48 
49   CHKERRQ(PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent));
50   PetscFunctionReturn(0);
51 };
52 
53 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user)
54 {
55   size_t         len;
56   PetscMPIInt    rank;
57 
58   PetscFunctionBegin;
59   CHKERRQ(PetscLogEventBegin(user->createMeshEvent,0,0,0,0));
60   CHKERRMPI(MPI_Comm_rank(comm, &rank));
61   CHKERRQ(PetscStrlen(user->input_file, &len));
62   if (len) {
63     if (user->debug) CHKERRQ(PetscPrintf(comm, "Loading mesh from file: %s and creating the coarse level DM object.\n",user->input_file));
64     CHKERRQ(DMMoabLoadFromFile(comm, user->dim, user->nghost, user->input_file, "", &user->dm));
65   } else {
66     if (user->debug) CHKERRQ(PetscPrintf(comm, "Creating a %D-dimensional structured %s mesh of %Dx%Dx%D in memory and creating a DM object.\n",user->dim,(user->simplex?"simplex":"regular"),user->nele,user->nele,user->nele));
67     CHKERRQ(DMMoabCreateBoxMesh(comm, user->dim, user->simplex, NULL, user->nele, user->nghost, &user->dm));
68   }
69   CHKERRQ(PetscObjectSetName((PetscObject)user->dm, "Coarse Mesh"));
70   CHKERRQ(PetscLogEventEnd(user->createMeshEvent,0,0,0,0));
71   PetscFunctionReturn(0);
72 }
73 
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   PetscErrorCode ierr;
84 
85   ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
86   comm = PETSC_COMM_WORLD;
87   createR = PETSC_FALSE;
88 
89   CHKERRQ(ProcessOptions(comm, &user));
90   CHKERRQ(CreateMesh(comm, &user));
91   CHKERRQ(DMSetFromOptions(user.dm));
92 
93   /* SetUp the data structures for DMMOAB */
94   CHKERRQ(DMSetUp(user.dm));
95 
96   CHKERRQ(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   PetscObjectReference((PetscObject)user.dm);
103 
104   if (user.nlevels) {
105     CHKERRQ(PetscMalloc1(user.nlevels, &degrees));
106     for (i=0; i < user.nlevels; i++) degrees[i] = user.degree;
107     if (user.debug) CHKERRQ(PetscPrintf(comm, "Generate the MOAB mesh hierarchy with %D levels.\n", user.nlevels));
108     CHKERRQ(DMMoabGenerateHierarchy(user.dm,user.nlevels,degrees));
109 
110     PetscBool usehierarchy=PETSC_FALSE;
111     if (usehierarchy) CHKERRQ(DMRefineHierarchy(user.dm,user.nlevels,&dmhierarchy[1]));
112     else {
113       if (user.debug) {
114         CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Level %D\n", 0));
115         CHKERRQ(DMView(user.dm, 0));
116       }
117       for (i=1; i<=user.nlevels; i++) {
118         if (user.debug) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Level %D\n", i));
119         CHKERRQ(DMRefine(dmhierarchy[i-1],MPI_COMM_NULL,&dmhierarchy[i]));
120         if (createR) CHKERRQ(DMCreateInterpolation(dmhierarchy[i-1],dmhierarchy[i],&R,NULL));
121         if (user.debug) {
122           CHKERRQ(DMView(dmhierarchy[i], 0));
123           if (createR) CHKERRQ(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) CHKERRQ(MatDestroy(&R));
129       }
130     }
131   }
132 
133   if (user.write_output) {
134     if (user.debug) CHKERRQ(PetscPrintf(comm, "Output mesh hierarchy to file: %s.\n",user.output_file));
135     CHKERRQ(DMMoabOutput(dmhierarchy[user.nlevels],(const char*)user.output_file,""));
136   }
137 
138   for (i=0; i<=user.nlevels; i++) CHKERRQ(DMDestroy(&dmhierarchy[i]));
139   CHKERRQ(PetscFree(degrees));
140   CHKERRQ(PetscFree(dmhierarchy));
141   CHKERRQ(DMDestroy(&user.dm));
142   ierr = PetscFinalize();
143   return 0;
144 }
145 
146 /*TEST
147 
148      build:
149        requires: moab
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