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