xref: /petsc/src/dm/impls/moab/tests/ex2.cxx (revision 58d68138c660dfb4e9f5b03334792cd4f2ffd7cc)
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   PetscBool     simplex; /* Use simplex elements */
13   PetscBool     interlace;
14   char          input_file[PETSC_MAX_PATH_LEN];  /* Import mesh from file */
15   char          output_file[PETSC_MAX_PATH_LEN]; /* Output mesh file name */
16   PetscBool     write_output;                    /* Write output mesh and data to file */
17   PetscInt      nfields;                         /* Number of fields */
18   char         *fieldnames[PETSC_MAX_PATH_LEN];  /* Name of a defined field on the mesh */
19 } AppCtx;
20 
21 PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {
22   PetscBool flg;
23 
24   PetscFunctionBegin;
25   options->debug         = PETSC_FALSE;
26   options->dim           = 2;
27   options->nele          = 5;
28   options->nfields       = 256;
29   options->simplex       = PETSC_FALSE;
30   options->write_output  = PETSC_FALSE;
31   options->interlace     = PETSC_FALSE;
32   options->input_file[0] = '\0';
33   PetscCall(PetscStrcpy(options->output_file, "ex2.h5m"));
34 
35   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMMOAB");
36   PetscCall(PetscOptionsBool("-debug", "Enable debug messages", "ex2.cxx", options->debug, &options->debug, NULL));
37   PetscCall(PetscOptionsBool("-interlace", "Use interlaced arrangement for the field data", "ex2.cxx", options->interlace, &options->interlace, NULL));
38   PetscCall(PetscOptionsBool("-simplex", "Create simplices instead of tensor product elements", "ex2.cxx", options->simplex, &options->simplex, NULL));
39   PetscCall(PetscOptionsRangeInt("-dim", "The topological mesh dimension", "ex2.cxx", options->dim, &options->dim, NULL, 1, 3));
40   PetscCall(PetscOptionsBoundedInt("-n", "The number of elements in each dimension", "ex2.cxx", options->nele, &options->nele, NULL, 1));
41   PetscCall(PetscOptionsString("-meshfile", "The input mesh file", "ex2.cxx", options->input_file, options->input_file, sizeof(options->input_file), NULL));
42   PetscCall(PetscOptionsString("-io", "Write out the mesh and solution that is defined on it (Default H5M format)", "ex2.cxx", options->output_file, options->output_file, sizeof(options->output_file), &options->write_output));
43   PetscCall(PetscOptionsStringArray("-fields", "The list of names of the field variables", "ex2.cxx", options->fieldnames, &options->nfields, &flg));
44   PetscOptionsEnd();
45 
46   if (options->debug) PetscCall(PetscPrintf(comm, "Total number of fields: %" PetscInt_FMT ".\n", options->nfields));
47   if (!flg) { /* if no field names were given by user, assign a default */
48     options->nfields = 1;
49     PetscCall(PetscStrallocpy("TestEX2Var", &options->fieldnames[0]));
50   }
51 
52   PetscCall(PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent));
53   PetscFunctionReturn(0);
54 }
55 
56 PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user) {
57   PetscInt    i;
58   size_t      len;
59   PetscMPIInt rank;
60 
61   PetscFunctionBegin;
62   PetscCall(PetscLogEventBegin(user->createMeshEvent, 0, 0, 0, 0));
63   PetscCallMPI(MPI_Comm_rank(comm, &rank));
64   PetscCall(PetscStrlen(user->input_file, &len));
65   if (len) {
66     if (user->debug) PetscCall(PetscPrintf(comm, "Loading mesh from file: %s and creating a DM object.\n", user->input_file));
67     PetscCall(DMMoabLoadFromFile(comm, user->dim, 1, user->input_file, "", &user->dm));
68   } else {
69     if (user->debug) {
70       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,
71                             user->nele, user->nele));
72     }
73     PetscCall(DMMoabCreateBoxMesh(comm, user->dim, user->simplex, NULL, user->nele, 1, &user->dm));
74   }
75 
76   if (user->debug) {
77     PetscPrintf(comm, "Setting field names to DM: \n");
78     for (i = 0; i < user->nfields; i++) PetscPrintf(comm, "\t Field{%" PetscInt_FMT "} = %s.\n", i, user->fieldnames[i]);
79   }
80   PetscCall(DMMoabSetFieldNames(user->dm, user->nfields, (const char **)user->fieldnames));
81   PetscCall(PetscObjectSetName((PetscObject)user->dm, "Structured Mesh"));
82   PetscCall(PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0));
83   PetscFunctionReturn(0);
84 }
85 
86 int main(int argc, char **argv) {
87   AppCtx      user; /* user-defined work context */
88   PetscRandom rctx;
89   Vec         solution;
90   Mat         system;
91   MPI_Comm    comm;
92   PetscInt    i;
93 
94   PetscFunctionBeginUser;
95   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
96   comm = PETSC_COMM_WORLD;
97   PetscCall(ProcessOptions(comm, &user));
98   PetscCall(CreateMesh(comm, &user));
99 
100   /* set block size */
101   PetscCall(DMMoabSetBlockSize(user.dm, (user.interlace ? user.nfields : 1)));
102   PetscCall(DMSetMatType(user.dm, MATAIJ));
103 
104   PetscCall(DMSetFromOptions(user.dm));
105 
106   /* SetUp the data structures for DMMOAB */
107   PetscCall(DMSetUp(user.dm));
108 
109   if (user.debug) PetscCall(PetscPrintf(comm, "Creating a global vector defined on DM and setting random data.\n"));
110   PetscCall(DMCreateGlobalVector(user.dm, &solution));
111   PetscCall(PetscRandomCreate(comm, &rctx));
112   PetscCall(VecSetRandom(solution, rctx));
113 
114   /* test if matrix allocation for the prescribed matrix type is done correctly */
115   if (user.debug) PetscCall(PetscPrintf(comm, "Creating a global matrix defined on DM with the right block structure.\n"));
116   PetscCall(DMCreateMatrix(user.dm, &system));
117 
118   if (user.write_output) {
119     PetscCall(DMMoabSetGlobalFieldVector(user.dm, solution));
120     if (user.debug) PetscPrintf(comm, "Output mesh and associated field data to file: %s.\n", user.output_file);
121     PetscCall(DMMoabOutput(user.dm, (const char *)user.output_file, ""));
122   }
123 
124   if (user.nfields) {
125     for (i = 0; i < user.nfields; i++) PetscCall(PetscFree(user.fieldnames[i]));
126   }
127   PetscCall(PetscRandomDestroy(&rctx));
128   PetscCall(VecDestroy(&solution));
129   PetscCall(MatDestroy(&system));
130   PetscCall(DMDestroy(&user.dm));
131   PetscCall(PetscFinalize());
132   return 0;
133 }
134 
135 /*TEST
136 
137    build:
138      requires: moab
139 
140    test:
141      args: -debug -fields v1,v2,v3
142 
143 TEST*/
144