Actual source code: meshpcice.c
petsc-3.4.5 2014-06-29
1: #include <petscdmmesh_formats.hh> /*I "petscmesh.h" I*/
5: PetscErrorCode WritePCICEVertices(DM dm, PetscViewer viewer)
6: {
7: ALE::Obj<PETSC_MESH_TYPE> m;
8: PetscErrorCode ierr;
10: DMMeshGetMesh(dm, m);
11: return ALE::PCICE::Viewer::writeVertices(m, viewer);
12: }
16: PetscErrorCode WritePCICEElements(DM dm, PetscViewer viewer)
17: {
18: ALE::Obj<PETSC_MESH_TYPE> m;
19: PetscErrorCode ierr;
21: DMMeshGetMesh(dm, m);
22: return ALE::PCICE::Viewer::writeElements(m, viewer);
23: }
27: PetscErrorCode WritePCICERestart(DM dm, PetscViewer viewer)
28: {
29: ALE::Obj<PETSC_MESH_TYPE> m;
30: PetscErrorCode ierr;
32: DMMeshGetMesh(dm, m);
33: return ALE::PCICE::Viewer::writeRestart(m, viewer);
34: }
38: /*@C
39: DMMeshCreatePCICE - Create a DMMesh from PCICE files.
41: Not Collective
43: Input Parameters:
44: + dim - The topological mesh dimension
45: . coordFilename - The file containing vertex coordinates
46: . adjFilename - The file containing the vertices for each element
47: . interpolate - The flag for construction of intermediate elements
48: . bcFilename - The file containing the boundary topology and conditions
49: . numBdFaces - The number of boundary faces (or edges)
50: - numBdVertices - The number of boundary vertices
52: Output Parameter:
53: . mesh - The DMMesh object
55: Level: beginner
57: .keywords: mesh, PCICE
58: .seealso: DMMeshCreate()
59: @*/
60: PetscErrorCode DMMeshCreatePCICE(MPI_Comm comm, const int dim, const char coordFilename[], const char adjFilename[], PetscBool interpolate, const char bcFilename[], DM *mesh)
61: {
62: ALE::Obj<PETSC_MESH_TYPE> m;
63: PetscInt debug = 0;
64: PetscBool flag;
65: PetscErrorCode ierr;
68: DMMeshCreate(comm, mesh);
69: PetscOptionsGetInt(NULL, "-debug", &debug, &flag);
70: try {
71: m = ALE::PCICE::Builder::readMesh(comm, dim, std::string(coordFilename), std::string(adjFilename), true, interpolate, debug);
72: if (debug) m->view("Mesh");
73: } catch(ALE::Exception e) {
74: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN, e.message());
75: }
76: if (bcFilename) ALE::PCICE::Builder::readBoundary(m, std::string(bcFilename));
77: DMMeshSetMesh(*mesh, m);
78: return(0);
79: }
83: /*@C
84: PCICERenumberBoundary - Change global element names into offsets
86: Collective on DMMesh
88: Input Parameters:
89: . mesh - the mesh
91: Level: advanced
93: .seealso: DMMeshCreate()
94: @*/
95: PetscErrorCode PCICERenumberBoundary(DM dm)
96: {
97: ALE::Obj<PETSC_MESH_TYPE> m;
98: PetscErrorCode ierr;
101: DMMeshGetMesh(dm, m);
102: try {
103: ALE::PCICE::fuseBoundary(m);
104: } catch(ALE::Exception e) {
105: SETERRQ(PETSC_COMM_SELF,100, e.msg().c_str());
106: }
107: return(0);
108: }
112: /*@C
113: BCSectionGetArray - Returns the array underlying the BCSection.
115: Not Collective
117: Input Parameters:
118: + mesh - The DMMesh object
119: - name - The section name
121: Output Parameters:
122: + numElements - The number of mesh element with values
123: . fiberDim - The number of values per element
124: - array - The array
126: Level: intermediate
128: .keywords: mesh, elements
129: .seealso: DMMeshCreate()
130: @*/
131: PetscErrorCode BCSectionGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscInt *array[])
132: {
133: ALE::Obj<PETSC_MESH_TYPE> m;
134: PetscErrorCode ierr;
137: DMMeshGetMesh(dm, m);
138: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& section = m->getIntSection(std::string(name));
139: if (!section->size()) {
140: *numElements = 0;
141: *fiberDim = 0;
142: *array = NULL;
143: return(0);
144: }
145: const PETSC_MESH_TYPE::int_section_type::chart_type& chart = section->getChart();
147: int fiberDimMin = section->getFiberDimension(*chart.begin());
148: int numElem = 0;
150: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
151: const int fiberDim = section->getFiberDimension(*c_iter);
153: if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
154: }
155: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
156: const int fiberDim = section->getFiberDimension(*c_iter);
158: numElem += fiberDim/fiberDimMin;
159: }
160: *numElements = numElem;
161: *fiberDim = fiberDimMin;
162: *array = (PetscInt*) section->restrictSpace();
163: return(0);
164: }
168: /*@C
169: BCSectionRealCreate - Creates a BCSection.
171: Not Collective
173: Input Parameters:
174: + mesh - The DMMesh object
175: . name - The section name
176: - fiberDim - The number of values per element
178: Level: intermediate
180: .keywords: mesh, elements
181: .seealso: DMMeshCreate()
182: @*/
183: PetscErrorCode BCSectionRealCreate(DM dm, const char name[], PetscInt fiberDim)
184: {
185: ALE::Obj<PETSC_MESH_TYPE> m;
186: PetscErrorCode ierr;
189: DMMeshGetMesh(dm, m);
190: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
191: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& ibc = m->getIntSection("IBC");
192: const PETSC_MESH_TYPE::int_section_type::chart_type& chart = ibc->getChart();
194: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
195: section->setFiberDimension(*p_iter, ibc->getFiberDimension(*p_iter));
196: }
197: m->allocate(section);
198: return(0);
199: }
203: /*@C
204: BCSectionRealGetArray - Returns the array underlying the BCSection.
206: Not Collective
208: Input Parameters:
209: + mesh - The DMMesh object
210: - name - The section name
212: Output Parameters:
213: + numElements - The number of mesh element with values
214: . fiberDim - The number of values per element
215: - array - The array
217: Level: intermediate
219: .keywords: mesh, elements
220: .seealso: DMMeshCreate()
221: @*/
222: PetscErrorCode BCSectionRealGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscReal *array[])
223: {
224: ALE::Obj<PETSC_MESH_TYPE> m;
225: PetscErrorCode ierr;
228: DMMeshGetMesh(dm, m);
229: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
230: if (!section->size()) {
231: *numElements = 0;
232: *fiberDim = 0;
233: *array = NULL;
234: return(0);
235: }
236: const PETSC_MESH_TYPE::real_section_type::chart_type& chart = section->getChart();
238: int fiberDimMin = section->getFiberDimension(*chart.begin());
239: int numElem = 0;
241: for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
242: const int fiberDim = section->getFiberDimension(*c_iter);
244: if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
245: }
246: for (PETSC_MESH_TYPE::real_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
247: const int fiberDim = section->getFiberDimension(*c_iter);
249: numElem += fiberDim/fiberDimMin;
250: }
251: *numElements = numElem;
252: *fiberDim = fiberDimMin;
253: *array = (PetscReal*) section->restrictSpace();
254: return(0);
255: }
259: PetscErrorCode BCFUNCGetArray(DM dm, PetscInt *numElements, PetscInt *fiberDim, PetscScalar *array[])
260: {
261: ALE::Obj<PETSC_MESH_TYPE> m;
262: PetscErrorCode ierr;
265: DMMeshGetMesh(dm, m);
266: #if 0
267: PETSC_MESH_TYPE::bc_values_type& bcValues = m->getBCValues();
268: *numElements = bcValues.size();
269: *fiberDim = 4;
270: *array = new PetscScalar[(*numElements)*(*fiberDim)];
271: for (int bcf = 1; bcf <= (int) bcValues.size(); ++bcf) {
272: (*array)[(bcf-1)*4+0] = bcValues[bcf].rho;
273: (*array)[(bcf-1)*4+1] = bcValues[bcf].u;
274: (*array)[(bcf-1)*4+2] = bcValues[bcf].v;
275: (*array)[(bcf-1)*4+3] = bcValues[bcf].p;
276: }
277: #else
278: *numElements = 0;
279: *fiberDim = 0;
280: *array = NULL;
281: #endif
282: return(0);
283: }
285: namespace ALE {
286: namespace PCICE {
287: //
288: // Builder methods
289: //
290: void Builder::readConnectivity(MPI_Comm comm, const std::string& filename, int& corners, const bool useZeroBase, int& numElements, int *vertices[])
291: {
292: PetscViewer viewer;
293: FILE *f;
294: PetscInt numCells, cellCount = 0;
295: PetscInt *verts;
296: char buf[2048];
297: PetscInt c;
298: PetscInt commRank;
301: MPI_Comm_rank(comm, &commRank);
303: if (commRank != 0) return;
304: PetscViewerCreate(PETSC_COMM_SELF, &viewer);
305: PetscViewerSetType(viewer, PETSCVIEWERASCII);
306: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
307: PetscViewerFileSetName(viewer, filename.c_str());
308: if (ierr) {
309: ostringstream txt;
310: txt << "Could not open PCICE connectivity file: " << filename;
311: throw ALE::Exception(txt.str().c_str());
312: }
313: PetscViewerASCIIGetPointer(viewer, &f);
314: if (!fgets(buf, 2048, f)) {
315: throw ALE::Exception("Invalid connectivity file: Missing number of elements");
316: }
317: const char *sizes = strtok(buf, " ");
318: numCells = atoi(sizes);
319: sizes = strtok(NULL, " ");
320: if (sizes != NULL) {
321: corners = atoi(sizes);
322: std::cout << "Reset corners to " << corners << std::endl;
323: }
324: PetscMalloc(numCells*corners * sizeof(PetscInt), &verts);
325: while (fgets(buf, 2048, f) != NULL) {
326: const char *v = strtok(buf, " ");
328: /* Ignore cell number */
329: v = strtok(NULL, " ");
330: for (c = 0; c < corners; c++) {
331: int vertex = atoi(v);
333: if (!useZeroBase) vertex -= 1;
334: verts[cellCount*corners+c] = vertex;
336: v = strtok(NULL, " ");
337: }
338: cellCount++;
339: }
340: PetscViewerDestroy(&viewer);
341: numElements = numCells;
342: *vertices = verts;
343: };
344: void Builder::readCoordinates(MPI_Comm comm, const std::string& filename, const int dim, int& numVertices, PetscReal *coordinates[])
345: {
346: PetscViewer viewer;
347: FILE *f;
348: PetscInt numVerts, vertexCount = 0;
349: PetscReal *coords;
350: char buf[2048];
351: PetscInt c;
352: PetscInt commRank;
355: MPI_Comm_rank(comm, &commRank);
357: if (commRank != 0) return;
358: PetscViewerCreate(PETSC_COMM_SELF, &viewer);
359: PetscViewerSetType(viewer, PETSCVIEWERASCII);
360: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
361: PetscViewerFileSetName(viewer, filename.c_str());
362: if (ierr) {
363: ostringstream txt;
364: txt << "Could not open PCICE coordinate file: " << filename;
365: throw ALE::Exception(txt.str().c_str());
366: }
367: PetscViewerASCIIGetPointer(viewer, &f);
368: numVerts = atoi(fgets(buf, 2048, f));
369: PetscMalloc(numVerts*dim * sizeof(PetscReal), &coords);
370: while (fgets(buf, 2048, f) != NULL) {
371: const char *x = strtok(buf, " ");
373: /* Ignore vertex number */
374: x = strtok(NULL, " ");
375: for (c = 0; c < dim; c++) {
376: coords[vertexCount*dim+c] = atof(x);
378: x = strtok(NULL, " ");
379: }
380: vertexCount++;
381: }
382: PetscViewerDestroy(&viewer);
383: numVertices = numVerts;
384: *coordinates = coords;
385: };
386: Obj<PETSC_MESH_TYPE> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& basename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0)
387: {
388: return readMesh(comm, dim, basename+".nodes", basename+".lcon", useZeroBase, interpolate, debug);
389: };
390: #if defined(PETSC_OPT_SIEVE)
391: Obj<PETSC_MESH_TYPE> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& coordFilename, const std::string& adjFilename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0)
392: {
393: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
394: Obj<Mesh> mesh = new Mesh(comm, dim, debug);
395: Obj<sieve_type> sieve = new sieve_type(comm, debug);
396: const Obj<FlexMesh> m = new FlexMesh(comm, dim, debug);
397: const Obj<FlexMesh::sieve_type> s = new FlexMesh::sieve_type(comm, debug);
398: int *cells = NULL;
399: PetscReal *coordinates = NULL;
400: int numCells = 0, numVertices = 0, numCorners = dim+1;
401: PetscErrorCode ierr;
403: ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
404: ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
405: ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
406: m->setSieve(s);
407: m->stratify();
408: mesh->setSieve(sieve);
409: std::map<Mesh::point_type,Mesh::point_type> renumbering;
410: ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering, false);
411: mesh->stratify();
412: ALE::ISieveConverter::convertOrientation(*s, *sieve, renumbering, m->getArrowSection("orientation").ptr());
413: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
414: if (cells) {PetscFree(cells);CHKERRXX(ierr);}
415: if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
416: return mesh;
417: };
418: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename)
419: {
420: throw ALE::Exception("Not implemented for optimized sieves");
421: };
422: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, PetscReal *coordinates[], const bool columnMajor)
423: {
424: const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
425: if (!coordSec->size()) {
426: *numVertices = 0;
427: *dim = 0;
428: *coordinates = NULL;
429: return;
430: }
431: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
432: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
433: int size = vertices->size();
434: int embedDim = coordSec->getFiberDimension(*vertices->begin());
435: PetscReal *coords;
436: PetscErrorCode ierr;
438: PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
439: for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
440: const Mesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);
441: const int row = vNumbering->getIndex(*v_iter);
443: if (columnMajor) {
444: for (int d = 0; d < embedDim; d++) {
445: coords[d*size + row] = array[d];
446: }
447: } else {
448: for (int d = 0; d < embedDim; d++) {
449: coords[row*embedDim + d] = array[d];
450: }
451: }
452: }
453: *numVertices = size;
454: *dim = embedDim;
455: *coordinates = coords;
456: };
457: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor)
458: {
459: if (!mesh->heightStratum(0)->size()) {
460: *numElements = 0;
461: *numCorners = 0;
462: *vertices = NULL;
463: return;
464: }
465: const Obj<Mesh::sieve_type>& sieve = mesh->getSieve();
466: const Obj<Mesh::label_sequence>& elements = mesh->heightStratum(0);
467: const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
468: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
469: int size = elements->size();
470: int corners = sieve->getConeSize(*elements->begin());
471: int *v;
472: PetscErrorCode ierr;
474: PetscMalloc(size*corners * sizeof(int), &v);CHKERRXX(ierr);
475: for (Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
476: const Obj<Mesh::sieve_type::coneSequence> cone = sieve->cone(*e_iter);
477: Mesh::sieve_type::coneSequence::const_iterator begin = cone->begin();
478: Mesh::sieve_type::coneSequence::const_iterator end = cone->end();
480: const int row = eNumbering->getIndex(*e_iter);
481: int c = -1;
482: if (columnMajor) {
483: for (Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
484: v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
485: }
486: } else {
487: for (Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
488: v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
489: }
490: }
491: }
492: *numElements = size;
493: *numCorners = corners;
494: *vertices = v;
495: };
496: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
497: {
498: throw ALE::Exception("Not implemented for optimized sieves");
499: };
500: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
501: {
502: throw ALE::Exception("Not implemented for optimized sieves");
503: };
504: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer)
505: {
506: throw ALE::Exception("Not implemented for optimized sieves");
507: };
508: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer)
509: {
510: throw ALE::Exception("Not implemented for optimized sieves");
511: };
512: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh)
513: {
514: throw ALE::Exception("Not implemented for optimized sieves");
515: };
516: #else
517: Obj<PETSC_MESH_TYPE> Builder::readMesh(MPI_Comm comm, const int dim, const std::string& coordFilename, const std::string& adjFilename, const bool useZeroBase = true, const bool interpolate = true, const int debug = 0)
518: {
519: Obj<Mesh> mesh = new DMMesh(comm, dim, debug);
520: Obj<sieve_type> sieve = new sieve_type(comm, debug);
521: int *cells = NULL;
522: double *coordinates = NULL;
523: int numCells = 0, numVertices = 0, numCorners = dim+1;
524: PetscErrorCode ierr;
526: ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
527: ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
528: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildTopology(sieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
529: mesh->setSieve(sieve);
530: mesh->stratify();
531: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
532: if (cells) {PetscFree(cells);CHKERRXX(ierr);}
533: if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
534: return mesh;
535: };
536: // Creates boundary sections:
537: // IBC[NBFS,2]: ALL
538: // BL[NBFS,1]:
539: // BNVEC[NBFS,2]:
540: // BCFUNC[NBCF,NV]: ALL
541: // IBNDFS[NBN,2]: STILL NEED 4-5
542: // BNNV[NBN,2]
543: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename)
544: {
545: PetscViewer viewer;
546: FILE *f;
547: char buf[2048];
550: const Obj<Mesh::int_section_type>& ibc = mesh->getIntSection("IBC");
551: const Obj<Mesh::int_section_type>& ibndfs = mesh->getIntSection("IBNDFS");
552: const Obj<Mesh::int_section_type>& ibcnum = mesh->getIntSection("IBCNUM");
553: const Obj<Mesh::int_section_type>& ibfcon = mesh->getIntSection("IBFCON");
554: const Obj<Mesh::real_section_type>& bl = mesh->getRealSection("BL");
555: const Obj<Mesh::real_section_type>& bnvec = mesh->getRealSection("BNVEC");
556: const Obj<Mesh::real_section_type>& bnnv = mesh->getRealSection("BNNV");
557: if (mesh->commRank() != 0) {
558: #if 0
559: mesh->distributeBCValues();
560: #endif
561: return;
562: }
563: PetscViewerCreate(PETSC_COMM_SELF, &viewer);CHKERRXX(ierr);
564: PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRXX(ierr);
565: PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRXX(ierr);
566: PetscViewerFileSetName(viewer, bcFilename.c_str());CHKERRXX(ierr);
567: PetscViewerASCIIGetPointer(viewer, &f);CHKERRXX(ierr);
568: // Create IBC section
569: int numBdFaces = atoi(strtok(fgets(buf, 2048, f), " "));
570: int *tmpIBC = new int[numBdFaces*4];
571: std::map<int,std::set<int> > elem2Idx;
572: std::map<int,int> bfReorder;
573: for (int bf = 0; bf < numBdFaces; bf++) {
574: const char *x = strtok(fgets(buf, 2048, f), " ");
576: // Ignore boundary face number
577: x = strtok(NULL, " "); tmpIBC[bf*4+0] = atoi(x);
578: x = strtok(NULL, " "); tmpIBC[bf*4+1] = atoi(x);
579: x = strtok(NULL, " "); tmpIBC[bf*4+2] = atoi(x);
580: x = strtok(NULL, " "); tmpIBC[bf*4+3] = atoi(x);
581: const int elem = tmpIBC[bf*4+0]-1;
583: ibc->addFiberDimension(elem, 4);
584: ibcnum->addFiberDimension(elem, 1);
585: ibfcon->addFiberDimension(elem, 2);
586: bl->addFiberDimension(elem, 1);
587: bnvec->addFiberDimension(elem, 2);
588: elem2Idx[elem].insert(bf);
589: }
590: mesh->allocate(ibc);
591: mesh->allocate(ibcnum);
592: mesh->allocate(ibfcon);
593: mesh->allocate(bl);
594: mesh->allocate(bnvec);
595: const DMMesh::int_section_type::chart_type& chart = ibc->getChart();
597: int num = 1;
599: for (Mesh::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
600: const int elem = *p_iter;
601: int bfNum[2];
602: int k = 0;
604: for (std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
605: bfReorder[(*i_iter)+1] = num;
606: bfNum[k++] = num;
607: num++;
608: }
609: ibcnum->updatePoint(elem, bfNum);
610: }
611: for (int bf = 0; bf < numBdFaces; bf++) {
612: const int elem = tmpIBC[bf*4]-1;
614: if (elem2Idx[elem].size() > 1) {
615: if (*elem2Idx[elem].begin() == bf) {
616: int values[8];
617: int k = 0;
619: for (std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
620: for (int v = 0; v < 4; ++v) {
621: values[k*4+v] = tmpIBC[*i_iter*4+v];
622: }
623: k++;
624: }
625: ibc->updatePoint(elem, values);
626: }
627: } else {
628: ibc->updatePoint(elem, &tmpIBC[bf*4]);
629: }
630: }
631: delete [] tmpIBC;
632: // Create BCFUNC section
633: int numBcFunc = atoi(strtok(fgets(buf, 2048, f), " "));
634: if (numBcFunc != 0) throw ALE::Exception("Cannot handle BCFUNCS after rewrite");
635: for (int bc = 0; bc < numBcFunc; bc++) {
636: #if 0
637: const char *x = strtok(fgets(buf, 2048, f), " ");
638: DMMesh::bc_value_type value;
640: // Ignore function number
641: x = strtok(NULL, " "); value.rho = atof(x);
642: x = strtok(NULL, " "); value.u = atof(x);
643: x = strtok(NULL, " "); value.v = atof(x);
644: x = strtok(NULL, " "); value.p = atof(x);
645: mesh->setBCValue(bc+1, value);
646: #endif
647: }
648: #if 0
649: mesh->distributeBCValues();
650: #endif
651: // Create IBNDFS section
652: int numBdVertices = atoi(strtok(fgets(buf, 2048, f), " "));
653: const int numElements = mesh->heightStratum(0)->size();
654: int *tmpIBNDFS = new int[numBdVertices*3];
656: for (int bv = 0; bv < numBdVertices; bv++) {
657: const char *x = strtok(fgets(buf, 2048, f), " ");
659: // Ignore boundary node number
660: x = strtok(NULL, " "); tmpIBNDFS[bv*3+0] = atoi(x);
661: x = strtok(NULL, " "); tmpIBNDFS[bv*3+1] = atoi(x);
662: x = strtok(NULL, " "); tmpIBNDFS[bv*3+2] = atoi(x);
663: ibndfs->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, 6);
664: }
665: mesh->allocate(ibndfs);
666: for (int bv = 0; bv < numBdVertices; bv++) {
667: int values[5];
669: values[0] = tmpIBNDFS[bv*3+0];
670: // Covert to new boundary face numbers
671: values[1] = bfReorder[tmpIBNDFS[bv*3+1]];
672: values[2] = bfReorder[tmpIBNDFS[bv*3+2]];
673: values[3] = 0;
674: values[4] = 0;
675: ibndfs->updatePoint(values[0]-1+numElements, values);
676: }
677: PetscViewerDestroy(&viewer);CHKERRXX(ierr);
678: // Create BNNV[NBN,2]
679: const int dim = mesh->getDimension();
681: for (int bv = 0; bv < numBdVertices; bv++) {
682: bnnv->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, dim);
683: }
684: mesh->allocate(bnnv);
685: delete [] tmpIBNDFS;
686: };
687: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, double *coordinates[], const bool columnMajor)
688: {
689: const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
690: if (!coordSec->size()) {
691: *numVertices = 0;
692: *dim = 0;
693: *coordinates = NULL;
694: return;
695: }
696: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
697: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
699: int size = vertices->size();
700: int embedDim = coordSec->getFiberDimension(*vertices->begin());
701: double *coords;
704: PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
705: for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
706: const DMMesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);
708: const int row = vNumbering->getIndex(*v_iter);
710: if (columnMajor) {
711: for (int d = 0; d < embedDim; d++) coords[d*size + row] = array[d];
712: } else {
713: for (int d = 0; d < embedDim; d++) coords[row*embedDim + d] = array[d];
714: }
715: }
716: *numVertices = size;
717: *dim = embedDim;
718: *coordinates = coords;
719: };
720: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor)
721: {
722: if (!mesh->heightStratum(0)->size()) {
723: *numElements = 0;
724: *numCorners = 0;
725: *vertices = NULL;
726: return;
727: }
728: const Obj<Mesh::sieve_type>& sieve = mesh->getSieve();
729: const Obj<Mesh::label_sequence>& elements = mesh->heightStratum(0);
730: const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
731: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
733: int size = elements->size();
734: int corners = sieve->getConeSize(*elements->begin());
735: int *v;
738: PetscMalloc(elements->size()*corners * sizeof(int), &v);CHKERRXX(ierr);
739: for (Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
740: const Obj<Mesh::sieve_type::traits::coneSequence> cone = sieve->cone(*e_iter);
741: DMMesh::sieve_type::traits::coneSequence::iterator begin = cone->begin();
742: DMMesh::sieve_type::traits::coneSequence::iterator end = cone->end();
744: const int row = eNumbering->getIndex(*e_iter);
745: int c = -1;
746: if (columnMajor) {
747: for (Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
748: v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
749: }
750: } else {
751: for (Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
752: v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
753: }
754: }
755: }
756: *numElements = size;
757: *numCorners = corners;
758: *vertices = v;
759: };
763: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
764: {
765: ALE::Obj<Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
766: #if 0
767: DMMesh::field_type::patch_type patch;
768: const double *array = coordinates->restrict (patch);
769: int numVertices;
770: PetscErrorCode ierr;
773: //FIX:
774: if (vertexBundle->getGlobalOffsets()) {
775: numVertices = vertexBundle->getGlobalOffsets()[mesh->commSize()];
776: } else {
777: numVertices = mesh->getTopology()->depthStratum(0)->size();
778: }
779: PetscViewerASCIIPrintf(viewer, "%D\n", numVertices);
780: if (mesh->commRank() == 0) {
781: int numLocalVertices = mesh->getTopology()->depthStratum(0)->size();
782: int embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
783: int vertexCount = 1;
785: for (int v = 0; v < numLocalVertices; v++) {
786: PetscViewerASCIIPrintf(viewer, "%7D ", vertexCount++);
787: for (int d = 0; d < embedDim; d++) {
788: if (d > 0) {
789: PetscViewerASCIIPrintf(viewer, " ");
790: }
791: PetscViewerASCIIPrintf(viewer, "% 12.5E", array[v*embedDim+d]);
792: }
793: PetscViewerASCIIPrintf(viewer, "\n");
794: }
795: for (int p = 1; p < mesh->commSize(); p++) {
796: double *remoteCoords;
797: MPI_Status status;
799: MPI_Recv(&numLocalVertices, 1, MPI_INT, p, 1, mesh->comm(), &status);
800: PetscMalloc(numLocalVertices*embedDim * sizeof(double), &remoteCoords);
801: MPI_Recv(remoteCoords, numLocalVertices*embedDim, MPI_DOUBLE, p, 1, mesh->comm(), &status);
802: for (int v = 0; v < numLocalVertices; v++) {
803: PetscViewerASCIIPrintf(viewer,"%7D ", vertexCount++);
804: for (int d = 0; d < embedDim; d++) {
805: if (d > 0) {
806: PetscViewerASCIIPrintf(viewer, " ");
807: }
808: PetscViewerASCIIPrintf(viewer, "% 12.5E", remoteCoords[v*embedDim+d]);
809: }
810: PetscViewerASCIIPrintf(viewer, "\n");
811: }
812: }
813: } else {
814: ALE::Obj<Mesh::bundle_type> globalOrder = coordinates->getGlobalOrder();
815: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = globalOrder->getPatch(patch);
816: const int *offsets = coordinates->getGlobalOffsets();
817: int embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
818: int numLocalVertices = (offsets[mesh->commRank()+1] - offsets[mesh->commRank()])/embedDim;
819: double *localCoords;
820: int k = 0;
822: PetscMalloc(numLocalVertices*embedDim * sizeof(double), &localCoords);
823: for (Mesh::bundle_type::order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
824: int dim = globalOrder->getFiberDimension(patch, *p_iter);
826: if (dim > 0) {
827: int offset = coordinates->getFiberOffset(patch, *p_iter);
829: for (int i = offset; i < offset+dim; ++i) localCoords[k++] = array[i];
830: }
831: }
832: if (k != numLocalVertices*embedDim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of coordinates to send %d should be %d", k, numLocalVertices*embedDim);
833: MPI_Send(&numLocalVertices, 1, MPI_INT, 0, 1, mesh->comm());
834: MPI_Send(localCoords, numLocalVertices*embedDim, MPI_DOUBLE, 0, 1, mesh->comm());
835: PetscFree(localCoords);
836: }
837: #endif
838: return(0);
839: };
843: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer)
844: {
845: #if 0
846: ALE::Obj<Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
847: ALE::Obj<Mesh::bundle_type> elementBundle = mesh->getBundle(topology->depth());
848: ALE::Obj<Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
849: ALE::Obj<Mesh::bundle_type> globalVertex = vertexBundle->getGlobalOrder();
850: ALE::Obj<Mesh::bundle_type> globalElement = elementBundle->getGlobalOrder();
851: DMMesh::bundle_type::patch_type patch;
852: std::string orderName("element");
853: int dim = mesh->getDimension();
854: int corners = topology->nCone(*elements->begin(), topology->depth())->size();
855: int numElements;
859: if (corners != dim+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "PCICE only supports simplicies");
860: if (!globalVertex) globalVertex = vertexBundle;
862: if (elementBundle->getGlobalOffsets()) {
863: numElements = elementBundle->getGlobalOffsets()[mesh->commSize()];
864: } else {
865: numElements = mesh->getTopology()->heightStratum(0)->size();
866: }
867: if (mesh->commRank() == 0) {
868: int elementCount = 1;
870: PetscViewerASCIIPrintf(viewer, "%d\n", numElements);
871: for (Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
872: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);
874: PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
875: for (Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
876: PetscViewerASCIIPrintf(viewer, " %7d", globalVertex->getIndex(patch, *c_itor).prefix);
877: }
878: PetscViewerASCIIPrintf(viewer, "\n");
879: }
880: for (int p = 1; p < mesh->commSize(); p++) {
881: int numLocalElements;
882: int *remoteVertices;
883: MPI_Status status;
885: MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
886: PetscMalloc(numLocalElements*corners * sizeof(int), &remoteVertices);
887: MPI_Recv(remoteVertices, numLocalElements*corners, MPI_INT, p, 1, mesh->comm(), &status);
888: for (int e = 0; e < numLocalElements; e++) {
889: PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
890: for (int c = 0; c < corners; c++) {
891: PetscViewerASCIIPrintf(viewer, " %7d", remoteVertices[e*corners+c]);
892: }
893: PetscViewerASCIIPrintf(viewer, "\n");
894: }
895: PetscFree(remoteVertices);
896: }
897: } else {
898: const int *offsets = elementBundle->getGlobalOffsets();
899: int numLocalElements = offsets[mesh->commRank()+1] - offsets[mesh->commRank()];
900: int *localVertices;
901: int k = 0;
903: PetscMalloc(numLocalElements*corners * sizeof(int), &localVertices);
904: for (Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
905: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);
907: if (globalElement->getFiberDimension(patch, *e_itor) > 0) {
908: for (Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
909: localVertices[k++] = globalVertex->getIndex(patch, *c_itor).prefix;
910: }
911: }
912: }
913: if (k != numLocalElements*corners) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numLocalElements*corners);
914: MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
915: MPI_Send(localVertices, numLocalElements*corners, MPI_INT, 0, 1, mesh->comm());
916: PetscFree(localVertices);
917: }
918: #endif
919: return(0);
920: };
924: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer)
925: {
926: Obj<Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
927: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
928: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
929: int embedDim = coordinates->getFiberDimension(*vertices->begin());
930: PetscErrorCode ierr;
933: PetscViewerASCIIPrintf(viewer, "%D\n", vertices->size());
934: for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
935: const DMMesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
937: PetscViewerASCIIPrintf(viewer, "%7D ", vNumbering->getIndex(*v_iter)+1);
938: for (int d = 0; d < embedDim; d++) {
939: if (d > 0) PetscViewerASCIIPrintf(viewer, " ");
940: PetscViewerASCIIPrintf(viewer, "% 12.5E", array[d]);
941: }
942: PetscViewerASCIIPrintf(viewer, "\n");
943: }
944: return(0);
945: };
949: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer)
950: {
951: const Obj<Mesh::real_section_type>& velocity = mesh->getRealSection("VELN");
952: const Obj<Mesh::real_section_type>& pressure = mesh->getRealSection("PN");
953: const Obj<Mesh::real_section_type>& temperature = mesh->getRealSection("TN");
954: const Obj<Mesh::numbering_type>& cNumbering = mesh->getFactory()->getNumbering(mesh, mesh->depth());
955: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getNumbering(mesh, 0);
956: const int numCells = cNumbering->getGlobalSize();
957: PetscErrorCode ierr;
960: int blen[2];
961: MPI_Aint indices[2];
962: MPI_Datatype oldtypes[2], newtype;
963: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
964: blen[1] = 4; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
966: MPI_Type_create_struct(2, blen, indices, oldtypes, &newtype);
967: MPI_Type_commit(&newtype);
969: if (mesh->commRank() == 0) {
970: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
972: for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
973: if (vNumbering->isLocal(*v_iter)) {
974: const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
975: const DMMesh::real_section_type::value_type *pn = pressure->restrictPoint(*v_iter);
976: const DMMesh::real_section_type::value_type *tn = temperature->restrictPoint(*v_iter);
978: PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", *v_iter-numCells+1, veln[0], veln[1], pn[0], tn[0]);
979: }
980: }
981: for (int p = 1; p < mesh->commSize(); p++) {
982: RestartType *remoteValues;
983: int numLocalElements;
984: MPI_Status status;
986: MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
987: PetscMalloc(numLocalElements * sizeof(RestartType), &remoteValues);
988: MPI_Recv(remoteValues, numLocalElements, newtype, p, 1, mesh->comm(), &status);
989: for (int e = 0; e < numLocalElements; e++) {
990: PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", remoteValues[e].vertex-numCells+1, remoteValues[e].veln_x, remoteValues[e].veln_y, remoteValues[e].pn, remoteValues[e].tn);
991: }
992: }
993: } else {
994: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
995: RestartType *localValues;
996: int numLocalElements = vNumbering->getLocalSize();
997: int k = 0;
999: PetscMalloc(numLocalElements * sizeof(RestartType), &localValues);
1000: for (Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
1001: if (vNumbering->isLocal(*v_iter)) {
1002: const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
1003: const DMMesh::real_section_type::value_type *pn = pressure->restrictPoint(*v_iter);
1004: const DMMesh::real_section_type::value_type *tn = temperature->restrictPoint(*v_iter);
1006: localValues[k].vertex = *v_iter;
1007: localValues[k].veln_x = veln[0];
1008: localValues[k].veln_y = veln[1];
1009: localValues[k].pn = pn[0];
1010: localValues[k].tn = tn[0];
1011: k++;
1012: }
1013: }
1014: if (k != numLocalElements) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of values to send for field, %d should be %d", k, numLocalElements);
1015: MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
1016: MPI_Send(localValues, numLocalElements, newtype, 0, 1, mesh->comm());
1017: PetscFree(localValues);
1018: }
1019: MPI_Type_free(&newtype);
1020: return(0);
1021: };
1023: // This class reconstructs the local pieces of the boundary that distributed PCICE needs.
1024: // The boundary along with the boundary conditions is encoded in a collection of sections
1025: // over the PCICE mesh. These sections contain a description of the boundary topology
1026: // using elements' global names. This is unacceptable for PCICE, since it interprets
1027: // elements of the connectivity data arrays as local offsets into (some other of) these arrays.
1028: // This subroutine performs the renumbering based on the local numbering on the distributed
1029: // mesh. Essentially, we replace each global element id by its local number.
1030: //
1031: // Note: Any vertex or element number from PCICE is 1-based, but in Sieve we are 0-based. Thus
1032: // we add and subtract 1 during conversion. Also, Sieve vertices are numbered after cells.
1033: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh)
1034: {
1035: // Extract PCICE boundary sections
1036: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCsec = mesh->getIntSection("IBC");
1037: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBNDFSsec = mesh->getIntSection("IBNDFS");
1038: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCNUMsec = mesh->getIntSection("IBCNUM");
1040: // Look at the sections, if debugging
1041: if (mesh->debug()) {
1042: IBCsec->view("IBCsec", mesh->comm());
1043: IBNDFSsec->view("IBNDFSsec", mesh->comm());
1044: }
1046: // Extract the local numberings
1047: ALE::Obj<PETSC_MESH_TYPE::numbering_type> vertexNum = mesh->getFactory()->getLocalNumbering(mesh, 0);
1048: ALE::Obj<PETSC_MESH_TYPE::numbering_type> elementNum = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
1049: const int numElements = mesh->getFactory()->getNumbering(mesh, mesh->depth())->getGlobalSize();
1050: std::map<int,int> bfMap;
1051: // Declare points to the extracted numbering data
1052: const PETSC_MESH_TYPE::numbering_type::value_type *vNum, *eNum;
1054: // Create map from serial bdFace numbers to local bdFace numbers
1055: {
1056: const PETSC_MESH_TYPE::int_section_type::chart_type chart = IBCNUMsec->getChart();
1057: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin = chart.begin();
1058: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end = chart.end();
1060: int num = 0;
1061: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1062: const int fiberDim = IBCNUMsec->getFiberDimension(*p_iter);
1063: const int *globalNum = IBCNUMsec->restrictPoint(*p_iter);
1065: for (int n = 0; n < fiberDim; ++n) bfMap[globalNum[n]] = ++num;
1066: }
1067: }
1068: // Renumber vertex section IBC
1069: {
1070: const PETSC_MESH_TYPE::int_section_type::chart_type IBCchart = IBCsec->getChart();
1071: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin = IBCchart.begin();
1072: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end = IBCchart.end();
1073: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1074: PETSC_MESH_TYPE::point_type p = *p_iter;
1075: const PETSC_MESH_TYPE::int_section_type::value_type *ibc_in = IBCsec->restrictPoint(p);
1076: int fiberDimension = IBCsec->getFiberDimension(p);
1077: PETSC_MESH_TYPE::int_section_type::value_type ibc_out[8];
1078: // k controls the update of edge connectivity for each containing element;
1079: // if fiberDimension is 4, only one boundary face is connected to the element, and that edge's data
1080: // are contained in entries 0 - 3 of the section over the element p;
1081: // if fiberDimension is 8, two boundary faces are connected to the element, so the second edge's data
1082: // are contained in entries 4 - 7
1083: for (int k = 0; k < fiberDimension/4; k++) {
1084: // Extract IBC entry 1 (entry kk*4) for edge kk connected to element p.
1085: // This is the entry that needs renumbering for renumbering (2,3 & 4 are invariant under distribution),
1086: // see IBC's description.
1087: // Here we assume that elementNum's domain contains all boundary elements found in IBC,
1088: // so we can restrict to the extracted entry.
1089: eNum = elementNum->restrictPoint((PETSC_MESH_TYPE::numbering_type::point_type) ibc_in[k*4]-1);
1090: ibc_out[k*4+0] = eNum[0]+1;
1091: // Copy the other entries right over
1092: ibc_out[k*4+1] = ibc_in[k*4+1];
1093: ibc_out[k*4+2] = ibc_in[k*4+2];
1094: ibc_out[k*4+3] = ibc_in[k*4+3];
1095: }
1096: // Update IBC
1097: IBCsec->updatePoint(p, ibc_out);
1098: }
1099: }
1100: {
1101: // Renumber vertex section IBNDFS
1102: const PETSC_MESH_TYPE::int_section_type::chart_type IBNDFSchart = IBNDFSsec->getChart();
1103: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin = IBNDFSchart.begin();
1104: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end = IBNDFSchart.end();
1105: for (PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1106: PETSC_MESH_TYPE::point_type p = *p_iter;
1107: const PETSC_MESH_TYPE::int_section_type::value_type *ibndfs_in = IBNDFSsec->restrictPoint(p);
1108: // Here we assume the fiber dimension is 5
1109: PETSC_MESH_TYPE::int_section_type::value_type ibndfs_out[5];
1110: // Renumber entries 1,2 & 3 (4 & 5 are invariant under distribution), see IBNDFS's description
1111: // Here we assume that vertexNum's domain contains all boundary verticies found in IBNDFS, so we can restrict to the first extracted entry
1112: vNum = vertexNum->restrictPoint((PETSC_MESH_TYPE::numbering_type::point_type)ibndfs_in[0]-1+numElements);
1113: ibndfs_out[0] = vNum[0]+1;
1114: // Map serial bdFace numbers to local bdFace numbers
1115: ibndfs_out[1] = bfMap[ibndfs_in[1]];
1116: ibndfs_out[2] = bfMap[ibndfs_in[2]];
1117: // Copy the other entries right over
1118: ibndfs_out[3] = ibndfs_in[3];
1119: ibndfs_out[4] = ibndfs_in[4];
1120: // Update IBNDFS
1121: IBNDFSsec->updatePoint(p,ibndfs_out);
1122: }
1123: }
1124: if (mesh->debug()) {
1125: IBCsec->view("Renumbered IBCsec", mesh->comm());
1126: IBNDFSsec->view("Renumbered IBNDFSsec", mesh->comm());
1127: }
1128: // It's not clear whether IBFCON needs to be renumbered (what does it mean that its entries are not "GLOBAL NODE NUMBER(s)" -- see IBFCON's descriptions
1129: };
1130: #endif
1131: };
1132: };