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
ProcessOptions(MPI_Comm comm,AppCtx * options)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(PetscStrncpy(options->output_file, "ex2.h5m", sizeof(options->output_file)));
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(PETSC_SUCCESS);
55 }
56
CreateMesh(MPI_Comm comm,AppCtx * user)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 } else {
71 if (user->debug) {
72 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,
73 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 PetscCall(PetscPrintf(comm, "Setting field names to DM: \n"));
80 for (i = 0; i < user->nfields; i++) PetscCall(PetscPrintf(comm, "\t Field{%" PetscInt_FMT "} = %s.\n", i, user->fieldnames[i]));
81 }
82 PetscCall(DMMoabSetFieldNames(user->dm, user->nfields, (const char **)user->fieldnames));
83 PetscCall(PetscObjectSetName((PetscObject)user->dm, "Structured Mesh"));
84 PetscCall(PetscLogEventEnd(user->createMeshEvent, 0, 0, 0, 0));
85 PetscFunctionReturn(PETSC_SUCCESS);
86 }
87
main(int argc,char ** argv)88 int main(int argc, char **argv)
89 {
90 AppCtx user; /* user-defined work context */
91 PetscRandom rctx;
92 Vec solution;
93 Mat system;
94 MPI_Comm comm;
95 PetscInt i;
96
97 PetscFunctionBeginUser;
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) PetscCall(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 !complex
142
143 test:
144 args: -debug -fields v1,v2,v3
145
146 TEST*/
147