Actual source code: meshexodus.c
petsc-3.4.5 2014-06-29
1: #include <petscdmmesh_formats.hh> /*I "petscdmmesh.h" I*/
2: #include <petsc-private/vecimpl.h>
4: #if defined(PETSC_HAVE_EXODUSII)
5: #include <netcdf.h>
6: #include <exodusII.h>
10: PetscErrorCode PetscReadExodusII(MPI_Comm comm, const char filename[], ALE::Obj<PETSC_MESH_TYPE>& mesh)
11: {
12: int exoid;
13: int CPU_word_size = 0, IO_word_size = 0;
14: const PetscMPIInt rank = mesh->commRank();
15: float version;
16: char title[PETSC_MAX_PATH_LEN+1], elem_type[PETSC_MAX_PATH_LEN+1];
17: int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets;
18: int ierr;
21: // Open EXODUS II file
22: exoid = ex_open(filename, EX_READ, &CPU_word_size, &IO_word_size, &version);CHKERRQ(!exoid);
23: // Read database parameters
24: ex_get_init(exoid, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets);
26: // Read vertex coordinates
27: float *x, *y, *z;
28: PetscMalloc3(num_nodes,float,&x,num_nodes,float,&y,num_nodes,float,&z);
29: ex_get_coord(exoid, x, y, z);
31: // Read element connectivity
32: int *eb_ids, *num_elem_in_block, *num_nodes_per_elem, *num_attr;
33: int **connect;
34: char **block_names;
35: if (num_elem_blk > 0) {
36: PetscMalloc5(num_elem_blk,int,&eb_ids,num_elem_blk,int,&num_elem_in_block,num_elem_blk,int,&num_nodes_per_elem,num_elem_blk,int,&num_attr,num_elem_blk,char*,&block_names);
37: ex_get_elem_blk_ids(exoid, eb_ids);
38: for (int eb = 0; eb < num_elem_blk; ++eb) {
39: PetscMalloc((PETSC_MAX_PATH_LEN+1) * sizeof(char), &block_names[eb]);
40: }
41: ex_get_names(exoid, EX_ELEM_BLOCK, block_names);
42: for (int eb = 0; eb < num_elem_blk; ++eb) {
43: ex_get_elem_block(exoid, eb_ids[eb], elem_type, &num_elem_in_block[eb], &num_nodes_per_elem[eb], &num_attr[eb]);
44: PetscFree(block_names[eb]);
45: }
46: PetscMalloc(num_elem_blk * sizeof(int*),&connect);
47: for (int eb = 0; eb < num_elem_blk; ++eb) {
48: if (num_elem_in_block[eb] > 0) {
49: PetscMalloc(num_nodes_per_elem[eb]*num_elem_in_block[eb] * sizeof(int),&connect[eb]);
50: ex_get_elem_conn(exoid, eb_ids[eb], connect[eb]);
51: }
52: }
53: }
55: // Read node sets
56: int *ns_ids, *num_nodes_in_set;
57: int **node_list;
58: if (num_node_sets > 0) {
59: PetscMalloc3(num_node_sets,int,&ns_ids,num_node_sets,int,&num_nodes_in_set,num_node_sets,int*,&node_list);
60: ex_get_node_set_ids(exoid, ns_ids);
61: for (int ns = 0; ns < num_node_sets; ++ns) {
62: int num_df_in_set;
63: ex_get_node_set_param (exoid, ns_ids[ns], &num_nodes_in_set[ns], &num_df_in_set);
64: PetscMalloc(num_nodes_in_set[ns] * sizeof(int), &node_list[ns]);
65: ex_get_node_set(exoid, ns_ids[ns], node_list[ns]);
66: }
67: }
68: ex_close(exoid);
70: // Build mesh topology
71: int numCorners = num_nodes_per_elem[0];
72: int *cells;
73: mesh->setDimension(num_dim);
74: PetscMalloc(numCorners*num_elem * sizeof(int), &cells);
75: for (int eb = 0, k = 0; eb < num_elem_blk; ++eb) {
76: for (int e = 0; e < num_elem_in_block[eb]; ++e, ++k) {
77: for (int c = 0; c < numCorners; ++c) {
78: cells[k*numCorners+c] = connect[eb][e*numCorners+c];
79: }
80: }
81: PetscFree(connect[eb]);
82: }
83: ALE::Obj<PETSC_MESH_TYPE::sieve_type> sieve = new PETSC_MESH_TYPE::sieve_type(mesh->comm(), mesh->debug());
84: bool interpolate = false;
86: try {
87: mesh->setSieve(sieve);
88: if (0 == rank) {
89: if (!interpolate) {
90: // Create the ISieve
91: sieve->setChart(PETSC_MESH_TYPE::sieve_type::chart_type(0, num_elem+num_nodes));
92: // Set cone and support sizes
93: for (int c = 0; c < num_elem; ++c) {
94: sieve->setConeSize(c, numCorners);
95: }
96: sieve->symmetrizeSizes(num_elem, numCorners, cells, num_elem - 1); /* Notice the -1 for 1-based indexing in cells[] */
97: // Allocate point storage
98: sieve->allocate();
99: // Fill up cones
100: int *cone = new int[numCorners];
101: int *coneO = new int[numCorners];
102: for (int v = 0; v < numCorners; ++v) coneO[v] = 1;
103: for (int c = 0; c < num_elem; ++c) {
104: for (int v = 0; v < numCorners; ++v) {
105: cone[v] = cells[c*numCorners+v]+num_elem - 1;
106: }
107: sieve->setCone(cone, c);
108: sieve->setConeOrientation(coneO, c);
109: } // for
110: delete[] cone; cone = 0;
111: delete[] coneO; coneO = 0;
112: // Symmetrize to fill up supports
113: sieve->symmetrize();
114: } else {
115: // Same old thing
116: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
117: ALE::Obj<FlexMesh::sieve_type> s = new FlexMesh::sieve_type(sieve->comm(), sieve->debug());
119: ALE::SieveBuilder<FlexMesh>::buildTopology(s, num_dim, num_elem, cells, num_nodes, interpolate, numCorners);
120: std::map<FlexMesh::point_type,FlexMesh::point_type> renumbering;
121: ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
122: }
123: if (!interpolate) {
124: // Optimized stratification
125: const ALE::Obj<PETSC_MESH_TYPE::label_type>& height = mesh->createLabel("height");
126: const ALE::Obj<PETSC_MESH_TYPE::label_type>& depth = mesh->createLabel("depth");
128: for (int c = 0; c < num_elem; ++c) {
129: height->setCone(0, c);
130: depth->setCone(1, c);
131: }
132: for (int v = num_elem; v < num_elem+num_nodes; ++v) {
133: height->setCone(1, v);
134: depth->setCone(0, v);
135: }
136: mesh->setHeight(1);
137: mesh->setDepth(1);
138: } else {
139: mesh->stratify();
140: }
141: } else {
142: mesh->getSieve()->setChart(PETSC_MESH_TYPE::sieve_type::chart_type());
143: mesh->getSieve()->allocate();
144: mesh->stratify();
145: }
146: } catch (ALE::Exception e) {
147: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB, e.msg().c_str());
148: }
149: PetscFree(cells);
151: // Build cell blocks
152: const ALE::Obj<PETSC_MESH_TYPE::label_type>& cellBlocks = mesh->createLabel("CellBlocks");
153: if (!rank) {
154: for (int eb = 0, k = 0; eb < num_elem_blk; ++eb) {
155: for (int e = 0; e < num_elem_in_block[eb]; ++e, ++k) {
156: mesh->setValue(cellBlocks, k, eb_ids[eb]);
157: }
158: }
159: }
160: if (num_elem_blk > 0) {
161: PetscFree(connect);
162: PetscFree5(eb_ids, num_elem_in_block, num_nodes_per_elem, num_attr, block_names);
163: }
165: // Build coordinates
166: double *coords;
167: PetscMalloc(num_dim*num_nodes * sizeof(double), &coords);
168: if (num_dim > 0) {
169: for (int v = 0; v < num_nodes; ++v) coords[v*num_dim+0] = x[v];
170: }
171: if (num_dim > 1) {
172: for (int v = 0; v < num_nodes; ++v) coords[v*num_dim+1] = y[v];
173: }
174: if (num_dim > 2) {
175: for (int v = 0; v < num_nodes; ++v) coords[v*num_dim+2] = z[v];
176: }
177: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, num_dim, coords);
178: PetscFree(coords);
179: PetscFree3(x, y, z);
181: // Build vertex sets
182: const ALE::Obj<PETSC_MESH_TYPE::label_type>& vertexSets = mesh->createLabel("VertexSets");
183: if (!rank) {
184: for (int ns = 0; ns < num_node_sets; ++ns) {
185: for (int n = 0; n < num_nodes_in_set[ns]; ++n) {
186: mesh->addValue(vertexSets, node_list[ns][n]+num_elem-1, ns_ids[ns]);
187: }
188: PetscFree(node_list[ns]);
189: }
190: }
191: if (num_node_sets > 0) {
192: PetscFree3(ns_ids,num_nodes_in_set,node_list);
193: }
194: return(0);
195: }
197: #endif // PETSC_HAVE_EXODUSII
201: /*@C
202: DMMeshCreateExodus - Create a DMMesh from an ExodusII file.
204: Not Collective
206: Input Parameters:
207: + comm - The MPI communicator
208: - filename - The ExodusII filename
210: Output Parameter:
211: . dm - The DMMesh object
213: Level: beginner
215: .keywords: mesh, ExodusII
216: .seealso: DMMeshCreate()
217: @*/
218: PetscErrorCode DMMeshCreateExodus(MPI_Comm comm, const char filename[], DM *dm)
219: {
220: PetscInt debug = 0;
221: PetscBool flag;
225: DMMeshCreate(comm, dm);
226: PetscOptionsGetInt(NULL, "-debug", &debug, &flag);
227: ALE::Obj<PETSC_MESH_TYPE> m = new PETSC_MESH_TYPE(comm, -1, debug);
228: #if defined(PETSC_HAVE_EXODUSII)
229: PetscReadExodusII(comm, filename, m);
230: #else
231: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
232: #endif
233: if (debug) m->view("Mesh");
234: DMMeshSetMesh(*dm, m);
235: return(0);
236: }
240: /*@
241: DMMeshCreateExodusNG - Create a Mesh from an ExodusII file.
243: Collective on comm
245: Input Parameters:
246: + comm - The MPI communicator
247: - exoid - The ExodusII id associated with a exodus file and obtained using ex_open
249: Output Parameter:
250: . dm - The DM object representing the mesh
252: ExodusII side sets are ignored
255: Interpolated meshes are not supported.
257: Level: beginner
259: .keywords: mesh,ExodusII
260: .seealso: MeshCreate() MeshCreateExodus()
261: @*/
262: PetscErrorCode DMMeshCreateExodusNG(MPI_Comm comm, PetscInt exoid,DM *dm)
263: {
264: #if defined(PETSC_HAVE_EXODUSII)
265: PetscBool debug = PETSC_FALSE;
266: PetscMPIInt num_proc,rank;
269: int num_dim,num_vertices = 0,num_cell = 0;
270: int num_cs = 0,num_vs = 0,num_fs = 0;
271: char title[PETSC_MAX_PATH_LEN+1];
272: char buffer[PETSC_MAX_PATH_LEN+1];
274: int *cs_id;
275: int num_cell_in_set,num_vertex_per_cell,num_attr;
276: int *cs_connect;
278: int *vs_id;
279: int num_vertex_in_set;
280: int *vs_vertex_list;
281: PetscReal *x,*y,*z;
282: PetscReal *coords;
284: PetscInt v,c,v_loc,c_loc,vs,cs;
286: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
287: typedef PETSC_MESH_TYPE::sieve_type sieve_type;
288: ALE::Obj<PETSC_MESH_TYPE> mesh;
289: ALE::Obj<FlexMesh::sieve_type> flexmesh_sieve;
290: ALE::Obj<FlexMesh> flexmesh;
291: ALE::Obj<PETSC_MESH_TYPE::sieve_type> sieve;
293: std::map<PETSC_MESH_TYPE::point_type,PETSC_MESH_TYPE::point_type> renumbering;
294: #endif
297: #if defined(PETSC_HAVE_EXODUSII)
298: MPI_Comm_rank(comm,&rank);
299: MPI_Comm_size(comm,&num_proc);
300: PetscOptionsGetBool(NULL,"-debug",&debug,NULL);
302: DMMeshCreate(comm,dm);
303: /*
304: I _really don't understand how this is needed
305: */
306: PetscObjectGetComm((PetscObject) *dm,&comm);
309: mesh = new PETSC_MESH_TYPE(comm,-1,debug);
310: DMMeshSetMesh(*dm,mesh);
312: sieve = new PETSC_MESH_TYPE::sieve_type(mesh->comm(),mesh->debug());
313: flexmesh_sieve = new FlexMesh::sieve_type(mesh->comm(),mesh->debug());
315: /*
316: Open EXODUS II file and read basic informations on rank 0,
317: then broadcast to all processors
318: */
319: if (!rank) {
320: ex_get_init(exoid,title,&num_dim,&num_vertices,&num_cell,&num_cs,&num_vs,&num_fs);
321: if (num_cs == 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Exodus file does not contain any cell set\n");
322: }
324: MPI_Bcast(&num_dim,1,MPIU_INT,0,comm);
325: MPI_Bcast(&num_cs,1,MPIU_INT,0,comm);
326: MPI_Bcast(&num_vs,1,MPIU_INT,0,comm);
327: MPI_Bcast(&num_fs,1,MPIU_INT,0,comm);
328: MPI_Bcast(title,PETSC_MAX_PATH_LEN+1,MPI_CHAR,0,comm);
330: PetscObjectSetName((PetscObject)*dm,title);
331: mesh->setDimension(num_dim);
333: /*
334: Read cell sets information then broadcast them
335: */
336: const ALE::Obj<PETSC_MESH_TYPE::label_type>& cellSets = mesh->createLabel("Cell Sets");
337: if (!rank) {
338: /*
339: Get cell sets IDs
340: */
341: PetscMalloc(num_cs*sizeof(int),&cs_id);
342: ex_get_elem_blk_ids(exoid,cs_id);
344: /*
345: Read the cell set connectivity table and build mesh topology
346: EXO standard requires that cells in cell sets be numbered sequentially
347: and be pairwise disjoint.
348: */
349: for (c=0,cs = 0; cs < num_cs; cs++) {
350: if (debug) {
351: PetscPrintf(comm,"Building cell set %i\n",cs);
352: }
353: ex_get_elem_block(exoid,cs_id[cs],buffer,&num_cell_in_set,&num_vertex_per_cell,&num_attr);
354: PetscMalloc(num_vertex_per_cell*num_cell_in_set*sizeof(int),&cs_connect);
355: ex_get_elem_conn(exoid,cs_id[cs],cs_connect);
356: /*
357: EXO uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices.
358: */
359: for (v = 0,c_loc = 0; c_loc < num_cell_in_set; c_loc++,c++) {
360: for (v_loc = 0; v_loc < num_vertex_per_cell; v_loc++,v++) {
361: if (debug) {
362: PetscPrintf(PETSC_COMM_SELF,"[%i]:\tinserting vertex \t%i in cell \t%i local vertex number \t%i\n",
363: rank,cs_connect[v]+num_cell-1,c,v_loc);
364: }
365: flexmesh_sieve->addArrow(sieve_type::point_type(cs_connect[v]+num_cell-1),
366: sieve_type::point_type(c),
367: v_loc);
368: }
369: mesh->setValue(cellSets,c,cs_id[cs]);
370: }
371: PetscFree(cs_connect);
372: }
373: PetscFree(cs_id);
374: }
375: /*
376: We do not interpolate and know that the numbering is compact (this is required in exo)
377: so no need to renumber (renumber = false)
378: */
379: ALE::ISieveConverter::convertSieve(*flexmesh_sieve,*sieve,renumbering,false);
380: mesh->setSieve(sieve);
381: mesh->stratify();
382: flexmesh = new FlexMesh(mesh->comm(),mesh->debug());
383: ALE::ISieveConverter::convertOrientation(*flexmesh_sieve,*sieve,renumbering,flexmesh->getArrowSection("orientation").ptr());
385: /*
386: Create Vertex set label
387: */
388: const ALE::Obj<PETSC_MESH_TYPE::label_type>& vertexSets = mesh->createLabel("Vertex Sets");
389: if (num_vs >0) {
390: if (!rank) {
391: /*
392: Get vertex set ids
393: */
394: PetscMalloc(num_vs*sizeof(int),&vs_id);
395: ex_get_node_set_ids(exoid,vs_id);
396: for (vs = 0; vs < num_vs; vs++) {
397: ex_get_node_set_param(exoid,vs_id[vs],&num_vertex_in_set,&num_attr);
398: PetscMalloc(num_vertex_in_set * sizeof(int),&vs_vertex_list);
399: ex_get_node_set(exoid,vs_id[vs],vs_vertex_list);
400: for (v = 0; v < num_vertex_in_set; v++) {
401: mesh->addValue(vertexSets,vs_vertex_list[v]+num_cell-1,vs_id[vs]);
402: }
403: PetscFree(vs_vertex_list);
404: }
405: PetscFree(vs_id);
406: }
407: }
408: /*
409: Read coordinates
410: */
411: PetscMalloc4(num_vertices,PetscReal,&x,
412: num_vertices,PetscReal,&y,
413: num_vertices,PetscReal,&z,
414: num_dim*num_vertices,PetscReal,&coords);
415: if (!rank) {
416: ex_get_coord(exoid,x,y,z);
417: PetscMalloc(num_dim*num_vertices * sizeof(PetscReal), &coords);
418: if (num_dim > 0) {
419: for (v = 0; v < num_vertices; ++v) coords[v*num_dim+0] = x[v];
420: }
421: if (num_dim > 1) {
422: for (v = 0; v < num_vertices; ++v) coords[v*num_dim+1] = y[v];
423: }
424: if (num_dim > 2) {
425: for (v = 0; v < num_vertices; ++v) coords[v*num_dim+2] = z[v];
426: }
427: }
428: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh,num_dim,coords);
429: if (!rank) {
430: PetscFree4(x,y,z,coords);
431: }
433: /*
434: if (!rank) {
435: ex_close(exoid);
436: }
437: */
438: #else
439: SETERRQ(comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
440: #endif
441: return(0);
442: }
446: /*@
447: DMMeshViewExodusSplit - Write a dmMesh geometry and topology into several ExodusII files.
449: Collective on comm
451: Input Parameters:
452: + comm - The MPI communicator
453: - filename - The ExodusII filename. Must be different on each processor
454: (suggest using prefix-<rank>.gen or prefix-<rank>.exo)
455: . dm - The DM object representing the body
457: Face Sets (Side Sets in Exodus terminology) are ignored
459: Interpolated meshes are not supported.
461: Level: beginner
463: .keywords: mesh,ExodusII
464: .seealso: MeshCreate()
465: @*/
466: PetscErrorCode DMMeshViewExodusSplit(DM dm,PetscInt exoid)
467: {
468: #if defined(PETSC_HAVE_EXODUSII)
469: PetscBool debug = PETSC_FALSE;
470: MPI_Comm comm;
473: int num_dim,num_vertices = 0,num_cells = 0;
474: int num_cs = 0,num_vs = 0;
475: int num_cs_global = 0,num_vs_global = 0;
476: const char *title;
477: IS csIS,vsIS;
478: IS csIS_global,vsIS_global;
479: const PetscInt *csID,*vsID,*labels;
481: PetscReal *coord;
483: PetscInt c,v,c_offset;
485: PetscInt set,num_cell_in_set,num_vertex_per_cell;
486: const PetscInt *cells;
487: IS cellsIS;
488: const char *cell_type;
489: int *elem_map,*cs_connect,num_cs_attr=0;
490: const PetscInt *elem_connect;
492: PetscInt num_vertex_in_set,num_vs_attr=0;
493: PetscInt *vertices;
494: IS verticesIS;
495: #endif
498: #if defined(PETSC_HAVE_EXODUSII)
499: /*
500: Extract mesh global properties from the DM
501: */
502: PetscObjectGetName((PetscObject)dm,&title);
503: DMMeshGetDimension(dm,&num_dim);
504: DMMeshGetStratumSize(dm,"height",0,&num_cells);
505: DMMeshGetStratumSize(dm,"depth",0,&num_vertices);
507: /*
508: Get the local and global number of sets
509: */
510: PetscObjectGetComm((PetscObject) dm,&comm);
511: DMMeshGetLabelSize(dm,"Cell Sets",&num_cs);
512: DMMeshGetLabelIdIS(dm,"Cell Sets",&csIS);
513: ISGetTotalIndices(csIS,&labels);
514: ISGetSize(csIS,&num_cs_global);
515: PetscSortRemoveDupsInt(&num_cs_global,(PetscInt*)labels);
516: ISCreateGeneral(comm,num_cs_global,labels,PETSC_COPY_VALUES,&csIS_global);
517: ISRestoreTotalIndices(csIS,&labels);
519: DMMeshGetLabelSize(dm,"Vertex Sets",&num_vs);
520: DMMeshGetLabelIdIS(dm,"Vertex Sets",&vsIS);
521: ISGetSize(vsIS,&num_vs_global);
522: if (num_vs_global > 0) {
523: ISGetTotalIndices(vsIS,&labels);
524: PetscSortRemoveDupsInt(&num_vs_global,(PetscInt*)labels);
525: ISCreateGeneral(comm,num_vs_global,labels,PETSC_COPY_VALUES,&vsIS_global);
526: ISRestoreTotalIndices(vsIS,&labels);
527: }
528: ex_put_init(exoid,title,num_dim,num_vertices,num_cells,num_cs_global,num_vs_global,0);
530: /*
531: Write coordinates
532: */
533: DMMeshGetCoordinates(dm,PETSC_TRUE,&num_vertices,&num_dim,&coord);
534: if (debug) {
535: for (c = 0; c < num_dim; c++) {
536: PetscPrintf(comm,"Coordinate %i:\n",c);
537: PetscRealView(num_vertices,&coord[c*num_vertices],PETSC_VIEWER_STDOUT_SELF);
538: }
539: }
540: ex_put_coord(exoid,coord,&coord[num_vertices],&coord[2*num_vertices]);
542: /*
543: Write cell set connectivity table and parameters
544: Compute the element number map
545: The element number map is not needed as long as the cell sets are of the type
546: 0 .. n1
547: n1+1 .. n2
548: n2+1 ..n3
549: which is the case when a mesh is read from an exo file then distributed, but one never knows
550: */
552: /*
553: The following loop should be based on csIS and not csIS_global, but EXO has no
554: way to write the block id's other than ex_put_elem_block
555: and ensight is bothered if all cell set ID's are not on all files...
556: Not a huge deal
557: */
559: PetscMalloc(num_cells*sizeof(int),&elem_map);
560: ISGetIndices(csIS_global,&csID);
561: for (c_offset = 0,set = 0; set < num_cs_global; set++) {
562: DMMeshGetStratumSize(dm,"Cell Sets",csID[set],&num_cell_in_set);
563: DMMeshGetStratumIS(dm,"Cell Sets",csID[set],&cellsIS);
564: ISGetIndices(cellsIS,&cells);
565: if (debug) {
566: PetscPrintf(PETSC_COMM_SELF,"Cell set %i: %i cells\n",csID[set],num_cell_in_set);
567: PetscIntView(num_cell_in_set,cells,PETSC_VIEWER_STDOUT_SELF);
568: }
569: /*
570: Add current block indices to elem_map. EXO uses fortran numbering
571: */
572: for (c = 0; c < num_cell_in_set; c++,c_offset++) {
573: elem_map[c_offset] = cells[c]+1;
574: }
575: /*
576: We make an educated guess as to the type of cell. This misses out quads in
577: a three-dimensional mesh.
578: This most likely needs to be fixed by calling again ex_put_elem_block with
579: the proper parameters
580: */
581: if (num_cell_in_set > 0) {
582: DMMeshGetConeSize(dm,cells[0],&num_vertex_per_cell);
583: if (num_vertex_per_cell == 2) {
584: cell_type = "BAR";
585: } else if (num_vertex_per_cell == 3) {
586: cell_type = "TRI";
587: } else if (num_vertex_per_cell == 8) {
588: cell_type = "HEX";
589: } else if (num_vertex_per_cell == 4 && num_dim == 2) {
590: cell_type = "QUAD";
591: } else if (num_vertex_per_cell == 4 && num_dim == 3) {
592: cell_type = "TET";
593: } else {
594: cell_type = "UNKNOWN";
595: }
596: }
597: ex_put_elem_block (exoid,csID[set],cell_type,num_cell_in_set,num_vertex_per_cell,num_cs_attr);
598: /*
599: Build connectivity table of the current block
600: */
601: PetscMalloc(num_cell_in_set*num_vertex_per_cell*sizeof(int),&cs_connect);
602: for (c = 0; c < num_cell_in_set; c++) {
603: DMMeshGetCone(dm,cells[c],&elem_connect);
604: for (v = 0; v < num_vertex_per_cell; v++) {
605: cs_connect[c*num_vertex_per_cell+v] = elem_connect[v]-num_cells+1;
606: }
607: }
608: if (num_cell_in_set > 0) {
609: ex_put_elem_conn(exoid,csID[set],cs_connect);
610: }
611: PetscFree(cs_connect);
612: ISRestoreIndices(cellsIS,&cells);
613: ISDestroy(&cellsIS);
614: }
615: ISRestoreIndices(csIS_global,&csID);
616: ex_put_elem_num_map(exoid,elem_map);
617: PetscFree(elem_map);
619: /*
620: Writing vertex sets
621: */
622: if (num_vs_global > 0) {
623: ISGetIndices(vsIS_global,&vsID);
624: for (set = 0; set < num_vs_global; set++) {
625: DMMeshGetStratumSize(dm,"Vertex Sets",vsID[set],&num_vertex_in_set);
626: ex_put_node_set_param(exoid,vsID[set],num_vertex_in_set,num_vs_attr);
628: DMMeshGetStratumIS(dm,"Vertex Sets",vsID[set],&verticesIS);
629: ISGetIndices(verticesIS,(const PetscInt**)&vertices);
631: for (v = 0; v < num_vertex_in_set; v++) {
632: vertices[v] -= num_cells-1;
633: }
635: if (num_vertex_in_set > 0) {
636: ex_put_node_set(exoid,vsID[set],vertices);
637: }
638: ISRestoreIndices(verticesIS,(const PetscInt**)&vertices);
639: ISDestroy(&verticesIS);
640: }
641: ISRestoreIndices(vsIS_global,&vsID);
642: }
643: ISDestroy(&csIS);
644: ISDestroy(&vsIS);
645: ISDestroy(&csIS_global);
646: #else
647: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
648: #endif
649: return(0);
650: }
654: /*@
655: DMMeshCreateScatterToZeroVertex - Creates the scatter required to scatter local (ghosted)
656: Vecs associated with a field defined at the vertices of a mesh to a Vec on cpu 0.
658: Input parameters:
659: . dm - the DMMesh representing the mesh
661: Output parameters:
662: . scatter: the scatter
664: Level: advanced
666: .keywords: mesh,ExodusII
667: .seealso DMMeshCreateScatterToZeroCell DMMeshCreateScatterToZeroVertexSet DMMeshCreateScatterToZeroCellSet
668: @*/
669: PetscErrorCode DMMeshCreateScatterToZeroVertex(DM dm,VecScatter *scatter)
670: {
671: PetscErrorCode ierr;
672: PetscInt i;
673: PetscInt *vertices;
674: PetscInt num_vertices,num_vertices_global,my_num_vertices_global;
675: PetscInt num_cells,num_cells_global,my_num_cells_global;
676: Vec v_local,v_zero;
677: MPI_Comm comm;
678: int rank;
679: IS is_local;
680: ALE::Obj<PETSC_MESH_TYPE> mesh;
681: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
682: std::map<PETSC_MESH_TYPE::point_type,PETSC_MESH_TYPE::point_type> renumbering;
686: DMMeshGetStratumSize(dm,"depth",0,&num_vertices);
687: DMMeshGetStratumSize(dm,"height",0,&num_cells);
688: PetscObjectGetComm((PetscObject) dm,&comm);
689: MPI_Comm_rank(comm,&rank);
691: DMMeshGetMesh(dm,mesh);
692: renumbering = mesh->getRenumbering();
694: my_num_cells_global = 0;
695: my_num_vertices_global = 0;
696: PetscMalloc(num_vertices*sizeof(PetscInt),&vertices);
697: /*
698: Get the total number of cells and vertices from the mapping, and build array of global indices of the local vertices
699: (TO array for the scatter) from the iterator
700: */
702: for (FlexMesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
703: /*
704: r_iter->first is a point number in the sequential (global) mesh
705: r_iter->second is the matching local point number in the distributed (local) mesh
706: */
707: if (r_iter->second > num_cells - 1 && r_iter->first > my_num_vertices_global) {
708: my_num_vertices_global = r_iter->first;
709: }
710: if (r_iter->second < num_cells && r_iter->first > my_num_cells_global) {
711: my_num_cells_global = r_iter->first;
712: }
713: if (r_iter->second > num_cells - 1) {
714: vertices[r_iter->second - num_cells] = r_iter->first;
715: }
716: }
717: my_num_cells_global++;
718: MPI_Allreduce(&my_num_cells_global,&num_cells_global,1,MPIU_INT,MPI_MAX,comm);
719: my_num_vertices_global -= num_cells_global-1;
720: MPI_Allreduce(&my_num_vertices_global,&num_vertices_global,1,MPIU_INT,MPI_MAX,comm);
722: /*
723: convert back vertices from point number to vertex numbers
724: */
725: for (i = 0; i < num_vertices; i++) {
726: vertices[i] -= num_cells_global;
727: }
729: /*
730: Build the IS and Vec required to create the VecScatter
731: A MUCH better way would be to use VecScatterCreateLocal and not create then destroy
732: Vecs, but I would need to understand better how it works
733: */
734: VecCreate(comm,&v_local);
735: VecSetSizes(v_local,num_vertices,PETSC_DETERMINE);
736: VecSetFromOptions(v_local);
737: VecCreate(comm,&v_zero);
738: VecSetSizes(v_zero,num_vertices_global*(rank==0),PETSC_DECIDE);
739: VecSetFromOptions(v_zero);
741: ISCreateGeneral(comm,num_vertices,vertices,PETSC_OWN_POINTER,&is_local);
742: VecScatterCreate(v_local,NULL,v_zero,is_local,scatter);
744: ISDestroy(&is_local);
745: VecDestroy(&v_local);
746: VecDestroy(&v_zero);
747: return(0);
748: }
752: /*@
753: DMMeshCreateScatterToZeroVertexSet - Creates the scatter required to scatter local (ghosted)
754: Vecs associated with a field defined at the a vertex set of a mesh to a Vec on cpu 0.
756: Input parameters:
757: . dm - the DMMesh representing the mesh
759: Output parameters:
760: . scatter - the scatter
762: Level: advanced
764: .keywords: mesh,ExodusII
765: .seealso DMMeshCreateScatterToZeroCell DMMeshCreateScatterToZeroVertex DMMeshCreateScatterToZeroCellSet
766: @*/
767: PetscErrorCode DMMeshCreateScatterToZeroVertexSet(DM dm,IS is_local,IS is_zero,VecScatter *scatter)
768: {
769: PetscErrorCode ierr;
770: const PetscInt *setvertices_local;
771: PetscInt *allvertices,*setvertices_zero,*setvertices_localtozero;
772: PetscInt num_vertices,num_vertices_global,my_num_vertices_global=0;
773: PetscInt num_cells,num_cells_global,my_num_cells_global=0;
774: PetscInt setsize_local,setsize_zero;
775: Vec v_local,v_zero;
776: MPI_Comm comm;
777: int rank,i,j,k,l,istart;
778: IS is_localtozero;
779: ALE::Obj<PETSC_MESH_TYPE> mesh;
780: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
781: std::map<PETSC_MESH_TYPE::point_type,PETSC_MESH_TYPE::point_type> renumbering;
785: DMMeshGetStratumSize(dm,"height",0,&num_cells);
786: DMMeshGetStratumSize(dm,"depth",0,&num_vertices);
787: PetscObjectGetComm((PetscObject) dm,&comm);
788: MPI_Comm_rank(comm,&rank);
790: DMMeshGetMesh(dm,mesh);
791: renumbering = mesh->getRenumbering();
793: PetscMalloc(num_vertices*sizeof(PetscInt),&allvertices);
794: /*
795: Build array of global indices of the local vertices (TO array for the scatter)
796: from the iterator
797: */
799: for (FlexMesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
800: /*
801: r_iter->first is a point number in the sequential (global) mesh
802: r_iter->second is the matching local point number in the distributed (local) mesh
803: */
804: if (r_iter->second > num_cells - 1 && r_iter->first > my_num_vertices_global) {
805: my_num_vertices_global = r_iter->first;
806: }
807: if (r_iter->second < num_cells && r_iter->first > my_num_cells_global) {
808: my_num_cells_global = r_iter->first;
809: }
810: if (r_iter->second > num_cells-1) {
811: allvertices[r_iter->second - num_cells] = r_iter->first;
812: }
813: }
814: my_num_cells_global++;
815: MPI_Allreduce(&my_num_cells_global,&num_cells_global,1,MPIU_INT,MPI_MAX,comm);
816: my_num_vertices_global -= num_cells_global-1;
817: MPI_Allreduce(&my_num_vertices_global,&num_vertices_global,1,MPIU_INT,MPI_MAX,comm);
819: ISGetSize(is_local,&setsize_local);
820: ISGetSize(is_zero,&setsize_zero);
821: PetscMalloc(setsize_local*sizeof(PetscInt),&setvertices_localtozero);
823: if (!rank) {
824: ISGetIndices(is_zero,(const PetscInt**)&setvertices_zero);
825: } else {
826: PetscMalloc(setsize_zero*sizeof(PetscInt),&setvertices_zero);
827: }
828: MPI_Bcast(setvertices_zero,setsize_zero,MPIU_INT,0,comm);
829: ISGetIndices(is_local,&setvertices_local);
832: istart = 0;
833: for (i = 0; i < setsize_local; i++) {
834: j = allvertices[setvertices_local[i]-num_cells];
835: /*
836: j is the cell number in the seq mesh of the i-th vertex of the local vertex set.
837: Search for the matching vertex in vertex set. Because vertex in the zero set are usually
838: numbered in ascending order, we start our search from the previous find.
839: ` */
841: for (l = 0, k = istart; l < setsize_zero; l++,k = (l+istart)%setsize_zero) {
842: if (setvertices_zero[k] == j) {
843: break;
844: }
845: }
846: setvertices_localtozero[i] = k;
847: istart = (k+1)%setsize_zero;
848: }
849: ISRestoreIndices(is_local,&setvertices_local);
850: if (!rank) {
851: ISRestoreIndices(is_zero,(const PetscInt**)&setvertices_zero);
852: } else {
853: PetscFree(setvertices_zero);
854: }
855: /*
856: Build the IS and Vec required to create the VecScatter
857: A MUCH better way would be to use VecScatterCreateLocal and not create then destroy
858: Vecs, but I would need to understand better how it works
859: */
860: VecCreate(comm,&v_local);
861: VecSetSizes(v_local,setsize_local,PETSC_DETERMINE);
862: VecSetFromOptions(v_local);
864: VecCreate(comm,&v_zero);
865: VecSetSizes(v_zero,setsize_zero*(rank==0),PETSC_DECIDE);
866: VecSetFromOptions(v_zero);
868: ISCreateGeneral(comm,setsize_local,setvertices_localtozero,PETSC_OWN_POINTER,&is_localtozero);
869: VecScatterCreate(v_local,NULL,v_zero,is_localtozero,scatter);
871: ISDestroy(&is_localtozero);
872: PetscFree(allvertices);
873: VecDestroy(&v_local);
874: VecDestroy(&v_zero);
875: return(0);
876: }
880: /*@
881: DMMeshCreateScatterToZeroCell - Creates the scatter required to scatter local (ghosted)
882: Vecs associated with a field defined at the cells of a mesh to a Vec on cpu 0.
884: Input parameters:
885: . dm - the DMMesh representing the mesh
887: Output parameters:
888: . scatter - the scatter
890: Level: advanced
892: .keywords: mesh,ExodusII
893: .seealso DMMeshCreateScatterToZeroVertex DMMeshCreateScatterToZeroVertexSet DMMeshCreateScatterToZeroCellSet
894: @*/
895: PetscErrorCode DMMeshCreateScatterToZeroCell(DM dm,VecScatter *scatter)
896: {
897: PetscErrorCode ierr;
898: PetscInt *cells;
899: PetscInt num_cells,num_cells_global,my_num_cells_global;
900: Vec v_local,v_zero;
901: MPI_Comm comm;
902: int rank;
903: IS is_local;
904: ALE::Obj<PETSC_MESH_TYPE> mesh;
905: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
906: std::map<PETSC_MESH_TYPE::point_type,PETSC_MESH_TYPE::point_type> renumbering;
910: DMMeshGetStratumSize(dm,"height",0,&num_cells);
911: PetscObjectGetComm((PetscObject) dm,&comm);
912: MPI_Comm_rank(comm,&rank);
914: DMMeshGetMesh(dm,mesh);
915: renumbering = mesh->getRenumbering();
917: my_num_cells_global = 0;
918: PetscMalloc(num_cells*sizeof(PetscInt),&cells);
919: /*
920: Get the total number of cells from the mapping, and build array of global indices of the local cells (TO array for the scatter)
921: from the iterator
922: */
924: for (FlexMesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
925: /*
926: r_iter->first is a point number in the sequential (global) mesh
927: r_iter->second is the matching local point number in the distributed (local) mesh
928: */
929: if (r_iter->second < num_cells && r_iter->first > my_num_cells_global) {
930: my_num_cells_global = r_iter->first;
931: }
932: if (r_iter->second < num_cells) {
933: cells[r_iter->second] = r_iter->first;
934: }
935: }
936: my_num_cells_global++;
937: MPI_Allreduce(&my_num_cells_global,&num_cells_global,1,MPIU_INT,MPI_MAX,comm);
939: /*
940: Build the IS and Vec required to create the VecScatter
941: A MUCH better way would be to use VecScatterCreateLocal and not create then destroy
942: Vecs, but I would need to understand better how it works
943: */
944: VecCreate(comm,&v_local);
945: VecSetSizes(v_local,num_cells,PETSC_DETERMINE);
946: VecSetFromOptions(v_local);
947: VecCreate(comm,&v_zero);
948: VecSetSizes(v_zero,num_cells_global*(rank==0),PETSC_DECIDE);
949: VecSetFromOptions(v_zero);
951: ISCreateGeneral(comm,num_cells,cells,PETSC_OWN_POINTER,&is_local);
952: VecScatterCreate(v_local,NULL,v_zero,is_local,scatter);
954: ISDestroy(&is_local);
955: VecDestroy(&v_local);
956: VecDestroy(&v_zero);
957: return(0);
958: }
962: /*@
963: DMMeshCreateScatterToZeroCellSet - Creates the scatter required to scatter local (ghosted)
964: Vecs associated with a field defined at a cell set of a mesh to a Vec on cpu 0.
966: Input parameters:
967: . dm - the DMMesh representing the mesh
969: Output parameters:
970: . scatter - the scatter
972: Level: advanced
974: .keywords: mesh,ExodusII
975: .seealso DMMeshCreateScatterToZeroCell DMMeshCreateScatterToZeroVertexSet DMMeshCreateScatterToZeroVertex
976: @*/
977: PetscErrorCode DMMeshCreateScatterToZeroCellSet(DM dm,IS is_local,IS is_zero,VecScatter *scatter)
978: {
979: PetscErrorCode ierr;
980: const PetscInt *setcells_local;
981: PetscInt *allcells,*setcells_zero,*setcells_localtozero;
982: PetscInt num_cells;
983: PetscInt setsize_local,setsize_zero;
984: Vec v_local,v_zero;
985: MPI_Comm comm;
986: int rank,i,j,k,l,istart;
987: IS is_localtozero;
988: ALE::Obj<PETSC_MESH_TYPE> mesh;
989: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
990: std::map<PETSC_MESH_TYPE::point_type,PETSC_MESH_TYPE::point_type> renumbering;
994: DMMeshGetStratumSize(dm,"height",0,&num_cells);
995: PetscObjectGetComm((PetscObject) dm,&comm);
996: MPI_Comm_rank(comm,&rank);
998: DMMeshGetMesh(dm,mesh);
999: renumbering = mesh->getRenumbering();
1001: PetscMalloc(num_cells*sizeof(PetscInt),&allcells);
1002: /*
1003: Build array of global indices of the local cells (TO array for the scatter)
1004: from the iterator
1005: */
1007: for (FlexMesh::renumbering_type::const_iterator r_iter = renumbering.begin(); r_iter != renumbering.end(); ++r_iter) {
1008: /*
1009: r_iter->first is a point number in the sequential (global) mesh
1010: r_iter->second is the matching local point number in the distributed (local) mesh
1011: */
1012: if (r_iter->second < num_cells) {
1013: allcells[r_iter->second] = r_iter->first;
1014: }
1015: }
1017: ISGetSize(is_local,&setsize_local);
1018: ISGetSize(is_zero,&setsize_zero);
1019: PetscMalloc(setsize_local*sizeof(PetscInt),&setcells_localtozero);
1021: if (!rank) {
1022: ISGetIndices(is_zero,(const PetscInt**)&setcells_zero);
1023: } else {
1024: PetscMalloc(setsize_zero*sizeof(PetscInt),&setcells_zero);
1025: }
1026: MPI_Bcast(setcells_zero,setsize_zero,MPIU_INT,0,comm);
1027: ISGetIndices(is_local,&setcells_local);
1029: istart = 0;
1030: for (i = 0; i < setsize_local; i++) {
1031: j = allcells[setcells_local[i]];
1032: /*
1033: j is the cell number in the seq mesh of the i-th cell of the local cell set.
1034: Search for the matching cell in cell set. Because cells in the zero set are usually
1035: numbered sequentially, we start our search from the previous find.
1036: ` */
1038: for (l = 0, k = istart; l < setsize_zero; l++,k = (l+istart)%setsize_zero) {
1039: if (setcells_zero[k] == j) {
1040: break;
1041: }
1042: }
1043: setcells_localtozero[i] = k;
1044: istart = (k+1)%setsize_zero;
1045: }
1046: ISRestoreIndices(is_local,&setcells_local);
1047: if (!rank) {
1048: ISRestoreIndices(is_zero,(const PetscInt**)&setcells_zero);
1049: } else {
1050: PetscFree(setcells_zero);
1051: }
1052: /*
1053: Build the IS and Vec required to create the VecScatter
1054: A MUCH better way would be to use VecScatterCreateLocal and not create then destroy
1055: Vecs, but I would need to understand better how it works
1056: */
1057: VecCreate(comm,&v_local);
1058: VecSetSizes(v_local,setsize_local,PETSC_DETERMINE);
1059: VecSetFromOptions(v_local);
1061: VecCreate(comm,&v_zero);
1062: VecSetSizes(v_zero,setsize_zero*(rank==0),PETSC_DECIDE);
1063: VecSetFromOptions(v_zero);
1065: ISCreateGeneral(comm,setsize_local,setcells_localtozero,PETSC_OWN_POINTER,&is_localtozero);
1066: VecScatterCreate(v_local,NULL,v_zero,is_localtozero,scatter);
1068: PetscFree(setcells_localtozero);
1069: PetscFree(allcells);
1070: VecDestroy(&v_local);
1071: VecDestroy(&v_zero);
1072: return(0);
1073: }
1077: /*@
1079: VecViewExodusVertex - Write a Vec representing nodal values of some field in an exodusII file.
1081: Collective on comm
1083: The behavior differs depending on the size of comm:
1084: if size(comm) == 1 each processor writes its local vector into a separate file
1085: if size(comm) > 1, the values are sent to cpu0, and written in a single file
1088: Input Parameters:
1089: + dm - the DMMesh representing the mesh
1090: . v - the LOCAL vector of values to be saved (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1091: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize()
1092: . exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1093: . step - the time step to write
1094: - exofield - the position in the exodus field of the first component
1096: Interpolated meshes are not supported.
1098: Level: beginner
1100: .keywords: mesh,ExodusII
1101: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecLoadExodusVertex()
1102: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCellSet() VecLoadExodusCellSet()
1103: VecViewExodusCell() VecLoadExodusCell()
1104: @*/
1105: PetscErrorCode VecViewExodusVertex(DM dm,Vec v,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1106: {
1107: #if defined(PETSC_HAVE_EXODUSII)
1108: PetscInt rank,num_proc;
1110: PetscInt c,num_vertices,num_vertices_zero,num_dof;
1111: Vec vdof,vdof_zero;
1112: PetscReal *vdof_array;
1113: VecScatter scatter;
1114: #endif
1117: #if defined(PETSC_HAVE_EXODUSII)
1118: MPI_Comm_size(comm,&num_proc);
1119: MPI_Comm_rank(comm,&rank);
1120: VecGetBlockSize(v,&num_dof);
1122: DMMeshGetStratumSize(dm,"depth",0,&num_vertices);
1124: VecCreate(comm,&vdof);
1125: VecSetSizes(vdof,num_vertices,PETSC_DETERMINE);
1126: VecSetFromOptions(vdof);
1128: if (num_proc == 1) {
1129: PetscMalloc(num_vertices*sizeof(PetscReal),&vdof_array);
1130: for (c = 0; c < num_dof; c++) {
1131: VecStrideGather(v,c,vdof,INSERT_VALUES);
1132: VecGetArray(vdof,&vdof_array);
1133: ex_put_nodal_var(exoid,step,exofield+c,num_vertices,vdof_array);
1134: VecRestoreArray(vdof,&vdof_array);
1135: if (ierr != 0) SETERRQ2(comm,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_nodal_var returned %i",exoid,ierr);
1136: }
1137: PetscFree(vdof_array);
1138: } else {
1139: /*
1140: Build the scatter sending one dof towards cpu 0
1141: */
1142: DMMeshCreateScatterToZeroVertex(dm,&scatter);
1143: /*
1144: Another an ugly hack to get the total number of vertices in the mesh. This has got to stop...
1145: */
1146: num_vertices_zero = scatter->to_n;
1148: VecCreate(comm,&vdof_zero);
1149: VecSetSizes(vdof_zero,num_vertices_zero,PETSC_DECIDE);
1150: VecSetFromOptions(vdof_zero);
1152: for (c = 0; c < num_dof; c++) {
1153: VecStrideGather(v,c,vdof,INSERT_VALUES);
1154: VecScatterBegin(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1155: VecScatterEnd(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1156: if (!rank) {
1157: VecGetArray(vdof_zero,&vdof_array);
1158: ex_put_nodal_var(exoid,step,exofield+c,num_vertices_zero,vdof_array);
1159: VecRestoreArray(vdof_zero,&vdof_array);
1160: }
1161: }
1162: /*
1163: Clean up
1164: */
1165: VecScatterDestroy(&scatter);
1166: VecDestroy(&vdof_zero);
1167: }
1168: VecDestroy(&vdof);
1169: #else
1170: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1171: #endif
1172: return(0);
1173: }
1177: /*@
1179: VecLoadExodusVertex - Loads a Vec representing nodal values of some field from an exodusII file.
1181: Collective on comm
1183: The behavior differs depending on the size of comm:
1184: if size(comm) == 1 each processor reads its local vector from a separate file
1185: if size(comm) > 1, the values are read by cpu 0 from a single file then scattered to each processor
1188: Input Parameters:
1189: + dm - the DMMesh representing the mesh
1190: . v - the LOCAL vector of values to be read (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1191: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize()
1192: . exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1193: . step - the time step to read
1194: - exofield - the position in the exodus field of the first component
1196: Interpolated meshes are not supported.
1198: Level: beginner
1200: .keywords: mesh,ExodusII
1201: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecViewExodusVertex()
1202: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCellSet() VecLoadExodusCellSet()
1203: VecViewExodusCell() VecLoadExodusCell()
1205: @*/
1206: PetscErrorCode VecLoadExodusVertex(DM dm,Vec v,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1207: {
1208: #if defined(PETSC_HAVE_EXODUSII)
1209: PetscInt rank,num_proc;
1211: PetscInt c,num_vertices,num_vertices_global,num_dof;
1212: Vec vdof,vzero;
1213: PetscReal *vdof_array;
1214: VecScatter scatter;
1215: #endif
1217: #if defined(PETSC_HAVE_EXODUSII)
1219: MPI_Comm_size(comm,&num_proc);
1220: MPI_Comm_rank(comm,&rank);
1221: VecGetBlockSize(v,&num_dof);
1223: DMMeshGetStratumSize(dm,"depth",0,&num_vertices);
1225: VecCreate(comm,&vdof);
1226: VecSetSizes(vdof,num_vertices,PETSC_DETERMINE);
1227: VecSetFromOptions(vdof);
1229: if (num_proc == 1) {
1230: for (c = 0; c < num_dof; c++) {
1231: VecGetArray(vdof,&vdof_array);
1232: ex_get_nodal_var(exoid,step,exofield+c,num_vertices,vdof_array);
1233: VecRestoreArray(vdof,&vdof_array);
1234: if (ierr != 0) SETERRQ2(comm,PETSC_ERR_FILE_READ,"Unable to read file id %i. ex_put_nodal_var returned %i",exoid,ierr);
1235: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1236: }
1237: } else {
1238: /*
1239: Build the scatter sending one dof towards cpu 0
1240: */
1241: DMMeshCreateScatterToZeroVertex(dm,&scatter);
1242: /*
1243: Another an ugly hack to get the total number of vertices in the mesh. This has got to stop...
1244: */
1245: num_vertices_global = scatter->to_n;
1246: MPI_Bcast(&num_vertices_global,1,MPIU_INT,0,comm);
1248: VecCreate(comm,&vzero);
1249: VecSetSizes(vzero,num_vertices_global*(rank==0),PETSC_DECIDE);
1250: VecSetFromOptions(vzero);
1252: for (c = 0; c < num_dof; c++) {
1253: if (!rank) {
1254: VecGetArray(vzero,&vdof_array);
1255: ex_get_nodal_var(exoid,step,exofield+c,num_vertices_global,vdof_array);
1256: VecRestoreArray(vzero,&vdof_array);
1257: }
1258: VecScatterBegin(scatter,vzero,vdof,INSERT_VALUES,SCATTER_REVERSE);
1259: VecScatterEnd(scatter,vzero,vdof,INSERT_VALUES,SCATTER_REVERSE);
1260: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1261: }
1262: /*
1263: Clean up
1264: */
1265: VecScatterDestroy(&scatter);
1266: VecDestroy(&vzero);
1267: }
1268: VecDestroy(&vdof);
1269: #else
1270: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1271: #endif
1272: return(0);
1273: }
1277: /*@
1279: VecViewExodusVertexSet - Write a Vec representing nodal values of some field at a vertex set in an exodusII file.
1281: Collective on comm
1283: The behavior differs depending on the size of comm:
1284: if size(comm) == 1 each processor writes its local vector into a separate file
1285: if size(comm) > 1, the values are sent to cpu0, and written in a single file
1288: Input Parameters:
1289: + dm - the DMMesh representing the mesh
1290: . v - the LOCAL vector of values to be saved (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1291: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize()
1292: . vsID - the vertex set ID
1293: . exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1294: . step - the time step to write
1295: - exofield - the position in the exodus field of the first component
1297: Interpolated meshes are not supported.
1299: Level: beginner
1301: .keywords: mesh,ExodusII
1302: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecLoadExodusVertexSet()
1303: VecViewExodusVertex() VecLoadExodusVertex() VecViewExodusCellSet() VecLoadExodusCellSet()
1304: VecViewExodusCell() VecLoadExodusCell()
1306: @*/
1307: PetscErrorCode VecViewExodusVertexSet(DM dm,Vec v,PetscInt vsID,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1308: {
1309: #if defined(PETSC_HAVE_EXODUSII)
1310: PetscInt rank,num_proc;
1312: PetscInt c,num_vertices_in_set=0,num_vertices_zero=0,num_cells_zero,num_dof;
1313: Vec vdof,vdof_zero;
1314: PetscReal *vdof_array;
1315: VecScatter scatter;
1316: PetscInt i,junk1,junk2,junk3,junk4,junk5;
1317: char junk[MAX_LINE_LENGTH+1];
1318: IS vsIS,vsIS_zero;
1319: PetscInt *vs_vertices_zero;
1320: MPI_Comm dmcomm;
1321: #endif
1323: #if defined(PETSC_HAVE_EXODUSII)
1325: MPI_Comm_size(comm,&num_proc);
1326: MPI_Comm_rank(comm,&rank);
1327: VecGetBlockSize(v,&num_dof);
1328: PetscObjectGetComm((PetscObject) dm,&dmcomm);
1330: DMMeshGetStratumSize(dm,"Vertex Sets",vsID,&num_vertices_in_set);
1331: VecCreateSeq(PETSC_COMM_SELF,num_vertices_in_set,&vdof);
1333: if (num_proc == 1) {
1334: if (num_vertices_in_set > 0) {
1335: for (c = 0; c < num_dof; c++) {
1336: VecStrideGather(v,c,vdof,INSERT_VALUES);
1337: VecGetArray(vdof,&vdof_array);
1338: ex_put_nset_var(exoid,step,exofield+c,vsID,num_vertices_in_set,vdof_array);
1339: if (ierr != 0) SETERRQ2(comm,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_nset_var returned %i",exoid,ierr);
1340: VecRestoreArray(vdof,&vdof_array);
1341: }
1342: }
1343: } else {
1344: /*
1345: Build the scatter sending one dof towards cpu 0
1346: */
1347: DMMeshGetStratumIS(dm,"Vertex Sets",vsID,&vsIS);
1348: if (!rank) {
1349: ex_get_node_set_param(exoid,vsID,&num_vertices_zero,&junk1);
1350: }
1351: PetscMalloc(num_vertices_zero*sizeof(PetscInt),&vs_vertices_zero);
1352: if (!rank) {
1353: ex_get_node_set(exoid,vsID,vs_vertices_zero);
1354: ex_get_init(exoid,junk,&junk1,&junk2,&num_cells_zero,&junk3,&junk4,&junk5);
1355: for (i = 0; i < num_vertices_zero; i++) {
1356: vs_vertices_zero[i] += num_cells_zero-1;
1357: }
1358: }
1359: ISCreateGeneral(dmcomm,num_vertices_zero,vs_vertices_zero,PETSC_OWN_POINTER,&vsIS_zero);
1360: DMMeshCreateScatterToZeroVertexSet(dm,vsIS,vsIS_zero,&scatter);
1361: ISDestroy(&vsIS);
1362: ISDestroy(&vsIS_zero);
1364: VecCreateSeq(PETSC_COMM_SELF,num_vertices_zero*(rank == 0),&vdof_zero);
1366: for (c = 0; c < num_dof; c++) {
1367: VecStrideGather(v,c,vdof,INSERT_VALUES);
1368: VecScatterBegin(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1369: VecScatterEnd(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1371: if (!rank) {
1372: VecGetArray(vdof_zero,&vdof_array);
1373: ex_put_nset_var(exoid,step,exofield+c,vsID,num_vertices_zero,vdof_array);
1374: if (ierr != 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_nset_var returned %i",exoid,ierr);
1375: VecRestoreArray(vdof_zero,&vdof_array);
1376: }
1377: }
1378: /*
1379: Clean up
1380: */
1381: VecScatterDestroy(&scatter);
1382: VecDestroy(&vdof_zero);
1383: }
1384: VecDestroy(&vdof);
1385: #else
1386: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1387: #endif
1388: return(0);
1389: }
1393: /*@
1395: VecLoadExodusVertexSet - Write a Vec representing nodal values of some field at a vertex set in an exodusII file.
1397: Collective on comm
1399: The behavior differs depending on the size of comm:
1400: if size(comm) == 1 each processor writes its local vector into a separate file
1401: if size(comm) > 1, the values are sent to cpu0, and written in a single file
1404: Input Parameters:
1405: + dm - the DMMesh representing the mesh
1406: . v - the LOCAL vector of values to be read (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1407: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize()
1408: . vsID - the vertex set ID
1409: . exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1410: . step - the time step to write
1411: - exofield - the position in the exodus field of the first component
1413: Interpolated meshes are not supported.
1415: Level: beginner
1417: .keywords: mesh,ExodusII
1418: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecLoadExodusVertex()
1420: @*/
1421: PetscErrorCode VecLoadExodusVertexSet(DM dm,Vec v,PetscInt vsID,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1422: {
1423: #if defined(PETSC_HAVE_EXODUSII)
1424: PetscInt rank,num_proc;
1426: PetscInt c,num_vertices_in_set=0,num_vertices_zero=0,num_cells_zero,num_dof;
1427: Vec vdof,vdof_zero;
1428: PetscReal *vdof_array;
1429: VecScatter scatter;
1430: PetscInt i,junk1,junk2,junk3,junk4,junk5;
1431: char junk[MAX_LINE_LENGTH+1];
1432: IS vsIS,vsIS_zero;
1433: PetscInt *vs_vertices_zero;
1434: MPI_Comm dmcomm;
1435: #endif
1437: #if defined(PETSC_HAVE_EXODUSII)
1439: MPI_Comm_size(comm,&num_proc);
1440: MPI_Comm_rank(comm,&rank);
1441: VecGetBlockSize(v,&num_dof);
1442: PetscObjectGetComm((PetscObject) dm,&dmcomm);
1444: DMMeshGetStratumSize(dm,"Vertex Sets",vsID,&num_vertices_in_set);
1445: VecCreateSeq(PETSC_COMM_SELF,num_vertices_in_set,&vdof);
1447: if (num_proc == 1) {
1448: if (num_vertices_in_set > 0) {
1449: for (c = 0; c < num_dof; c++) {
1450: VecGetArray(vdof,&vdof_array);
1451: ex_get_nset_var(exoid,step,exofield+c,vsID,num_vertices_in_set,vdof_array);
1452: if (ierr != 0) SETERRQ2(comm,PETSC_ERR_FILE_READ,"Unable to read from file id %i. ex_put_nset_var returned %i",exoid,ierr);
1453: VecRestoreArray(vdof,&vdof_array);
1454: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1455: }
1456: }
1457: } else {
1458: /*
1459: Build the scatter sending one dof towards cpu 0
1460: */
1461: DMMeshGetStratumIS(dm,"Vertex Sets",vsID,&vsIS);
1462: if (!rank) {
1463: ex_get_node_set_param(exoid,vsID,&num_vertices_zero,&junk1);
1464: }
1465: PetscMalloc(num_vertices_zero*sizeof(PetscInt),&vs_vertices_zero);
1466: if (!rank) {
1467: ex_get_node_set(exoid,vsID,vs_vertices_zero);
1468: ex_get_init(exoid,junk,&junk1,&junk2,&num_cells_zero,&junk3,&junk4,&junk5);
1469: for (i = 0; i < num_vertices_zero; i++) {
1470: vs_vertices_zero[i] += num_cells_zero-1;
1471: }
1472: }
1473: ISCreateGeneral(dmcomm,num_vertices_zero,vs_vertices_zero,PETSC_OWN_POINTER,&vsIS_zero);
1474: DMMeshCreateScatterToZeroVertexSet(dm,vsIS,vsIS_zero,&scatter);
1475: ISDestroy(&vsIS);
1476: ISDestroy(&vsIS_zero);
1478: VecCreateSeq(PETSC_COMM_SELF,num_vertices_zero*(rank == 0),&vdof_zero);
1480: for (c = 0; c < num_dof; c++) {
1481: if (!rank) {
1482: VecGetArray(vdof_zero,&vdof_array);
1483: ex_get_nset_var(exoid,step,exofield+c,vsID,num_vertices_zero,vdof_array);
1484: if (ierr != 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read from file id %i. ex_get_nset_var returned %i",exoid,ierr);
1485: VecRestoreArray(vdof_zero,&vdof_array);
1486: }
1487: VecScatterBegin(scatter,vdof_zero,vdof,INSERT_VALUES,SCATTER_REVERSE);
1488: VecScatterEnd(scatter,vdof_zero,vdof,INSERT_VALUES,SCATTER_REVERSE);
1490: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1491: }
1492: /*
1493: Clean up
1494: */
1495: VecScatterDestroy(&scatter);
1496: VecDestroy(&vdof_zero);
1497: }
1498: VecDestroy(&vdof);
1499: #else
1500: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1501: #endif
1502: return(0);
1503: }
1507: /*@
1509: VecViewExodusCell - Write a Vec representing the values of a field at all cells to an exodusII file.
1511: Input Parameters:
1512: + dm - the DMMesh representing the mesh
1513: . v - the LOCAL vector of values to be saved (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1514: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize().
1515: . comm - the communicator associated to the exo file
1516: + if size(comm) == 1 each processor writes its local vector into a separate file
1517: . if size(comm) > 1, the values are sent to cpu 0, and written in a single file
1518: - exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1519: . step - the time step to write
1520: . exofield - the position in the exodus field of the first component
1522: - Interpolated meshes are not supported.
1524: Level: beginner
1526: .keywords: mesh,ExodusII
1527: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecLoadExodusCell()
1528: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCellSet() VecLoadExodusCellSet()
1529: VecViewExodusVertex() VecLoadExodusVertex()
1531: @*/
1532: PetscErrorCode VecViewExodusCell(DM dm,Vec v,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1533: {
1534: #if defined(PETSC_HAVE_EXODUSII)
1536: PetscInt rank,num_proc,num_dof,num_cells,num_cells_in_set,num_cells_zero=0;
1537: PetscInt num_cs_zero,num_cs,set,c,istart;
1538: PetscInt *setsID_zero;
1539: const PetscInt *setsID;
1540: IS setIS,csIS,setIS_zero;
1541: VecScatter setscatter,setscattertozero;
1542: Vec vdof,vdof_set,vdof_set_zero;
1543: PetscReal *vdof_set_array;
1544: PetscInt junk1,junk2,junk3,junk4,junk5;
1545: char junk[MAX_LINE_LENGTH+1];
1546: MPI_Comm dmcomm;
1547: #endif
1549: #if defined(PETSC_HAVE_EXODUSII)
1551: MPI_Comm_size(comm,&num_proc);
1552: MPI_Comm_rank(comm,&rank);
1553: VecGetBlockSize(v,&num_dof);
1554: PetscObjectGetComm((PetscObject) dm,&dmcomm);
1555: DMMeshGetLabelSize(dm,"Cell Sets",&num_cs);
1556: DMMeshGetStratumSize(dm,"height",0,&num_cells);
1558: VecCreateSeq(PETSC_COMM_SELF,num_cells,&vdof);
1560: if (num_proc == 1) {
1561: DMMeshGetLabelSize(dm,"Cell Sets",&num_cs);
1562: DMMeshGetLabelIdIS(dm,"Cell Sets",&csIS);
1563: ISGetIndices(csIS,&setsID);
1565: for (set = 0; set < num_cs; set++) {
1566: /*
1567: Get the IS for the current set, the Vec containing a single dof and create
1568: the scatter for the restriction to the current cell set
1569: */
1570: DMMeshGetStratumIS(dm,"Cell Sets",setsID[set],&setIS);
1571: DMMeshGetStratumSize(dm,"Cell Sets",setsID[set],&num_cells_in_set);
1572: VecCreateSeq(PETSC_COMM_SELF,num_cells_in_set,&vdof_set);
1573: VecScatterCreate(vdof,setIS,vdof_set,NULL,&setscatter);
1574: for (c = 0; c < num_dof; c++) {
1575: VecStrideGather(v,c,vdof,INSERT_VALUES);
1577: /*
1578: Get the restriction of vdof to the current set
1579: */
1580: VecScatterBegin(setscatter,vdof,vdof_set,INSERT_VALUES,SCATTER_FORWARD);
1581: VecScatterEnd(setscatter,vdof,vdof_set,INSERT_VALUES,SCATTER_FORWARD);
1583: /*
1584: Write array to disk
1585: */
1586: VecGetArray(vdof_set,&vdof_set_array);
1587: ex_put_elem_var(exoid,step,exofield+c,setsID[set],num_cells_in_set,vdof_set_array);
1588: if (ierr != 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_elem_var returned %i",exoid,ierr);
1589: VecRestoreArray(vdof_set,&vdof_set_array);
1590: }
1591: VecDestroy(&vdof_set);
1592: VecScatterDestroy(&setscatter);
1593: ISDestroy(&setIS);
1594: }
1595: ISRestoreIndices(csIS,&setsID);
1596: ISDestroy(&csIS);
1597: } else {
1598: /*
1599: Get the number of blocks and the list of ID directly from the file. This is easier
1600: than trying to reconstruct the global lists from all individual IS
1601: */
1602: if (!rank) {
1603: ex_get_init(exoid,junk,&junk1,&junk2,&junk3,&num_cs_zero,&junk4,&junk5);
1604: }
1605: MPI_Bcast(&num_cs_zero,1,MPIU_INT,0,dmcomm);
1606: PetscMalloc(num_cs_zero * sizeof(PetscInt),&setsID_zero);
1607: if (!rank) {
1608: ex_get_elem_blk_ids(exoid,setsID_zero);
1609: }
1610: MPI_Bcast(setsID_zero,num_cs_zero,MPIU_INT,0,dmcomm);
1612: istart = 0;
1613: for (set = 0; set < num_cs_zero; set++) {
1614: /*
1615: Get the size of the size of the set on cpu 0 and create a Vec to receive values
1616: */
1617: if (!rank) {
1618: ex_get_elem_block(exoid,setsID_zero[set],junk,&num_cells_zero,&junk1,&junk2);
1619: }
1620: VecCreateSeq(PETSC_COMM_SELF,num_cells_zero,&vdof_set_zero);
1622: /*
1623: Get the IS for the current set, the Vec containing a single dof and create
1624: the scatter for the restriction to the current cell set
1625: */
1626: DMMeshGetStratumIS(dm,"Cell Sets",setsID_zero[set],&setIS);
1627: DMMeshGetStratumSize(dm,"Cell Sets",setsID_zero[set],&num_cells_in_set);
1628: VecCreateSeq(PETSC_COMM_SELF,num_cells_in_set,&vdof_set);
1629: VecScatterCreate(vdof,setIS,vdof_set,NULL,&setscatter);
1630: /*
1631: Get the scatter to send the values of a single dof at a single block to cpu 0
1632: */
1633: ISCreateStride(dmcomm,num_cells_zero,istart,1,&setIS_zero);
1634: DMMeshCreateScatterToZeroCellSet(dm,setIS,setIS_zero,&setscattertozero);
1635: ISDestroy(&setIS_zero);
1637: for (c = 0; c < num_dof; c++) {
1638: VecStrideGather(v,c,vdof,INSERT_VALUES);
1640: /*
1641: Get the restriction of vdof to the current set
1642: */
1643: VecScatterBegin(setscatter,vdof,vdof_set,INSERT_VALUES,SCATTER_FORWARD);
1644: VecScatterEnd(setscatter,vdof,vdof_set,INSERT_VALUES,SCATTER_FORWARD);
1646: /*
1647: Scatter vdof_set to cpu 0
1648: */
1649: VecScatterBegin(setscattertozero,vdof_set,vdof_set_zero,INSERT_VALUES,SCATTER_FORWARD);
1650: VecScatterEnd(setscattertozero,vdof_set,vdof_set_zero,INSERT_VALUES,SCATTER_FORWARD);
1652: /*
1653: Write array to disk
1654: */
1655: if (!rank) {
1656: VecGetArray(vdof_set_zero,&vdof_set_array);
1657: ex_put_elem_var(exoid,step,exofield+c,setsID_zero[set],num_cells_zero,vdof_set_array);
1658: VecRestoreArray(vdof_set_zero,&vdof_set_array);
1659: }
1660: }
1661: istart += num_cells_zero;
1662: VecDestroy(&vdof_set_zero);
1663: VecDestroy(&vdof_set);
1664: /*
1665: Does this really need to be protected with this test?
1666: */
1667: if (num_cells_in_set > 0) {
1668: /*
1669: Does this really need to be protected with this test?
1670: */
1671: VecScatterDestroy(&setscatter);
1672: }
1673: VecScatterDestroy(&setscattertozero);
1674: ISDestroy(&setIS);
1675: }
1676: PetscFree(setsID_zero);
1677: }
1679: VecDestroy(&vdof);
1680: #else
1681: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1682: #endif
1683: return(0);
1684: }
1688: /*@
1690: VecLoadExodusCell - Read a Vec representing the values of a field at all cells from an exodusII file.
1692: Input Parameters:
1693: + dm - the DMMesh representing the mesh
1694: . v - the LOCAL vector of values to be read (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1695: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize().
1696: . comm - the communicator associated to the exo file
1697: + if size(comm) == 1 each processor writes its local vector into a separate file
1698: . if size(comm) > 1, the values are sent to cpu 0, and written in a single file
1699: - exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1700: . step - the time step to write
1701: . exofield - the position in the exodus field of the first component
1703: - Interpolated meshes are not supported.
1705: Level: beginner
1707: .keywords: mesh,ExodusII
1708: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecViewExodusCell()
1709: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCellSet() VecLoadExodusCellSet()
1710: VecViewExodusVertex() VecLoadExodusVertex()
1712: @*/
1713: PetscErrorCode VecLoadExodusCell(DM dm,Vec v,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1714: {
1715: #if defined(PETSC_HAVE_EXODUSII)
1717: PetscInt rank,num_proc,num_dof,num_cells,num_cells_in_set,num_cells_zero=0;
1718: PetscInt num_cs_zero,num_cs,set,c,istart;
1719: PetscInt *setsID_zero;
1720: const PetscInt *setsID;
1721: IS setIS,csIS,setIS_zero;
1722: VecScatter setscatter,setscattertozero;
1723: Vec vdof,vdof_set,vdof_set_zero;
1724: PetscReal *vdof_set_array;
1725: PetscInt junk1,junk2,junk3,junk4,junk5;
1726: char junk[MAX_LINE_LENGTH+1];
1727: MPI_Comm dmcomm;
1728: #endif
1730: #if defined(PETSC_HAVE_EXODUSII)
1732: MPI_Comm_size(comm,&num_proc);
1733: MPI_Comm_rank(comm,&rank);
1734: VecGetBlockSize(v,&num_dof);
1735: PetscObjectGetComm((PetscObject) dm,&dmcomm);
1736: DMMeshGetLabelSize(dm,"Cell Sets",&num_cs);
1737: DMMeshGetStratumSize(dm,"height",0,&num_cells);
1739: VecCreateSeq(PETSC_COMM_SELF,num_cells,&vdof);
1741: if (num_proc == 1) {
1742: DMMeshGetLabelSize(dm,"Cell Sets",&num_cs);
1743: DMMeshGetLabelIdIS(dm,"Cell Sets",&csIS);
1744: ISGetIndices(csIS,&setsID);
1746: for (c = 0; c < num_dof; c++) {
1747: for (set = 0; set < num_cs; set++) {
1748: /*
1749: Get the IS for the current set, the Vec containing a single dof and create
1750: the scatter for the restriction to the current cell set
1751: */
1752: DMMeshGetStratumIS(dm,"Cell Sets",setsID[set],&setIS);
1753: DMMeshGetStratumSize(dm,"Cell Sets",setsID[set],&num_cells_in_set);
1754: VecCreateSeq(PETSC_COMM_SELF,num_cells_in_set,&vdof_set);
1755: VecScatterCreate(vdof,setIS,vdof_set,NULL,&setscatter);
1756: /*
1757: Read array from disk
1758: */
1759: VecGetArray(vdof_set,&vdof_set_array);
1760: ex_get_elem_var(exoid,step,exofield+c,setsID[set],num_cells_in_set,vdof_set_array);
1761: VecRestoreArray(vdof_set,&vdof_set_array);
1762: /*
1763: Copy values to full dof vec
1764: */
1765: VecScatterBegin(setscatter,vdof_set,vdof,INSERT_VALUES,SCATTER_REVERSE);
1766: VecScatterEnd(setscatter,vdof_set,vdof,INSERT_VALUES,SCATTER_REVERSE);
1768: }
1769: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1770: VecDestroy(&vdof_set);
1771: VecScatterDestroy(&setscatter);
1772: ISDestroy(&setIS);
1773: }
1774: ISRestoreIndices(csIS,&setsID);
1775: ISDestroy(&csIS);
1776: } else {
1777: /*
1778: Get the number of blocks and the list of ID directly from the file. This is easier
1779: than trying to reconstruct the global lists from all individual IS
1780: */
1781: if (!rank) {
1782: ex_get_init(exoid,junk,&junk1,&junk2,&junk3,&num_cs_zero,&junk4,&junk5);
1783: }
1784: MPI_Bcast(&num_cs_zero,1,MPIU_INT,0,dmcomm);
1785: PetscMalloc(num_cs_zero * sizeof(PetscInt),&setsID_zero);
1786: if (!rank) {
1787: ex_get_elem_blk_ids(exoid,setsID_zero);
1788: }
1789: MPI_Bcast(setsID_zero,num_cs_zero,MPIU_INT,0,dmcomm);
1791: for (c = 0; c < num_dof; c++) {
1792: istart = 0;
1793: for (set = 0; set < num_cs_zero; set++) {
1794: /*
1795: Get the size of the size of the set on cpu 0 and create a Vec to receive values
1796: */
1797: ex_get_elem_block(exoid,setsID_zero[set],junk,&num_cells_zero,&junk1,&junk2);
1798: VecCreateSeq(PETSC_COMM_SELF,num_cells_zero,&vdof_set_zero);
1800: /*
1801: Get the IS for the current set, the Vec containing a single dof and create
1802: the scatter for the restriction to the current cell set
1803: */
1804: DMMeshGetStratumIS(dm,"Cell Sets",setsID_zero[set],&setIS);
1805: DMMeshGetStratumSize(dm,"Cell Sets",setsID_zero[set],&num_cells_in_set);
1806: VecCreateSeq(PETSC_COMM_SELF,num_cells_in_set,&vdof_set);
1807: VecScatterCreate(vdof,setIS,vdof_set,NULL,&setscatter);
1808: /*
1809: Get the scatter to send the values of a single dof at a single block to cpu 0
1810: */
1811: ISCreateStride(dmcomm,num_cells_zero,istart,1,&setIS_zero);
1812: DMMeshCreateScatterToZeroCellSet(dm,setIS,setIS_zero,&setscattertozero);
1813: ISDestroy(&setIS_zero);
1815: /*
1816: Read array from disk
1817: */
1818: if (!rank) {
1819: VecGetArray(vdof_set_zero,&vdof_set_array);
1820: ex_get_elem_var(exoid,step,exofield+c,setsID_zero[set],num_cells_zero,vdof_set_array);
1821: VecRestoreArray(vdof_set_zero,&vdof_set_array);
1822: }
1823: /*
1824: Send values back to their owning cpu
1825: */
1826: VecScatterBegin(setscattertozero,vdof_set_zero,vdof_set,INSERT_VALUES,SCATTER_REVERSE);
1827: VecScatterEnd(setscattertozero,vdof_set_zero,vdof_set,INSERT_VALUES,SCATTER_REVERSE);
1828: /*
1829: Copy the values associated to the current set back into the component vec
1830: */
1831: VecScatterBegin(setscatter,vdof_set,vdof,INSERT_VALUES,SCATTER_REVERSE);
1832: VecScatterEnd(setscatter,vdof_set,vdof,INSERT_VALUES,SCATTER_REVERSE);
1833: istart += num_cells_zero;
1834: MPI_Bcast(&istart,1,MPIU_INT,0,dmcomm);
1835: }
1836: /*
1837: Copy the component back into v
1838: */
1839: VecStrideScatter(vdof,c,v,INSERT_VALUES);
1840: VecDestroy(&vdof_set_zero);
1841: VecDestroy(&vdof_set);
1842: VecScatterDestroy(&setscatter);
1843: VecScatterDestroy(&setscattertozero);
1844: ISDestroy(&setIS);
1845: }
1846: PetscFree(setsID_zero);
1847: VecDestroy(&vdof);
1848: }
1850: VecDestroy(&vdof);
1851: #else
1852: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1853: #endif
1854: return(0);
1855: }
1859: /*@
1861: VecViewExodusCellSet - Write a Vec representing the values of a field at a cell set to an exodusII file.
1863: Input Parameters:
1864: + dm - the DMMesh representing the mesh
1865: . v - the LOCAL vector of values to be saved (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1866: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize().
1867: . csID - the ID of the cell set
1868: . comm - the communicator associated to the exo file
1869: + if size(comm) == 1 each processor writes its local vector into a separate file
1870: . if size(comm) > 1, the values are sent to cpu 0, and written in a single file
1871: - exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1872: . step - the time step to write
1873: . exofield - the position in the exodus field of the first component
1875: - Interpolated meshes are not supported.
1877: Level: beginner
1879: .keywords: mesh,ExodusII
1880: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecLoadExodusCellSet()
1881: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCell() VecLoadExodusCell()
1882: VecViewExodusVertex() VecLoadExodusVertex()
1884: @*/
1885: PetscErrorCode VecViewExodusCellSet(DM dm,Vec v,PetscInt csID,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
1886: {
1887: #if defined(PETSC_HAVE_EXODUSII)
1889: PetscInt rank,num_proc,num_dof,num_cells,num_cells_zero=0;
1890: PetscInt i,c;
1891: PetscReal *vdof_array;
1892: IS csIS,csIS_zero;
1893: Vec vdof,vdof_zero;
1894: VecScatter scatter;
1895: PetscInt istart = 0;
1897: PetscInt junk1,junk2,junk3,junk4,junk5;
1898: char junk[MAX_LINE_LENGTH+1];
1899: PetscInt num_cs_global;
1900: PetscInt *csIDs;
1901: MPI_Comm dmcomm;
1902: #endif
1904: #if defined(PETSC_HAVE_EXODUSII)
1906: MPI_Comm_size(comm,&num_proc);
1907: MPI_Comm_rank(comm,&rank);
1908: VecGetBlockSize(v,&num_dof);
1909: PetscObjectGetComm((PetscObject) dm,&dmcomm);
1910: DMMeshGetStratumIS(dm,"Cell Sets",csID,&csIS);
1911: ISGetSize(csIS,&num_cells);
1913: VecCreateSeq(PETSC_COMM_SELF,num_cells,&vdof);
1915: if (num_proc == 1) {
1916: for (c = 0; c < num_dof; c++) {
1917: VecStrideGather(v,c,vdof,INSERT_VALUES);
1918: VecGetArray(vdof,&vdof_array);
1919: ex_put_elem_var(exoid,step,exofield+c,csID,num_cells,vdof_array);
1920: VecRestoreArray(vdof,&vdof_array);
1921: if (ierr != 0) SETERRQ2(comm,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_elem_var returned %i",exoid,ierr);
1922: }
1923: } else {
1924: /*
1925: Get the global size of the cell set from the file because we can
1926: There is no direct way to get the index of the first cell in a cell set from exo.
1927: (which depends on the order on the file) so we need to seek the cell set in the file and
1928: accumulate offsets...
1929: */
1930: if (!rank) {
1931: istart = 0;
1932: ex_get_init(exoid,junk,&junk1,&junk2,&junk3,&num_cs_global,&junk4,&junk5);
1933: PetscMalloc(num_cs_global * sizeof(PetscInt),&csIDs);
1934: ex_get_elem_blk_ids(exoid,csIDs);
1935: for (i = 0; i < num_cs_global; i++) {
1936: ex_get_elem_block(exoid,csIDs[i],junk,&num_cells_zero,&junk1,&junk2);
1937: if (csIDs[i] == csID) break;
1938: else {
1939: istart += num_cells_zero;
1940: }
1941: }
1942: PetscFree(csIDs);
1943: ex_get_elem_block(exoid,csID,junk,&num_cells_zero,&junk1,&junk2);
1944: }
1945: ISCreateStride(dmcomm,num_cells_zero,istart,1,&csIS_zero);
1946: DMMeshCreateScatterToZeroCellSet(dm,csIS,csIS_zero,&scatter);
1947: ISDestroy(&csIS_zero);
1949: VecCreateSeq(PETSC_COMM_SELF,num_cells_zero,&vdof_zero);
1951: for (c = 0; c < num_dof; c++) {
1952: VecStrideGather(v,c,vdof,INSERT_VALUES);
1953: VecScatterBegin(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1954: VecScatterEnd(scatter,vdof,vdof_zero,INSERT_VALUES,SCATTER_FORWARD);
1955: if (!rank) {
1956: VecGetArray(vdof_zero,&vdof_array);
1957: ex_put_elem_var(exoid,step,exofield+c,csID,num_cells_zero,vdof_array);
1958: VecRestoreArray(vdof_zero,&vdof_array);
1959: }
1960: }
1961: VecDestroy(&vdof_zero);
1962: VecScatterDestroy(&scatter);
1963: }
1964: ISDestroy(&csIS);
1965: VecDestroy(&vdof);
1966: #else
1967: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
1968: #endif
1969: return(0);
1970: }
1974: /*@
1976: VecLoadExodusCellSet - Read a Vec representing the values of a field at a cell set from an exodusII file.
1978: Input Parameters:
1979: + dm - the DMMesh representing the mesh
1980: . v - the LOCAL vector of values to be saved (i.e. with ghost values) obtained with SectionRealCreateLocalVector.
1981: if v represents a field with several components, the block size must be set accordingly using VecSetBlockSize().
1982: . csID - the ID of the cell set
1983: . comm - the communicator associated to the exo file
1984: + if size(comm) == 1 each processor writes its local vector into a separate file
1985: . if size(comm) > 1, the values are sent to cpu 0, and written in a single file
1986: - exoid - the id of the exodusII file (obtained with ex_open or ex_create)
1987: . step - the time step to write
1988: . exofield - the position in the exodus field of the first component
1990: - Interpolated meshes are not supported.
1992: Level: beginner
1994: .keywords: mesh,ExodusII
1995: .seealso: MeshCreate() MeshCreateExodus() SectionRealCreateLocalVector() VecViewExodusCellSet()
1996: VecViewExodusVertexSet() VecLoadExodusVertexSet() VecViewExodusCell() VecLoadExodusCell()
1997: VecViewExodusVertex() VecLoadExodusVertex()
1999: @*/
2000: PetscErrorCode VecLoadExodusCellSet(DM dm,Vec v,PetscInt csID,MPI_Comm comm,PetscInt exoid,PetscInt step,PetscInt exofield)
2001: {
2002: #if defined(PETSC_HAVE_EXODUSII)
2004: PetscInt rank,num_proc,num_dof,num_cells,num_cells_zero=0;
2005: PetscInt i,c;
2006: PetscReal *vdof_array;
2007: IS csIS,csIS_zero;
2008: Vec vdof,vdof_zero;
2009: VecScatter scatter;
2010: PetscInt istart = 0;
2012: PetscInt junk1,junk2,junk3,junk4,junk5;
2013: PetscInt num_cs_global;
2014: PetscInt *csIDs;
2015: char junk[MAX_LINE_LENGTH+1];
2016: MPI_Comm dmcomm;
2017: #endif
2019: #if defined(PETSC_HAVE_EXODUSII)
2021: MPI_Comm_size(comm,&num_proc);
2022: MPI_Comm_rank(comm,&rank);
2023: VecGetBlockSize(v,&num_dof);
2024: PetscObjectGetComm((PetscObject) dm,&dmcomm);
2026: DMMeshGetStratumIS(dm,"Cell Sets",csID,&csIS);
2027: ISGetSize(csIS,&num_cells);
2029: VecCreateSeq(PETSC_COMM_SELF,num_cells,&vdof);
2031: if (num_proc == 1) {
2032: for (c = 0; c < num_dof; c++) {
2033: VecGetArray(vdof,&vdof_array);
2034: ex_get_elem_var(exoid,step,exofield+c,csID,num_cells,vdof_array);
2035: VecRestoreArray(vdof,&vdof_array);
2036: if (ierr) SETERRQ2(comm,PETSC_ERR_FILE_WRITE,"Unable to write to file id %i. ex_put_elem_var returned %i",exoid,ierr);
2037: VecStrideScatter(vdof,c,v,INSERT_VALUES);
2038: }
2039: } else {
2040: /*
2041: Get the global size of the cell set from the file because we can
2042: There is no direct way to get the index of the first cell in a cell set from exo.
2043: (which depends on the order on the file) so we need to seek the cell set in the file and
2044: accumulate offsets...
2045: */
2046: if (!rank) {
2047: istart = 0;
2048: ex_get_init(exoid,junk,&junk1,&junk2,&junk3,&num_cs_global,&junk4,&junk5);
2049: PetscMalloc(num_cs_global * sizeof(PetscInt),&csIDs);
2050: ex_get_elem_blk_ids(exoid,csIDs);
2051: for (i = 0; i < num_cs_global; i++) {
2052: ex_get_elem_block(exoid,csIDs[i],junk,&num_cells_zero,&junk1,&junk2);
2053: if (csIDs[i] == csID) break;
2054: else istart += num_cells_zero;
2055: }
2056: PetscFree(csIDs);
2057: ex_get_elem_block(exoid,csID,junk,&num_cells_zero,&junk1,&junk2);
2058: }
2059: ISCreateStride(dmcomm,num_cells_zero,istart,1,&csIS_zero);
2060: DMMeshCreateScatterToZeroCellSet(dm,csIS,csIS_zero,&scatter);
2061: ISDestroy(&csIS_zero);
2063: VecCreateSeq(PETSC_COMM_SELF,num_cells_zero,&vdof_zero);
2065: for (c = 0; c < num_dof; c++) {
2066: if (!rank) {
2067: VecGetArray(vdof_zero,&vdof_array);
2068: ex_get_elem_var(exoid,step,exofield+c,csID,num_cells_zero,vdof_array);
2069: VecRestoreArray(vdof_zero,&vdof_array);
2070: }
2071: VecScatterBegin(scatter,vdof_zero,vdof,INSERT_VALUES,SCATTER_REVERSE);
2072: VecScatterEnd(scatter,vdof_zero,vdof,INSERT_VALUES,SCATTER_REVERSE);
2073: VecStrideScatter(vdof,c,v,INSERT_VALUES);
2074: }
2075: VecDestroy(&vdof_zero);
2076: VecScatterDestroy(&scatter);
2077: }
2078: ISDestroy(&csIS);
2079: VecDestroy(&vdof);
2080: #else
2081: SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "This method requires ExodusII support. Reconfigure using --download-exodusii");
2082: #endif
2083: return(0);
2084: }