1 #define PETSCDM_DLL
2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
3
4 #if defined(PETSC_HAVE_EXODUSII)
5 #include <netcdf.h>
6 #include <exodusII.h>
7 #endif
8
9 #include <petsc/private/viewerimpl.h>
10 #include <petsc/private/viewerexodusiiimpl.h>
11 #if defined(PETSC_HAVE_EXODUSII)
12 /*
13 PETSC_VIEWER_EXODUSII_ - Creates an ExodusII PetscViewer shared by all processors in a communicator.
14
15 Collective
16
17 Input Parameter:
18 . comm - the MPI communicator to share the ExodusII PetscViewer
19
20 Level: intermediate
21
22 Notes:
23 misses Fortran bindings
24
25 Notes:
26 Unlike almost all other PETSc routines, PETSC_VIEWER_EXODUSII_ does not return
27 an error code. The GLVIS PetscViewer is usually used in the form
28 $ XXXView(XXX object, PETSC_VIEWER_EXODUSII_(comm));
29
30 .seealso: PetscViewerExodusIIOpen(), PetscViewerType, PetscViewerCreate(), PetscViewerDestroy()
31 */
PETSC_VIEWER_EXODUSII_(MPI_Comm comm)32 PetscViewer PETSC_VIEWER_EXODUSII_(MPI_Comm comm)
33 {
34 PetscViewer viewer;
35 PetscErrorCode ierr;
36
37 PetscFunctionBegin;
38 ierr = PetscViewerExodusIIOpen(comm, "mesh.exo", FILE_MODE_WRITE, &viewer);
39 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);}
40 ierr = PetscObjectRegisterDestroy((PetscObject) viewer);
41 if (ierr) {PetscError(PETSC_COMM_SELF,__LINE__,"PETSC_VIEWER_EXODUSII_",__FILE__,PETSC_ERR_PLIB,PETSC_ERROR_INITIAL," ");PetscFunctionReturn(0);}
42 PetscFunctionReturn(viewer);
43 }
44
PetscViewerView_ExodusII(PetscViewer v,PetscViewer viewer)45 static PetscErrorCode PetscViewerView_ExodusII(PetscViewer v, PetscViewer viewer)
46 {
47 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
48 PetscErrorCode ierr;
49
50 PetscFunctionBegin;
51 if (exo->filename) {ierr = PetscViewerASCIIPrintf(viewer, "Filename: %s\n", exo->filename);CHKERRQ(ierr);}
52 PetscFunctionReturn(0);
53 }
54
PetscViewerSetFromOptions_ExodusII(PetscOptionItems * PetscOptionsObject,PetscViewer v)55 static PetscErrorCode PetscViewerSetFromOptions_ExodusII(PetscOptionItems *PetscOptionsObject, PetscViewer v)
56 {
57 PetscErrorCode ierr;
58
59 PetscFunctionBegin;
60 ierr = PetscOptionsHead(PetscOptionsObject, "ExodusII PetscViewer Options");CHKERRQ(ierr);
61 ierr = PetscOptionsTail();CHKERRQ(ierr);
62 PetscFunctionReturn(0);
63 }
64
PetscViewerSetUp_ExodusII(PetscViewer viewer)65 static PetscErrorCode PetscViewerSetUp_ExodusII(PetscViewer viewer)
66 {
67 PetscFunctionBegin;
68 PetscFunctionReturn(0);
69 }
70
PetscViewerDestroy_ExodusII(PetscViewer viewer)71 static PetscErrorCode PetscViewerDestroy_ExodusII(PetscViewer viewer)
72 {
73 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
74 PetscErrorCode ierr;
75
76 PetscFunctionBegin;
77 if (exo->exoid >= 0) {ierr = ex_close(exo->exoid);CHKERRQ(ierr);}
78 ierr = PetscFree(exo);CHKERRQ(ierr);
79 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetName_C",NULL);CHKERRQ(ierr);
80 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileGetName_C",NULL);CHKERRQ(ierr);
81 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerFileSetMode_C",NULL);CHKERRQ(ierr);
82 ierr = PetscObjectComposeFunction((PetscObject)viewer,"PetscViewerExodusIIGetId",NULL);CHKERRQ(ierr);
83 PetscFunctionReturn(0);
84 }
85
PetscViewerFileSetName_ExodusII(PetscViewer viewer,const char name[])86 static PetscErrorCode PetscViewerFileSetName_ExodusII(PetscViewer viewer, const char name[])
87 {
88 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
89 PetscMPIInt rank;
90 int CPU_word_size, IO_word_size, EXO_mode;
91 PetscErrorCode ierr;
92
93
94 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) viewer), &rank);CHKERRQ(ierr);
95 CPU_word_size = sizeof(PetscReal);
96 IO_word_size = sizeof(PetscReal);
97
98 PetscFunctionBegin;
99 if (exo->exoid >= 0) {ex_close(exo->exoid); exo->exoid = -1;}
100 if (exo->filename) {ierr = PetscFree(exo->filename);CHKERRQ(ierr);}
101 ierr = PetscStrallocpy(name, &exo->filename);CHKERRQ(ierr);
102 /* Create or open the file collectively */
103 switch (exo->btype) {
104 case FILE_MODE_READ:
105 EXO_mode = EX_CLOBBER;
106 break;
107 case FILE_MODE_APPEND:
108 EXO_mode = EX_CLOBBER;
109 break;
110 case FILE_MODE_WRITE:
111 EXO_mode = EX_CLOBBER;
112 break;
113 default:
114 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER, "Must call PetscViewerFileSetMode() before PetscViewerFileSetName()");
115 }
116 #if defined(PETSC_USE_64BIT_INDICES)
117 EXO_mode += EX_ALL_INT64_API;
118 #endif
119 exo->exoid = ex_create(name, EXO_mode, &CPU_word_size, &IO_word_size);
120 if (exo->exoid < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_create failed for %s", name);
121 PetscFunctionReturn(0);
122 }
123
PetscViewerFileGetName_ExodusII(PetscViewer viewer,const char ** name)124 static PetscErrorCode PetscViewerFileGetName_ExodusII(PetscViewer viewer, const char **name)
125 {
126 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
127
128 PetscFunctionBegin;
129 *name = exo->filename;
130 PetscFunctionReturn(0);
131 }
132
PetscViewerFileSetMode_ExodusII(PetscViewer viewer,PetscFileMode type)133 static PetscErrorCode PetscViewerFileSetMode_ExodusII(PetscViewer viewer, PetscFileMode type)
134 {
135 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
136
137 PetscFunctionBegin;
138 exo->btype = type;
139 PetscFunctionReturn(0);
140 }
141
PetscViewerFileGetMode_ExodusII(PetscViewer viewer,PetscFileMode * type)142 static PetscErrorCode PetscViewerFileGetMode_ExodusII(PetscViewer viewer, PetscFileMode *type)
143 {
144 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
145
146 PetscFunctionBegin;
147 *type = exo->btype;
148 PetscFunctionReturn(0);
149 }
150
PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer,int * exoid)151 static PetscErrorCode PetscViewerExodusIIGetId_ExodusII(PetscViewer viewer, int *exoid)
152 {
153 PetscViewer_ExodusII *exo = (PetscViewer_ExodusII *) viewer->data;
154
155 PetscFunctionBegin;
156 *exoid = exo->exoid;
157 PetscFunctionReturn(0);
158 }
159
160 /*@C
161 PetscViewerExodusIIOpen - Opens a file for ExodusII input/output.
162
163 Collective
164
165 Input Parameters:
166 + comm - MPI communicator
167 . name - name of file
168 - type - type of file
169 $ FILE_MODE_WRITE - create new file for binary output
170 $ FILE_MODE_READ - open existing file for binary input
171 $ FILE_MODE_APPEND - open existing file for binary output
172
173 Output Parameter:
174 . exo - PetscViewer for Exodus II input/output to use with the specified file
175
176 Level: beginner
177
178 Note:
179 This PetscViewer should be destroyed with PetscViewerDestroy().
180
181
182 .seealso: PetscViewerPushFormat(), PetscViewerDestroy(),
183 DMLoad(), PetscFileMode, PetscViewer, PetscViewerSetType(), PetscViewerFileSetMode(), PetscViewerFileSetName()
184 @*/
PetscViewerExodusIIOpen(MPI_Comm comm,const char name[],PetscFileMode type,PetscViewer * exo)185 PetscErrorCode PetscViewerExodusIIOpen(MPI_Comm comm, const char name[], PetscFileMode type, PetscViewer *exo)
186 {
187 PetscErrorCode ierr;
188
189 PetscFunctionBegin;
190 ierr = PetscViewerCreate(comm, exo);CHKERRQ(ierr);
191 ierr = PetscViewerSetType(*exo, PETSCVIEWEREXODUSII);CHKERRQ(ierr);
192 ierr = PetscViewerFileSetMode(*exo, type);CHKERRQ(ierr);
193 ierr = PetscViewerFileSetName(*exo, name);CHKERRQ(ierr);
194 ierr = PetscViewerSetFromOptions(*exo);CHKERRQ(ierr);
195 PetscFunctionReturn(0);
196 }
197
198 /*MC
199 PETSCVIEWEREXODUSII - A viewer that writes to an Exodus II file
200
201
202 .seealso: PetscViewerExodusIIOpen(), PetscViewerCreate(), PETSCVIEWERBINARY, PETSCVIEWERHDF5, DMView(),
203 PetscViewerFileSetName(), PetscViewerFileSetMode(), PetscViewerFormat, PetscViewerType, PetscViewerSetType()
204
205 Level: beginner
206 M*/
207
PetscViewerCreate_ExodusII(PetscViewer v)208 PETSC_EXTERN PetscErrorCode PetscViewerCreate_ExodusII(PetscViewer v)
209 {
210 PetscViewer_ExodusII *exo;
211 PetscErrorCode ierr;
212
213 PetscFunctionBegin;
214 ierr = PetscNewLog(v,&exo);CHKERRQ(ierr);
215
216 v->data = (void*) exo;
217 v->ops->destroy = PetscViewerDestroy_ExodusII;
218 v->ops->setfromoptions = PetscViewerSetFromOptions_ExodusII;
219 v->ops->setup = PetscViewerSetUp_ExodusII;
220 v->ops->view = PetscViewerView_ExodusII;
221 v->ops->flush = 0;
222 exo->btype = (PetscFileMode) -1;
223 exo->filename = 0;
224 exo->exoid = -1;
225
226 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetName_C",PetscViewerFileSetName_ExodusII);CHKERRQ(ierr);
227 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetName_C",PetscViewerFileGetName_ExodusII);CHKERRQ(ierr);
228 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileSetMode_C",PetscViewerFileSetMode_ExodusII);CHKERRQ(ierr);
229 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerFileGetMode_C",PetscViewerFileGetMode_ExodusII);CHKERRQ(ierr);
230 ierr = PetscObjectComposeFunction((PetscObject)v,"PetscViewerExodusIIGetId_C",PetscViewerExodusIIGetId_ExodusII);CHKERRQ(ierr);
231 PetscFunctionReturn(0);
232 }
233
234 /*
235 EXOGetVarIndex - Locate a result in an exodus file based on its name
236
237 Not collective
238
239 Input Parameters:
240 + exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
241 . obj_type - the type of entity for instance EX_NODAL, EX_ELEM_BLOCK
242 - name - the name of the result
243
244 Output Parameters:
245 . varIndex - the location in the exodus file of the result
246
247 Notes:
248 The exodus variable index is obtained by comparing name and the
249 names of zonal variables declared in the exodus file. For instance if name is "V"
250 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
251 amongst all variables of type obj_type.
252
253 Level: beginner
254
255 .seealso: EXOGetVarIndex(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadNodal_PlexEXO(),VecLoadZonal_PlexEXO()
256 */
EXOGetVarIndex_Private(int exoid,ex_entity_type obj_type,const char name[],int * varIndex)257 static PetscErrorCode EXOGetVarIndex_Private(int exoid, ex_entity_type obj_type, const char name[], int *varIndex)
258 {
259 int num_vars, i, j;
260 char ext_name[MAX_STR_LENGTH+1], var_name[MAX_STR_LENGTH+1];
261 const int num_suffix = 5;
262 char *suffix[5];
263 PetscBool flg;
264 PetscErrorCode ierr;
265
266 PetscFunctionBegin;
267 suffix[0] = (char *) "";
268 suffix[1] = (char *) "_X";
269 suffix[2] = (char *) "_XX";
270 suffix[3] = (char *) "_1";
271 suffix[4] = (char *) "_11";
272
273 *varIndex = 0;
274 PetscStackCallStandard(ex_get_variable_param,(exoid, obj_type, &num_vars));
275 for (i = 0; i < num_vars; ++i) {
276 PetscStackCallStandard(ex_get_variable_name,(exoid, obj_type, i+1, var_name));
277 for (j = 0; j < num_suffix; ++j){
278 ierr = PetscStrncpy(ext_name, name, MAX_STR_LENGTH);CHKERRQ(ierr);
279 ierr = PetscStrlcat(ext_name, suffix[j], MAX_STR_LENGTH);CHKERRQ(ierr);
280 ierr = PetscStrcasecmp(ext_name, var_name, &flg);
281 if (flg) {
282 *varIndex = i+1;
283 PetscFunctionReturn(0);
284 }
285 }
286 }
287 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to locate variable %s in ExodusII file.", name);
288 PetscFunctionReturn(-1);
289 }
290
291 /*
292 DMPlexView_ExodusII_Internal - Write a DM to disk in exodus format
293
294 Collective on dm
295
296 Input Parameters:
297 + dm - The dm to be written
298 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
299 - degree - the degree of the interpolation space
300
301 Notes:
302 Not all DM can be written to disk this way. For instance, exodus assume that element blocks (mapped to "Cell sets" labels)
303 consists of sequentially numbered cells. If this is not the case, the exodus file will be corrupted.
304
305 If the dm has been distributed and exoid points to different files on each MPI rank, only the "local" part of the DM
306 (including "ghost" cells and vertices) will be written. If all exoid points to the same file, the resulting file will
307 probably be corrupted.
308
309 DMPlex only represents geometry while most post-processing software expect that a mesh also provides information
310 on the discretization space. This function assumes that the file represents Lagrange finite elements of order 1 or 2.
311 It should be extended to use PetscFE objects.
312
313 This function will only handle TRI, TET, QUAD and HEX cells.
314 Level: beginner
315
316 .seealso: EXOGetVarIndex_Private(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadNodal_PlexEXO(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
317 */
DMPlexView_ExodusII_Internal(DM dm,int exoid,PetscInt degree)318 PetscErrorCode DMPlexView_ExodusII_Internal(DM dm, int exoid, PetscInt degree)
319 {
320 enum ElemType {TRI, QUAD, TET, HEX};
321 MPI_Comm comm;
322 /* Connectivity Variables */
323 PetscInt cellsNotInConnectivity;
324 /* Cell Sets */
325 DMLabel csLabel;
326 IS csIS;
327 const PetscInt *csIdx;
328 PetscInt num_cs, cs;
329 enum ElemType *type;
330 PetscBool hasLabel;
331 /* Coordinate Variables */
332 DM cdm;
333 PetscSection section;
334 Vec coord;
335 PetscInt **nodes;
336 PetscInt depth, d, dim, skipCells = 0;
337 PetscInt pStart, pEnd, p, cStart, cEnd, numCells, vStart, vEnd, numVertices, eStart, eEnd, numEdges, fStart, fEnd, numFaces, numNodes;
338 PetscInt num_vs, num_fs;
339 PetscMPIInt rank, size;
340 const char *dmName;
341 PetscInt nodesTriP1[4] = {3,0,0,0};
342 PetscInt nodesTriP2[4] = {3,3,0,0};
343 PetscInt nodesQuadP1[4] = {4,0,0,0};
344 PetscInt nodesQuadP2[4] = {4,4,0,1};
345 PetscInt nodesTetP1[4] = {4,0,0,0};
346 PetscInt nodesTetP2[4] = {4,6,0,0};
347 PetscInt nodesHexP1[4] = {8,0,0,0};
348 PetscInt nodesHexP2[4] = {8,12,6,1};
349 PetscErrorCode ierr;
350
351 PetscFunctionBegin;
352 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr);
353 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
354 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
355 /* --- Get DM info --- */
356 ierr = PetscObjectGetName((PetscObject) dm, &dmName);CHKERRQ(ierr);
357 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
358 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
359 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
360 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
361 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
362 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
363 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
364 numCells = cEnd - cStart;
365 numEdges = eEnd - eStart;
366 numVertices = vEnd - vStart;
367 if (depth == 3) {numFaces = fEnd - fStart;}
368 else {numFaces = 0;}
369 ierr = DMGetLabelSize(dm, "Cell Sets", &num_cs);CHKERRQ(ierr);
370 ierr = DMGetLabelSize(dm, "Vertex Sets", &num_vs);CHKERRQ(ierr);
371 ierr = DMGetLabelSize(dm, "Face Sets", &num_fs);CHKERRQ(ierr);
372 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr);
373 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
374 if (num_cs > 0) {
375 ierr = DMGetLabel(dm, "Cell Sets", &csLabel);CHKERRQ(ierr);
376 ierr = DMLabelGetValueIS(csLabel, &csIS);CHKERRQ(ierr);
377 ierr = ISGetIndices(csIS, &csIdx);CHKERRQ(ierr);
378 }
379 ierr = PetscMalloc1(num_cs, &nodes);CHKERRQ(ierr);
380 /* Set element type for each block and compute total number of nodes */
381 ierr = PetscMalloc1(num_cs, &type);CHKERRQ(ierr);
382 numNodes = numVertices;
383 if (degree == 2) {numNodes += numEdges;}
384 cellsNotInConnectivity = numCells;
385 for (cs = 0; cs < num_cs; ++cs) {
386 IS stratumIS;
387 const PetscInt *cells;
388 PetscScalar *xyz = NULL;
389 PetscInt csSize, closureSize;
390
391 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr);
392 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr);
393 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr);
394 ierr = DMPlexVecGetClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr);
395 switch (dim) {
396 case 2:
397 if (closureSize == 3*dim) {type[cs] = TRI;}
398 else if (closureSize == 4*dim) {type[cs] = QUAD;}
399 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
400 break;
401 case 3:
402 if (closureSize == 4*dim) {type[cs] = TET;}
403 else if (closureSize == 8*dim) {type[cs] = HEX;}
404 else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of vertices %D in dimension %D has no ExodusII type", closureSize/dim, dim);
405 break;
406 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D not handled by ExodusII viewer", dim);
407 }
408 if ((degree == 2) && (type[cs] == QUAD)) {numNodes += csSize;}
409 if ((degree == 2) && (type[cs] == HEX)) {numNodes += csSize; numNodes += numFaces;}
410 ierr = DMPlexVecRestoreClosure(cdm, NULL, coord, cells[0], &closureSize, &xyz);CHKERRQ(ierr);
411 /* Set nodes and Element type */
412 if (type[cs] == TRI) {
413 if (degree == 1) nodes[cs] = nodesTriP1;
414 else if (degree == 2) nodes[cs] = nodesTriP2;
415 } else if (type[cs] == QUAD) {
416 if (degree == 1) nodes[cs] = nodesQuadP1;
417 else if (degree == 2) nodes[cs] = nodesQuadP2;
418 } else if (type[cs] == TET) {
419 if (degree == 1) nodes[cs] = nodesTetP1;
420 else if (degree == 2) nodes[cs] = nodesTetP2;
421 } else if (type[cs] == HEX) {
422 if (degree == 1) nodes[cs] = nodesHexP1;
423 else if (degree == 2) nodes[cs] = nodesHexP2;
424 }
425 /* Compute the number of cells not in the connectivity table */
426 cellsNotInConnectivity -= nodes[cs][3]*csSize;
427
428 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr);
429 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
430 }
431 if (num_cs > 0) {PetscStackCallStandard(ex_put_init,(exoid, dmName, dim, numNodes, numCells, num_cs, num_vs, num_fs));}
432 /* --- Connectivity --- */
433 for (cs = 0; cs < num_cs; ++cs) {
434 IS stratumIS;
435 const PetscInt *cells;
436 PetscInt *connect, off = 0;
437 PetscInt edgesInClosure = 0, facesInClosure = 0, verticesInClosure = 0;
438 PetscInt csSize, c, connectSize, closureSize;
439 char *elem_type = NULL;
440 char elem_type_tri3[] = "TRI3", elem_type_quad4[] = "QUAD4";
441 char elem_type_tri6[] = "TRI6", elem_type_quad9[] = "QUAD9";
442 char elem_type_tet4[] = "TET4", elem_type_hex8[] = "HEX8";
443 char elem_type_tet10[] = "TET10", elem_type_hex27[] = "HEX27";
444
445 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr);
446 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr);
447 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr);
448 /* Set Element type */
449 if (type[cs] == TRI) {
450 if (degree == 1) elem_type = elem_type_tri3;
451 else if (degree == 2) elem_type = elem_type_tri6;
452 } else if (type[cs] == QUAD) {
453 if (degree == 1) elem_type = elem_type_quad4;
454 else if (degree == 2) elem_type = elem_type_quad9;
455 } else if (type[cs] == TET) {
456 if (degree == 1) elem_type = elem_type_tet4;
457 else if (degree == 2) elem_type = elem_type_tet10;
458 } else if (type[cs] == HEX) {
459 if (degree == 1) elem_type = elem_type_hex8;
460 else if (degree == 2) elem_type = elem_type_hex27;
461 }
462 connectSize = nodes[cs][0] + nodes[cs][1] + nodes[cs][2] + nodes[cs][3];
463 ierr = PetscMalloc1(PetscMax(27,connectSize)*csSize, &connect);CHKERRQ(ierr);
464 PetscStackCallStandard(ex_put_block,(exoid, EX_ELEM_BLOCK, csIdx[cs], elem_type, csSize, connectSize, 0, 0, 1));
465 /* Find number of vertices, edges, and faces in the closure */
466 verticesInClosure = nodes[cs][0];
467 if (depth > 1) {
468 if (dim == 2) {
469 ierr = DMPlexGetConeSize(dm, cells[0], &edgesInClosure);CHKERRQ(ierr);
470 } else if (dim == 3) {
471 PetscInt *closure = NULL;
472
473 ierr = DMPlexGetConeSize(dm, cells[0], &facesInClosure);CHKERRQ(ierr);
474 ierr = DMPlexGetTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
475 edgesInClosure = closureSize - facesInClosure - 1 - verticesInClosure;
476 ierr = DMPlexRestoreTransitiveClosure(dm, cells[0], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
477 }
478 }
479 /* Get connectivity for each cell */
480 for (c = 0; c < csSize; ++c) {
481 PetscInt *closure = NULL;
482 PetscInt temp, i;
483
484 ierr = DMPlexGetTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
485 for (i = 0; i < connectSize; ++i) {
486 if (i < nodes[cs][0]) {/* Vertices */
487 connect[i+off] = closure[(i+edgesInClosure+facesInClosure+1)*2] + 1;
488 connect[i+off] -= cellsNotInConnectivity;
489 } else if (i < nodes[cs][0]+nodes[cs][1]) { /* Edges */
490 connect[i+off] = closure[(i-verticesInClosure+facesInClosure+1)*2] + 1;
491 if (nodes[cs][2] == 0) connect[i+off] -= numFaces;
492 connect[i+off] -= cellsNotInConnectivity;
493 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]){ /* Cells */
494 connect[i+off] = closure[0] + 1;
495 connect[i+off] -= skipCells;
496 } else if (i < nodes[cs][0]+nodes[cs][1]+nodes[cs][3]+nodes[cs][2]){ /* Faces */
497 connect[i+off] = closure[(i-edgesInClosure-verticesInClosure)*2] + 1;
498 connect[i+off] -= cellsNotInConnectivity;
499 } else {
500 connect[i+off] = -1;
501 }
502 }
503 /* Tetrahedra are inverted */
504 if (type[cs] == TET) {
505 temp = connect[0+off]; connect[0+off] = connect[1+off]; connect[1+off] = temp;
506 if (degree == 2) {
507 temp = connect[5+off]; connect[5+off] = connect[6+off]; connect[6+off] = temp;
508 temp = connect[7+off]; connect[7+off] = connect[8+off]; connect[8+off] = temp;
509 }
510 }
511 /* Hexahedra are inverted */
512 if (type[cs] == HEX) {
513 temp = connect[1+off]; connect[1+off] = connect[3+off]; connect[3+off] = temp;
514 if (degree == 2) {
515 temp = connect[8+off]; connect[8+off] = connect[11+off]; connect[11+off] = temp;
516 temp = connect[9+off]; connect[9+off] = connect[10+off]; connect[10+off] = temp;
517 temp = connect[16+off]; connect[16+off] = connect[17+off]; connect[17+off] = temp;
518 temp = connect[18+off]; connect[18+off] = connect[19+off]; connect[19+off] = temp;
519
520 temp = connect[12+off]; connect[12+off] = connect[16+off]; connect[16+off] = temp;
521 temp = connect[13+off]; connect[13+off] = connect[17+off]; connect[17+off] = temp;
522 temp = connect[14+off]; connect[14+off] = connect[18+off]; connect[18+off] = temp;
523 temp = connect[15+off]; connect[15+off] = connect[19+off]; connect[19+off] = temp;
524
525 temp = connect[23+off]; connect[23+off] = connect[26+off]; connect[26+off] = temp;
526 temp = connect[24+off]; connect[24+off] = connect[25+off]; connect[25+off] = temp;
527 temp = connect[25+off]; connect[25+off] = connect[26+off]; connect[26+off] = temp;
528 }
529 }
530 off += connectSize;
531 ierr = DMPlexRestoreTransitiveClosure(dm, cells[c], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
532 }
533 PetscStackCallStandard(ex_put_conn,(exoid, EX_ELEM_BLOCK, csIdx[cs], connect, 0, 0));
534 skipCells += (nodes[cs][3] == 0)*csSize;
535 ierr = PetscFree(connect);CHKERRQ(ierr);
536 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr);
537 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
538 }
539 ierr = PetscFree(type);CHKERRQ(ierr);
540 /* --- Coordinates --- */
541 ierr = PetscSectionCreate(comm, §ion);CHKERRQ(ierr);
542 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr);
543 if (num_cs) {
544 for (d = 0; d < depth; ++d) {
545 ierr = DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);CHKERRQ(ierr);
546 for (p = pStart; p < pEnd; ++p) {
547 ierr = PetscSectionSetDof(section, p, nodes[0][d] > 0);CHKERRQ(ierr);
548 }
549 }
550 }
551 for (cs = 0; cs < num_cs; ++cs) {
552 IS stratumIS;
553 const PetscInt *cells;
554 PetscInt csSize, c;
555
556 ierr = DMLabelGetStratumIS(csLabel, csIdx[cs], &stratumIS);CHKERRQ(ierr);
557 ierr = ISGetIndices(stratumIS, &cells);CHKERRQ(ierr);
558 ierr = ISGetSize(stratumIS, &csSize);CHKERRQ(ierr);
559 for (c = 0; c < csSize; ++c) {
560 ierr = PetscSectionSetDof(section, cells[c], nodes[cs][3] > 0);CHKERRQ(ierr);
561 }
562 ierr = ISRestoreIndices(stratumIS, &cells);CHKERRQ(ierr);
563 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
564 }
565 if (num_cs > 0) {
566 ierr = ISRestoreIndices(csIS, &csIdx);CHKERRQ(ierr);
567 ierr = ISDestroy(&csIS);CHKERRQ(ierr);
568 }
569 ierr = PetscFree(nodes);CHKERRQ(ierr);
570 ierr = PetscSectionSetUp(section);CHKERRQ(ierr);
571 if (numNodes > 0) {
572 const char *coordNames[3] = {"x", "y", "z"};
573 PetscScalar *coords, *closure;
574 PetscReal *cval;
575 PetscInt hasDof, n = 0;
576
577 /* There can't be more than 24 values in the closure of a point for the coord section */
578 ierr = PetscCalloc3(numNodes*3, &coords, dim, &cval, 24, &closure);CHKERRQ(ierr);
579 ierr = DMGetCoordinatesLocal(dm, &coord);CHKERRQ(ierr);
580 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
581 for (p = pStart; p < pEnd; ++p) {
582 ierr = PetscSectionGetDof(section, p, &hasDof);
583 if (hasDof) {
584 PetscInt closureSize = 24, j;
585
586 ierr = DMPlexVecGetClosure(cdm, NULL, coord, p, &closureSize, &closure);CHKERRQ(ierr);
587 for (d = 0; d < dim; ++d) {
588 cval[d] = 0.0;
589 for (j = 0; j < closureSize/dim; j++) cval[d] += closure[j*dim+d];
590 coords[d*numNodes+n] = cval[d] * dim / closureSize;
591 }
592 ++n;
593 }
594 }
595 PetscStackCallStandard(ex_put_coord,(exoid, &coords[0*numNodes], &coords[1*numNodes], &coords[2*numNodes]));
596 ierr = PetscFree3(coords, cval, closure);CHKERRQ(ierr);
597 PetscStackCallStandard(ex_put_coord_names,(exoid, (char **) coordNames));
598 }
599 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr);
600 /* --- Node Sets/Vertex Sets --- */
601 DMHasLabel(dm, "Vertex Sets", &hasLabel);
602 if (hasLabel) {
603 PetscInt i, vs, vsSize;
604 const PetscInt *vsIdx, *vertices;
605 PetscInt *nodeList;
606 IS vsIS, stratumIS;
607 DMLabel vsLabel;
608 ierr = DMGetLabel(dm, "Vertex Sets", &vsLabel);CHKERRQ(ierr);
609 ierr = DMLabelGetValueIS(vsLabel, &vsIS);CHKERRQ(ierr);
610 ierr = ISGetIndices(vsIS, &vsIdx);CHKERRQ(ierr);
611 for (vs=0; vs<num_vs; ++vs) {
612 ierr = DMLabelGetStratumIS(vsLabel, vsIdx[vs], &stratumIS);CHKERRQ(ierr);
613 ierr = ISGetIndices(stratumIS, &vertices);CHKERRQ(ierr);
614 ierr = ISGetSize(stratumIS, &vsSize);CHKERRQ(ierr);
615 ierr = PetscMalloc1(vsSize, &nodeList);
616 for (i=0; i<vsSize; ++i) {
617 nodeList[i] = vertices[i] - skipCells + 1;
618 }
619 PetscStackCallStandard(ex_put_set_param,(exoid, EX_NODE_SET, vsIdx[vs], vsSize, 0));
620 PetscStackCallStandard(ex_put_set,(exoid, EX_NODE_SET, vsIdx[vs], nodeList, NULL));
621 ierr = ISRestoreIndices(stratumIS, &vertices);CHKERRQ(ierr);
622 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
623 ierr = PetscFree(nodeList);CHKERRQ(ierr);
624 }
625 ierr = ISRestoreIndices(vsIS, &vsIdx);CHKERRQ(ierr);
626 ierr = ISDestroy(&vsIS);CHKERRQ(ierr);
627 }
628 /* --- Side Sets/Face Sets --- */
629 ierr = DMHasLabel(dm, "Face Sets", &hasLabel);CHKERRQ(ierr);
630 if (hasLabel) {
631 PetscInt i, j, fs, fsSize;
632 const PetscInt *fsIdx, *faces;
633 IS fsIS, stratumIS;
634 DMLabel fsLabel;
635 PetscInt numPoints, *points;
636 PetscInt elem_list_size = 0;
637 PetscInt *elem_list, *elem_ind, *side_list;
638
639 ierr = DMGetLabel(dm, "Face Sets", &fsLabel);CHKERRQ(ierr);
640 /* Compute size of Node List and Element List */
641 ierr = DMLabelGetValueIS(fsLabel, &fsIS);CHKERRQ(ierr);
642 ierr = ISGetIndices(fsIS, &fsIdx);CHKERRQ(ierr);
643 for (fs=0; fs<num_fs; ++fs) {
644 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr);
645 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr);
646 elem_list_size += fsSize;
647 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
648 }
649 ierr = PetscMalloc3(num_fs, &elem_ind,
650 elem_list_size, &elem_list,
651 elem_list_size, &side_list);CHKERRQ(ierr);
652 elem_ind[0] = 0;
653 for (fs=0; fs<num_fs; ++fs) {
654 ierr = DMLabelGetStratumIS(fsLabel, fsIdx[fs], &stratumIS);CHKERRQ(ierr);
655 ierr = ISGetIndices(stratumIS, &faces);CHKERRQ(ierr);
656 ierr = ISGetSize(stratumIS, &fsSize);CHKERRQ(ierr);
657 /* Set Parameters */
658 PetscStackCallStandard(ex_put_set_param,(exoid, EX_SIDE_SET, fsIdx[fs], fsSize, 0));
659 /* Indices */
660 if (fs<num_fs-1) {
661 elem_ind[fs+1] = elem_ind[fs] + fsSize;
662 }
663
664 for (i=0; i<fsSize; ++i) {
665 /* Element List */
666 points = NULL;
667 ierr = DMPlexGetTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr);
668 elem_list[elem_ind[fs] + i] = points[2] +1;
669 ierr = DMPlexRestoreTransitiveClosure(dm, faces[i], PETSC_FALSE, &numPoints, &points);CHKERRQ(ierr);
670
671 /* Side List */
672 points = NULL;
673 ierr = DMPlexGetTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
674 for (j=1; j<numPoints; ++j) {
675 if (points[j*2]==faces[i]) {break;}
676 }
677 /* Convert HEX sides */
678 if (numPoints == 27) {
679 if (j == 1) {j = 5;}
680 else if (j == 2) {j = 6;}
681 else if (j == 3) {j = 1;}
682 else if (j == 4) {j = 3;}
683 else if (j == 5) {j = 2;}
684 else if (j == 6) {j = 4;}
685 }
686 /* Convert TET sides */
687 if (numPoints == 15) {
688 --j;
689 if (j == 0) {j = 4;}
690 }
691 side_list[elem_ind[fs] + i] = j;
692 ierr = DMPlexRestoreTransitiveClosure(dm, elem_list[elem_ind[fs] + i]-1, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
693
694 }
695 ierr = ISRestoreIndices(stratumIS, &faces);CHKERRQ(ierr);
696 ierr = ISDestroy(&stratumIS);CHKERRQ(ierr);
697 }
698 ierr = ISRestoreIndices(fsIS, &fsIdx);CHKERRQ(ierr);
699 ierr = ISDestroy(&fsIS);CHKERRQ(ierr);
700
701 /* Put side sets */
702 for (fs=0; fs<num_fs; ++fs) {
703 PetscStackCallStandard(ex_put_set,(exoid, EX_SIDE_SET, fsIdx[fs], &elem_list[elem_ind[fs]], &side_list[elem_ind[fs]]));
704 }
705 ierr = PetscFree3(elem_ind,elem_list,side_list);CHKERRQ(ierr);
706 }
707 PetscFunctionReturn(0);
708 }
709
710 /*
711 VecViewPlex_ExodusII_Nodal_Internal - Write a Vec corresponding to a nodal field to an exodus file
712
713 Collective on v
714
715 Input Parameters:
716 + v - The vector to be written
717 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
718 - step - the time step to write at (exodus steps are numbered starting from 1)
719
720 Notes:
721 The exodus result nodal variable index is obtained by comparing the Vec name and the
722 names of nodal variables declared in the exodus file. For instance for a Vec named "V"
723 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
724 amongst all nodal variables.
725
726 Level: beginner
727
728 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecLoadNodal_PlexEXO(),VecViewZonal_PlexEXO(),VecLoadZonal_PlexEXO()
729 @*/
VecViewPlex_ExodusII_Nodal_Internal(Vec v,int exoid,int step)730 PetscErrorCode VecViewPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
731 {
732 MPI_Comm comm;
733 PetscMPIInt size;
734 DM dm;
735 Vec vNatural, vComp;
736 const PetscScalar *varray;
737 const char *vecname;
738 PetscInt xs, xe, bs;
739 PetscBool useNatural;
740 int offset;
741 PetscErrorCode ierr;
742
743 PetscFunctionBegin;
744 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr);
745 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
746 ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
747 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr);
748 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
749 if (useNatural) {
750 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr);
751 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr);
752 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr);
753 } else {
754 vNatural = v;
755 }
756 /* Get the location of the variable in the exodus file */
757 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr);
758 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr);
759 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file\n", vecname);
760 /* Write local chunk of the result in the exodus file
761 exodus stores each component of a vector-valued field as a separate variable.
762 We assume that they are stored sequentially */
763 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr);
764 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr);
765 if (bs == 1) {
766 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr);
767 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
768 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr);
769 } else {
770 IS compIS;
771 PetscInt c;
772
773 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);
774 for (c = 0; c < bs; ++c) {
775 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr);
776 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
777 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr);
778 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
779 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr);
780 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
781 }
782 ierr = ISDestroy(&compIS);CHKERRQ(ierr);
783 }
784 if (useNatural) {ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);}
785 PetscFunctionReturn(0);
786 }
787
788 /*
789 VecLoadPlex_ExodusII_Nodal_Internal - Read a Vec corresponding to a nodal field from an exodus file
790
791 Collective on v
792
793 Input Parameters:
794 + v - The vector to be written
795 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
796 - step - the time step to read at (exodus steps are numbered starting from 1)
797
798 Notes:
799 The exodus result nodal variable index is obtained by comparing the Vec name and the
800 names of nodal variables declared in the exodus file. For instance for a Vec named "V"
801 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
802 amongst all nodal variables.
803
804 Level: beginner
805
806 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecViewZonal_PlexEXO(), VecLoadZonal_PlexEXO()
807 */
VecLoadPlex_ExodusII_Nodal_Internal(Vec v,int exoid,int step)808 PetscErrorCode VecLoadPlex_ExodusII_Nodal_Internal(Vec v, int exoid, int step)
809 {
810 MPI_Comm comm;
811 PetscMPIInt size;
812 DM dm;
813 Vec vNatural, vComp;
814 PetscScalar *varray;
815 const char *vecname;
816 PetscInt xs, xe, bs;
817 PetscBool useNatural;
818 int offset;
819 PetscErrorCode ierr;
820
821 PetscFunctionBegin;
822 ierr = PetscObjectGetComm((PetscObject) v, &comm);CHKERRQ(ierr);
823 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
824 ierr = VecGetDM(v,&dm);CHKERRQ(ierr);
825 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr);
826 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
827 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);}
828 else {vNatural = v;}
829 /* Get the location of the variable in the exodus file */
830 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr);
831 ierr = EXOGetVarIndex_Private(exoid, EX_NODAL, vecname, &offset);CHKERRQ(ierr);
832 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate nodal variable %s in exodus file", vecname);
833 /* Read local chunk from the file */
834 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr);
835 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr);
836 if (bs == 1) {
837 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr);
838 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset, 1, xs+1, xe-xs, varray));
839 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr);
840 } else {
841 IS compIS;
842 PetscInt c;
843
844 ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);
845 for (c = 0; c < bs; ++c) {
846 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr);
847 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
848 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr);
849 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_NODAL, offset+c, 1, xs/bs+1, (xe-xs)/bs, varray));
850 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr);
851 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
852 }
853 ierr = ISDestroy(&compIS);CHKERRQ(ierr);
854 }
855 if (useNatural) {
856 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr);
857 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr);
858 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);
859 }
860 PetscFunctionReturn(0);
861 }
862
863 /*
864 VecViewPlex_ExodusII_Zonal_Internal - Write a Vec corresponding to a zonal (cell based) field to an exodus file
865
866 Collective on v
867
868 Input Parameters:
869 + v - The vector to be written
870 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
871 - step - the time step to write at (exodus steps are numbered starting from 1)
872
873 Notes:
874 The exodus result zonal variable index is obtained by comparing the Vec name and the
875 names of zonal variables declared in the exodus file. For instance for a Vec named "V"
876 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
877 amongst all zonal variables.
878
879 Level: beginner
880
881 .seealso: EXOGetVarIndex_Private(),DMPlexView_ExodusII_Internal(),VecViewPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Nodal_Internal(),VecLoadPlex_ExodusII_Zonal_Internal()
882 */
VecViewPlex_ExodusII_Zonal_Internal(Vec v,int exoid,int step)883 PetscErrorCode VecViewPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
884 {
885 MPI_Comm comm;
886 PetscMPIInt size;
887 DM dm;
888 Vec vNatural, vComp;
889 const PetscScalar *varray;
890 const char *vecname;
891 PetscInt xs, xe, bs;
892 PetscBool useNatural;
893 int offset;
894 IS compIS;
895 PetscInt *csSize, *csID;
896 PetscInt numCS, set, csxs = 0;
897 PetscErrorCode ierr;
898
899 PetscFunctionBegin;
900 ierr = PetscObjectGetComm((PetscObject)v, &comm);CHKERRQ(ierr);
901 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
902 ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
903 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr);
904 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
905 if (useNatural) {
906 ierr = DMGetGlobalVector(dm, &vNatural);CHKERRQ(ierr);
907 ierr = DMPlexGlobalToNaturalBegin(dm, v, vNatural);CHKERRQ(ierr);
908 ierr = DMPlexGlobalToNaturalEnd(dm, v, vNatural);CHKERRQ(ierr);
909 } else {
910 vNatural = v;
911 }
912 /* Get the location of the variable in the exodus file */
913 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr);
914 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr);
915 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
916 /* Write local chunk of the result in the exodus file
917 exodus stores each component of a vector-valued field as a separate variable.
918 We assume that they are stored sequentially
919 Zonal variables are accessed one element block at a time, so we loop through the cell sets,
920 but once the vector has been reordered to natural size, we cannot use the label informations
921 to figure out what to save where. */
922 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
923 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr);
924 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
925 for (set = 0; set < numCS; ++set) {
926 ex_block block;
927
928 block.id = csID[set];
929 block.type = EX_ELEM_BLOCK;
930 PetscStackCallStandard(ex_get_block_param,(exoid, &block));
931 csSize[set] = block.num_entry;
932 }
933 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr);
934 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr);
935 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);}
936 for (set = 0; set < numCS; set++) {
937 PetscInt csLocalSize, c;
938
939 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
940 local slice of zonal values: xs/bs,xm/bs-1
941 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
942 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
943 if (bs == 1) {
944 ierr = VecGetArrayRead(vNatural, &varray);CHKERRQ(ierr);
945 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
946 ierr = VecRestoreArrayRead(vNatural, &varray);CHKERRQ(ierr);
947 } else {
948 for (c = 0; c < bs; ++c) {
949 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr);
950 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
951 ierr = VecGetArrayRead(vComp, &varray);CHKERRQ(ierr);
952 PetscStackCallStandard(ex_put_partial_var,(exoid, step, EX_ELEM_BLOCK, offset+c, csID[set], PetscMax(xs/bs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs/bs)]));
953 ierr = VecRestoreArrayRead(vComp, &varray);CHKERRQ(ierr);
954 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
955 }
956 }
957 csxs += csSize[set];
958 }
959 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr);
960 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);}
961 if (useNatural) {ierr = DMRestoreGlobalVector(dm,&vNatural);CHKERRQ(ierr);}
962 PetscFunctionReturn(0);
963 }
964
965 /*
966 VecLoadPlex_ExodusII_Zonal_Internal - Read a Vec corresponding to a zonal (cell based) field from an exodus file
967
968 Collective on v
969
970 Input Parameters:
971 + v - The vector to be written
972 . exoid - the exodus id of a file (obtained from ex_open or ex_create for instance)
973 - step - the time step to read at (exodus steps are numbered starting from 1)
974
975 Notes:
976 The exodus result zonal variable index is obtained by comparing the Vec name and the
977 names of zonal variables declared in the exodus file. For instance for a Vec named "V"
978 the location in the exodus file will be the first match of "V", "V_X", "V_XX", "V_1", or "V_11"
979 amongst all zonal variables.
980
981 Level: beginner
982
983 .seealso: EXOGetVarIndex_Private(), DMPlexView_ExodusII_Internal(), VecViewPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Nodal_Internal(), VecLoadPlex_ExodusII_Zonal_Internal()
984 */
VecLoadPlex_ExodusII_Zonal_Internal(Vec v,int exoid,int step)985 PetscErrorCode VecLoadPlex_ExodusII_Zonal_Internal(Vec v, int exoid, int step)
986 {
987 MPI_Comm comm;
988 PetscMPIInt size;
989 DM dm;
990 Vec vNatural, vComp;
991 PetscScalar *varray;
992 const char *vecname;
993 PetscInt xs, xe, bs;
994 PetscBool useNatural;
995 int offset;
996 IS compIS;
997 PetscInt *csSize, *csID;
998 PetscInt numCS, set, csxs = 0;
999 PetscErrorCode ierr;
1000
1001 PetscFunctionBegin;
1002 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
1003 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
1004 ierr = VecGetDM(v, &dm);CHKERRQ(ierr);
1005 ierr = DMGetUseNatural(dm, &useNatural);CHKERRQ(ierr);
1006 useNatural = useNatural && size > 1 ? PETSC_TRUE : PETSC_FALSE;
1007 if (useNatural) {ierr = DMGetGlobalVector(dm,&vNatural);CHKERRQ(ierr);}
1008 else {vNatural = v;}
1009 /* Get the location of the variable in the exodus file */
1010 ierr = PetscObjectGetName((PetscObject) v, &vecname);CHKERRQ(ierr);
1011 ierr = EXOGetVarIndex_Private(exoid, EX_ELEM_BLOCK, vecname, &offset);CHKERRQ(ierr);
1012 if (!offset) SETERRQ1(comm, PETSC_ERR_FILE_UNEXPECTED, "Cannot locate zonal variable %s in exodus file\n", vecname);
1013 /* Read local chunk of the result in the exodus file
1014 exodus stores each component of a vector-valued field as a separate variable.
1015 We assume that they are stored sequentially
1016 Zonal variables are accessed one element block at a time, so we loop through the cell sets,
1017 but once the vector has been reordered to natural size, we cannot use the label informations
1018 to figure out what to save where. */
1019 numCS = ex_inquire_int(exoid, EX_INQ_ELEM_BLK);
1020 ierr = PetscMalloc2(numCS, &csID, numCS, &csSize);CHKERRQ(ierr);
1021 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, csID));
1022 for (set = 0; set < numCS; ++set) {
1023 ex_block block;
1024
1025 block.id = csID[set];
1026 block.type = EX_ELEM_BLOCK;
1027 PetscStackCallStandard(ex_get_block_param,(exoid, &block));
1028 csSize[set] = block.num_entry;
1029 }
1030 ierr = VecGetOwnershipRange(vNatural, &xs, &xe);CHKERRQ(ierr);
1031 ierr = VecGetBlockSize(vNatural, &bs);CHKERRQ(ierr);
1032 if (bs > 1) {ierr = ISCreateStride(comm, (xe-xs)/bs, xs, bs, &compIS);CHKERRQ(ierr);}
1033 for (set = 0; set < numCS; ++set) {
1034 PetscInt csLocalSize, c;
1035
1036 /* range of indices for set setID[set]: csxs:csxs + csSize[set]-1
1037 local slice of zonal values: xs/bs,xm/bs-1
1038 intersection: max(xs/bs,csxs),min(xm/bs-1,csxs + csSize[set]-1) */
1039 csLocalSize = PetscMax(0, PetscMin(xe/bs, csxs+csSize[set]) - PetscMax(xs/bs, csxs));
1040 if (bs == 1) {
1041 ierr = VecGetArray(vNatural, &varray);CHKERRQ(ierr);
1042 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset, csID[set], PetscMax(xs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs)]));
1043 ierr = VecRestoreArray(vNatural, &varray);CHKERRQ(ierr);
1044 } else {
1045 for (c = 0; c < bs; ++c) {
1046 ierr = ISStrideSetStride(compIS, (xe-xs)/bs, xs+c, bs);CHKERRQ(ierr);
1047 ierr = VecGetSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
1048 ierr = VecGetArray(vComp, &varray);CHKERRQ(ierr);
1049 PetscStackCallStandard(ex_get_partial_var,(exoid, step, EX_ELEM_BLOCK, offset+c, csID[set], PetscMax(xs/bs-csxs, 0)+1, csLocalSize, &varray[PetscMax(0, csxs-xs/bs)]));
1050 ierr = VecRestoreArray(vComp, &varray);CHKERRQ(ierr);
1051 ierr = VecRestoreSubVector(vNatural, compIS, &vComp);CHKERRQ(ierr);
1052 }
1053 }
1054 csxs += csSize[set];
1055 }
1056 ierr = PetscFree2(csID, csSize);CHKERRQ(ierr);
1057 if (bs > 1) {ierr = ISDestroy(&compIS);CHKERRQ(ierr);}
1058 if (useNatural) {
1059 ierr = DMPlexNaturalToGlobalBegin(dm, vNatural, v);CHKERRQ(ierr);
1060 ierr = DMPlexNaturalToGlobalEnd(dm, vNatural, v);CHKERRQ(ierr);
1061 ierr = DMRestoreGlobalVector(dm, &vNatural);CHKERRQ(ierr);
1062 }
1063 PetscFunctionReturn(0);
1064 }
1065 #endif
1066
1067 /*@
1068 PetscViewerExodusIIGetId - Get the file id of the ExodusII file
1069
1070 Logically Collective on PetscViewer
1071
1072 Input Parameter:
1073 . viewer - the PetscViewer
1074
1075 Output Parameter:
1076 - exoid - The ExodusII file id
1077
1078 Level: intermediate
1079
1080 .seealso: PetscViewerFileSetMode(), PetscViewerCreate(), PetscViewerSetType(), PetscViewerBinaryOpen()
1081 @*/
PetscViewerExodusIIGetId(PetscViewer viewer,int * exoid)1082 PetscErrorCode PetscViewerExodusIIGetId(PetscViewer viewer, int *exoid)
1083 {
1084 PetscErrorCode ierr;
1085
1086 PetscFunctionBegin;
1087 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 1);
1088 ierr = PetscTryMethod(viewer, "PetscViewerExodusIIGetId_C",(PetscViewer,int*),(viewer,exoid));CHKERRQ(ierr);
1089 PetscFunctionReturn(0);
1090 }
1091
1092 /*@C
1093 DMPlexCreateExodusFromFile - Create a DMPlex mesh from an ExodusII file.
1094
1095 Collective
1096
1097 Input Parameters:
1098 + comm - The MPI communicator
1099 . filename - The name of the ExodusII file
1100 - interpolate - Create faces and edges in the mesh
1101
1102 Output Parameter:
1103 . dm - The DM object representing the mesh
1104
1105 Level: beginner
1106
1107 .seealso: DMPLEX, DMCreate(), DMPlexCreateExodus()
1108 @*/
DMPlexCreateExodusFromFile(MPI_Comm comm,const char filename[],PetscBool interpolate,DM * dm)1109 PetscErrorCode DMPlexCreateExodusFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1110 {
1111 PetscMPIInt rank;
1112 PetscErrorCode ierr;
1113 #if defined(PETSC_HAVE_EXODUSII)
1114 int CPU_word_size = sizeof(PetscReal), IO_word_size = 0, exoid = -1;
1115 float version;
1116 #endif
1117
1118 PetscFunctionBegin;
1119 PetscValidCharPointer(filename, 2);
1120 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1121 #if defined(PETSC_HAVE_EXODUSII)
1122 if (!rank) {
1123 exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);
1124 if (exoid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "ex_open(\"%s\",...) did not return a valid file ID", filename);
1125 }
1126 ierr = DMPlexCreateExodus(comm, exoid, interpolate, dm);CHKERRQ(ierr);
1127 if (!rank) {PetscStackCallStandard(ex_close,(exoid));}
1128 #else
1129 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1130 #endif
1131 PetscFunctionReturn(0);
1132 }
1133
1134 #if defined(PETSC_HAVE_EXODUSII)
ExodusGetCellType_Private(const char * elem_type,DMPolytopeType * ct)1135 static PetscErrorCode ExodusGetCellType_Private(const char *elem_type, DMPolytopeType *ct)
1136 {
1137 PetscBool flg;
1138 PetscErrorCode ierr;
1139
1140 PetscFunctionBegin;
1141 *ct = DM_POLYTOPE_UNKNOWN;
1142 ierr = PetscStrcmp(elem_type, "TRI", &flg);CHKERRQ(ierr);
1143 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;}
1144 ierr = PetscStrcmp(elem_type, "TRI3", &flg);CHKERRQ(ierr);
1145 if (flg) {*ct = DM_POLYTOPE_TRIANGLE; goto done;}
1146 ierr = PetscStrcmp(elem_type, "QUAD", &flg);CHKERRQ(ierr);
1147 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1148 ierr = PetscStrcmp(elem_type, "QUAD4", &flg);CHKERRQ(ierr);
1149 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1150 ierr = PetscStrcmp(elem_type, "SHELL4", &flg);CHKERRQ(ierr);
1151 if (flg) {*ct = DM_POLYTOPE_QUADRILATERAL; goto done;}
1152 ierr = PetscStrcmp(elem_type, "TETRA", &flg);CHKERRQ(ierr);
1153 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;}
1154 ierr = PetscStrcmp(elem_type, "TET4", &flg);CHKERRQ(ierr);
1155 if (flg) {*ct = DM_POLYTOPE_TETRAHEDRON; goto done;}
1156 ierr = PetscStrcmp(elem_type, "WEDGE", &flg);CHKERRQ(ierr);
1157 if (flg) {*ct = DM_POLYTOPE_TRI_PRISM; goto done;}
1158 ierr = PetscStrcmp(elem_type, "HEX", &flg);CHKERRQ(ierr);
1159 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1160 ierr = PetscStrcmp(elem_type, "HEX8", &flg);CHKERRQ(ierr);
1161 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1162 ierr = PetscStrcmp(elem_type, "HEXAHEDRON", &flg);CHKERRQ(ierr);
1163 if (flg) {*ct = DM_POLYTOPE_HEXAHEDRON; goto done;}
1164 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized element type %s", elem_type);
1165 done:
1166 PetscFunctionReturn(0);
1167 }
1168 #endif
1169
1170 /*@
1171 DMPlexCreateExodus - Create a DMPlex mesh from an ExodusII file ID.
1172
1173 Collective
1174
1175 Input Parameters:
1176 + comm - The MPI communicator
1177 . exoid - The ExodusII id associated with a exodus file and obtained using ex_open
1178 - interpolate - Create faces and edges in the mesh
1179
1180 Output Parameter:
1181 . dm - The DM object representing the mesh
1182
1183 Level: beginner
1184
1185 .seealso: DMPLEX, DMCreate()
1186 @*/
DMPlexCreateExodus(MPI_Comm comm,PetscInt exoid,PetscBool interpolate,DM * dm)1187 PetscErrorCode DMPlexCreateExodus(MPI_Comm comm, PetscInt exoid, PetscBool interpolate, DM *dm)
1188 {
1189 #if defined(PETSC_HAVE_EXODUSII)
1190 PetscMPIInt num_proc, rank;
1191 PetscSection coordSection;
1192 Vec coordinates;
1193 PetscScalar *coords;
1194 PetscInt coordSize, v;
1195 PetscErrorCode ierr;
1196 /* Read from ex_get_init() */
1197 char title[PETSC_MAX_PATH_LEN+1];
1198 int dim = 0, dimEmbed = 0, numVertices = 0, numCells = 0, numHybridCells = 0;
1199 int num_cs = 0, num_vs = 0, num_fs = 0;
1200 #endif
1201
1202 PetscFunctionBegin;
1203 #if defined(PETSC_HAVE_EXODUSII)
1204 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1205 ierr = MPI_Comm_size(comm, &num_proc);CHKERRQ(ierr);
1206 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
1207 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
1208 /* Open EXODUS II file and read basic informations on rank 0, then broadcast to all processors */
1209 if (!rank) {
1210 ierr = PetscMemzero(title,PETSC_MAX_PATH_LEN+1);CHKERRQ(ierr);
1211 PetscStackCallStandard(ex_get_init,(exoid, title, &dimEmbed, &numVertices, &numCells, &num_cs, &num_vs, &num_fs));
1212 if (!num_cs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
1213 }
1214 ierr = MPI_Bcast(title, PETSC_MAX_PATH_LEN+1, MPI_CHAR, 0, comm);CHKERRQ(ierr);
1215 ierr = MPI_Bcast(&dim, 1, MPI_INT, 0, comm);CHKERRQ(ierr);
1216 ierr = PetscObjectSetName((PetscObject) *dm, title);CHKERRQ(ierr);
1217 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr);
1218 /* We do not want this label automatically computed, instead we compute it here */
1219 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr);
1220
1221 /* Read cell sets information */
1222 if (!rank) {
1223 PetscInt *cone;
1224 int c, cs, ncs, c_loc, v, v_loc;
1225 /* Read from ex_get_elem_blk_ids() */
1226 int *cs_id, *cs_order;
1227 /* Read from ex_get_elem_block() */
1228 char buffer[PETSC_MAX_PATH_LEN+1];
1229 int num_cell_in_set, num_vertex_per_cell, num_hybrid, num_attr;
1230 /* Read from ex_get_elem_conn() */
1231 int *cs_connect;
1232
1233 /* Get cell sets IDs */
1234 ierr = PetscMalloc2(num_cs, &cs_id, num_cs, &cs_order);CHKERRQ(ierr);
1235 PetscStackCallStandard(ex_get_ids,(exoid, EX_ELEM_BLOCK, cs_id));
1236 /* Read the cell set connectivity table and build mesh topology
1237 EXO standard requires that cells in cell sets be numbered sequentially and be pairwise disjoint. */
1238 /* Check for a hybrid mesh */
1239 for (cs = 0, num_hybrid = 0; cs < num_cs; ++cs) {
1240 DMPolytopeType ct;
1241 char elem_type[PETSC_MAX_PATH_LEN];
1242
1243 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr);
1244 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type));
1245 ierr = ExodusGetCellType_Private(elem_type, &ct);CHKERRQ(ierr);
1246 dim = PetscMax(dim, DMPolytopeTypeGetDim(ct));
1247 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1248 switch (ct) {
1249 case DM_POLYTOPE_TRI_PRISM:
1250 cs_order[cs] = cs;
1251 numHybridCells += num_cell_in_set;
1252 ++num_hybrid;
1253 break;
1254 default:
1255 for (c = cs; c > cs-num_hybrid; --c) cs_order[c] = cs_order[c-1];
1256 cs_order[cs-num_hybrid] = cs;
1257 }
1258 }
1259 /* First set sizes */
1260 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1261 DMPolytopeType ct;
1262 char elem_type[PETSC_MAX_PATH_LEN];
1263 const PetscInt cs = cs_order[ncs];
1264
1265 ierr = PetscArrayzero(elem_type, sizeof(elem_type));CHKERRQ(ierr);
1266 PetscStackCallStandard(ex_get_elem_type,(exoid, cs_id[cs], elem_type));
1267 ierr = ExodusGetCellType_Private(elem_type, &ct);CHKERRQ(ierr);
1268 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set,&num_vertex_per_cell, 0, 0, &num_attr));
1269 for (c_loc = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1270 ierr = DMPlexSetConeSize(*dm, c, num_vertex_per_cell);CHKERRQ(ierr);
1271 ierr = DMPlexSetCellType(*dm, c, ct);CHKERRQ(ierr);
1272 }
1273 }
1274 for (v = numCells; v < numCells+numVertices; ++v) {ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);}
1275 ierr = DMSetUp(*dm);CHKERRQ(ierr);
1276 for (ncs = 0, c = 0; ncs < num_cs; ++ncs) {
1277 const PetscInt cs = cs_order[ncs];
1278 PetscStackCallStandard(ex_get_block,(exoid, EX_ELEM_BLOCK, cs_id[cs], buffer, &num_cell_in_set, &num_vertex_per_cell, 0, 0, &num_attr));
1279 ierr = PetscMalloc2(num_vertex_per_cell*num_cell_in_set,&cs_connect,num_vertex_per_cell,&cone);CHKERRQ(ierr);
1280 PetscStackCallStandard(ex_get_conn,(exoid, EX_ELEM_BLOCK, cs_id[cs], cs_connect,NULL,NULL));
1281 /* EXO uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
1282 for (c_loc = 0, v = 0; c_loc < num_cell_in_set; ++c_loc, ++c) {
1283 DMPolytopeType ct;
1284
1285 for (v_loc = 0; v_loc < num_vertex_per_cell; ++v_loc, ++v) {
1286 cone[v_loc] = cs_connect[v]+numCells-1;
1287 }
1288 ierr = DMPlexGetCellType(*dm, c, &ct);CHKERRQ(ierr);
1289 ierr = DMPlexInvertCell(ct, cone);CHKERRQ(ierr);
1290 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr);
1291 ierr = DMSetLabelValue(*dm, "Cell Sets", c, cs_id[cs]);CHKERRQ(ierr);
1292 }
1293 ierr = PetscFree2(cs_connect,cone);CHKERRQ(ierr);
1294 }
1295 ierr = PetscFree2(cs_id, cs_order);CHKERRQ(ierr);
1296 }
1297 {
1298 PetscInt ints[] = {dim, dimEmbed};
1299
1300 ierr = MPI_Bcast(ints, 2, MPIU_INT, 0, comm);CHKERRQ(ierr);
1301 ierr = DMSetDimension(*dm, ints[0]);CHKERRQ(ierr);
1302 ierr = DMSetCoordinateDim(*dm, ints[1]);CHKERRQ(ierr);
1303 dim = ints[0];
1304 dimEmbed = ints[1];
1305 }
1306 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
1307 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
1308 if (interpolate) {
1309 DM idm;
1310
1311 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
1312 ierr = DMDestroy(dm);CHKERRQ(ierr);
1313 *dm = idm;
1314 }
1315
1316 /* Create vertex set label */
1317 if (!rank && (num_vs > 0)) {
1318 int vs, v;
1319 /* Read from ex_get_node_set_ids() */
1320 int *vs_id;
1321 /* Read from ex_get_node_set_param() */
1322 int num_vertex_in_set;
1323 /* Read from ex_get_node_set() */
1324 int *vs_vertex_list;
1325
1326 /* Get vertex set ids */
1327 ierr = PetscMalloc1(num_vs, &vs_id);CHKERRQ(ierr);
1328 PetscStackCallStandard(ex_get_ids,(exoid, EX_NODE_SET, vs_id));
1329 for (vs = 0; vs < num_vs; ++vs) {
1330 PetscStackCallStandard(ex_get_set_param,(exoid, EX_NODE_SET, vs_id[vs], &num_vertex_in_set, NULL));
1331 ierr = PetscMalloc1(num_vertex_in_set, &vs_vertex_list);CHKERRQ(ierr);
1332 PetscStackCallStandard(ex_get_set,(exoid, EX_NODE_SET, vs_id[vs], vs_vertex_list, NULL));
1333 for (v = 0; v < num_vertex_in_set; ++v) {
1334 ierr = DMSetLabelValue(*dm, "Vertex Sets", vs_vertex_list[v]+numCells-1, vs_id[vs]);CHKERRQ(ierr);
1335 }
1336 ierr = PetscFree(vs_vertex_list);CHKERRQ(ierr);
1337 }
1338 ierr = PetscFree(vs_id);CHKERRQ(ierr);
1339 }
1340 /* Read coordinates */
1341 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
1342 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
1343 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr);
1344 ierr = PetscSectionSetChart(coordSection, numCells, numCells + numVertices);CHKERRQ(ierr);
1345 for (v = numCells; v < numCells+numVertices; ++v) {
1346 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr);
1347 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr);
1348 }
1349 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
1350 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
1351 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
1352 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
1353 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
1354 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr);
1355 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
1356 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
1357 if (!rank) {
1358 PetscReal *x, *y, *z;
1359
1360 ierr = PetscMalloc3(numVertices,&x,numVertices,&y,numVertices,&z);CHKERRQ(ierr);
1361 PetscStackCallStandard(ex_get_coord,(exoid, x, y, z));
1362 if (dimEmbed > 0) {
1363 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+0] = x[v];
1364 }
1365 if (dimEmbed > 1) {
1366 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+1] = y[v];
1367 }
1368 if (dimEmbed > 2) {
1369 for (v = 0; v < numVertices; ++v) coords[v*dimEmbed+2] = z[v];
1370 }
1371 ierr = PetscFree3(x,y,z);CHKERRQ(ierr);
1372 }
1373 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
1374 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
1375 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
1376
1377 /* Create side set label */
1378 if (!rank && interpolate && (num_fs > 0)) {
1379 int fs, f, voff;
1380 /* Read from ex_get_side_set_ids() */
1381 int *fs_id;
1382 /* Read from ex_get_side_set_param() */
1383 int num_side_in_set;
1384 /* Read from ex_get_side_set_node_list() */
1385 int *fs_vertex_count_list, *fs_vertex_list;
1386 /* Read side set labels */
1387 char fs_name[MAX_STR_LENGTH+1];
1388
1389 /* Get side set ids */
1390 ierr = PetscMalloc1(num_fs, &fs_id);CHKERRQ(ierr);
1391 PetscStackCallStandard(ex_get_ids,(exoid, EX_SIDE_SET, fs_id));
1392 for (fs = 0; fs < num_fs; ++fs) {
1393 PetscStackCallStandard(ex_get_set_param,(exoid, EX_SIDE_SET, fs_id[fs], &num_side_in_set, NULL));
1394 ierr = PetscMalloc2(num_side_in_set,&fs_vertex_count_list,num_side_in_set*4,&fs_vertex_list);CHKERRQ(ierr);
1395 PetscStackCallStandard(ex_get_side_set_node_list,(exoid, fs_id[fs], fs_vertex_count_list, fs_vertex_list));
1396 /* Get the specific name associated with this side set ID. */
1397 int fs_name_err = ex_get_name(exoid, EX_SIDE_SET, fs_id[fs], fs_name);
1398 for (f = 0, voff = 0; f < num_side_in_set; ++f) {
1399 const PetscInt *faces = NULL;
1400 PetscInt faceSize = fs_vertex_count_list[f], numFaces;
1401 PetscInt faceVertices[4], v;
1402
1403 if (faceSize > 4) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "ExodusII side cannot have %d > 4 vertices", faceSize);
1404 for (v = 0; v < faceSize; ++v, ++voff) {
1405 faceVertices[v] = fs_vertex_list[voff]+numCells-1;
1406 }
1407 ierr = DMPlexGetFullJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr);
1408 if (numFaces != 1) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "Invalid ExodusII side %d in set %d maps to %d faces", f, fs, numFaces);
1409 ierr = DMSetLabelValue(*dm, "Face Sets", faces[0], fs_id[fs]);CHKERRQ(ierr);
1410 /* Only add the label if one has been detected for this side set. */
1411 if (!fs_name_err) {
1412 ierr = DMSetLabelValue(*dm, fs_name, faces[0], fs_id[fs]);CHKERRQ(ierr);
1413 }
1414 ierr = DMPlexRestoreJoin(*dm, faceSize, faceVertices, &numFaces, &faces);CHKERRQ(ierr);
1415 }
1416 ierr = PetscFree2(fs_vertex_count_list,fs_vertex_list);CHKERRQ(ierr);
1417 }
1418 ierr = PetscFree(fs_id);CHKERRQ(ierr);
1419 }
1420 #else
1421 SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1422 #endif
1423 PetscFunctionReturn(0);
1424 }
1425