Actual source code: meshpcice.c
petsc-3.3-p7 2013-05-11
1: #include <petscdmmesh_formats.hh> /*I "petscmesh.h" I*/
5: PetscErrorCode WritePCICEVertices(DM dm, PetscViewer viewer)
6: {
7: ALE::Obj<PETSC_MESH_TYPE> m;
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;
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;
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(PETSC_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) {
77: ALE::PCICE::Builder::readBoundary(m, std::string(bcFilename));
78: }
79: DMMeshSetMesh(*mesh, m);
80: return(0);
81: }
85: /*@C
86: PCICERenumberBoundary - Change global element names into offsets
88: Collective on DMMesh
90: Input Parameters:
91: . mesh - the mesh
93: Level: advanced
95: .seealso: DMMeshCreate()
96: @*/
97: PetscErrorCode PCICERenumberBoundary(DM dm)
98: {
99: ALE::Obj<PETSC_MESH_TYPE> m;
100: PetscErrorCode ierr;
103: DMMeshGetMesh(dm, m);
104: try {
105: ALE::PCICE::fuseBoundary(m);
106: } catch(ALE::Exception e) {
107: SETERRQ(PETSC_COMM_SELF,100, e.msg().c_str());
108: }
109: return(0);
110: }
114: /*@C
115: BCSectionGetArray - Returns the array underlying the BCSection.
117: Not Collective
119: Input Parameters:
120: + mesh - The DMMesh object
121: - name - The section name
123: Output Parameters:
124: + numElements - The number of mesh element with values
125: . fiberDim - The number of values per element
126: - array - The array
128: Level: intermediate
130: .keywords: mesh, elements
131: .seealso: DMMeshCreate()
132: @*/
133: PetscErrorCode BCSectionGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscInt *array[])
134: {
135: ALE::Obj<PETSC_MESH_TYPE> m;
136: PetscErrorCode ierr;
139: DMMeshGetMesh(dm, m);
140: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& section = m->getIntSection(std::string(name));
141: if (!section->size()) {
142: *numElements = 0;
143: *fiberDim = 0;
144: *array = NULL;
145: return(0);
146: }
147: const PETSC_MESH_TYPE::int_section_type::chart_type& chart = section->getChart();
148: int fiberDimMin = section->getFiberDimension(*chart.begin());
149: int numElem = 0;
151: for(PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
152: const int fiberDim = section->getFiberDimension(*c_iter);
154: if (fiberDim < fiberDimMin) fiberDimMin = fiberDim;
155: }
156: for(PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
157: const int fiberDim = section->getFiberDimension(*c_iter);
159: numElem += fiberDim/fiberDimMin;
160: }
161: *numElements = numElem;
162: *fiberDim = fiberDimMin;
163: *array = (PetscInt *) section->restrictSpace();
164: return(0);
165: }
169: /*@C
170: BCSectionRealCreate - Creates a BCSection.
172: Not Collective
174: Input Parameters:
175: + mesh - The DMMesh object
176: . name - The section name
177: - fiberDim - The number of values per element
179: Level: intermediate
181: .keywords: mesh, elements
182: .seealso: DMMeshCreate()
183: @*/
184: PetscErrorCode BCSectionRealCreate(DM dm, const char name[], PetscInt fiberDim)
185: {
186: ALE::Obj<PETSC_MESH_TYPE> m;
187: PetscErrorCode ierr;
190: DMMeshGetMesh(dm, m);
191: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
192: const ALE::Obj<PETSC_MESH_TYPE::int_section_type>& ibc = m->getIntSection("IBC");
193: const PETSC_MESH_TYPE::int_section_type::chart_type& chart = ibc->getChart();
195: for(PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
196: section->setFiberDimension(*p_iter, ibc->getFiberDimension(*p_iter));
197: }
198: m->allocate(section);
199: return(0);
200: }
204: /*@C
205: BCSectionRealGetArray - Returns the array underlying the BCSection.
207: Not Collective
209: Input Parameters:
210: + mesh - The DMMesh object
211: - name - The section name
213: Output Parameters:
214: + numElements - The number of mesh element with values
215: . fiberDim - The number of values per element
216: - array - The array
218: Level: intermediate
220: .keywords: mesh, elements
221: .seealso: DMMeshCreate()
222: @*/
223: PetscErrorCode BCSectionRealGetArray(DM dm, const char name[], PetscInt *numElements, PetscInt *fiberDim, PetscReal *array[])
224: {
225: ALE::Obj<PETSC_MESH_TYPE> m;
226: PetscErrorCode ierr;
229: DMMeshGetMesh(dm, m);
230: const ALE::Obj<PETSC_MESH_TYPE::real_section_type>& section = m->getRealSection(std::string(name));
231: if (!section->size()) {
232: *numElements = 0;
233: *fiberDim = 0;
234: *array = NULL;
235: return(0);
236: }
237: 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: PetscViewer viewer;
292: FILE *f;
293: PetscInt numCells, cellCount = 0;
294: PetscInt *verts;
295: char buf[2048];
296: PetscInt c;
297: PetscInt commRank;
300: MPI_Comm_rank(comm, &commRank);
302: if (commRank != 0) return;
303: PetscViewerCreate(PETSC_COMM_SELF, &viewer);
304: PetscViewerSetType(viewer, PETSCVIEWERASCII);
305: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
306: PetscViewerFileSetName(viewer, filename.c_str());
307: if (ierr) {
308: ostringstream txt;
309: txt << "Could not open PCICE connectivity file: " << filename;
310: throw ALE::Exception(txt.str().c_str());
311: }
312: PetscViewerASCIIGetPointer(viewer, &f);
313: if (fgets(buf, 2048, f) == NULL) {
314: throw ALE::Exception("Invalid connectivity file: Missing number of elements");
315: }
316: const char *sizes = strtok(buf, " ");
317: numCells = atoi(sizes);
318: sizes = strtok(NULL, " ");
319: if (sizes != NULL) {
320: corners = atoi(sizes);
321: std::cout << "Reset corners to " << corners << std::endl;
322: }
323: PetscMalloc(numCells*corners * sizeof(PetscInt), &verts);
324: while(fgets(buf, 2048, f) != NULL) {
325: const char *v = strtok(buf, " ");
327: /* Ignore cell number */
328: v = strtok(NULL, " ");
329: for(c = 0; c < corners; c++) {
330: int vertex = atoi(v);
332: if (!useZeroBase) vertex -= 1;
333: verts[cellCount*corners+c] = vertex;
334: v = strtok(NULL, " ");
335: }
336: cellCount++;
337: }
338: PetscViewerDestroy(&viewer);
339: numElements = numCells;
340: *vertices = verts;
341: };
342: void Builder::readCoordinates(MPI_Comm comm, const std::string& filename, const int dim, int& numVertices, PetscReal *coordinates[]) {
343: PetscViewer viewer;
344: FILE *f;
345: PetscInt numVerts, vertexCount = 0;
346: PetscReal *coords;
347: char buf[2048];
348: PetscInt c;
349: PetscInt commRank;
352: MPI_Comm_rank(comm, &commRank);
354: if (commRank != 0) return;
355: PetscViewerCreate(PETSC_COMM_SELF, &viewer);
356: PetscViewerSetType(viewer, PETSCVIEWERASCII);
357: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
358: PetscViewerFileSetName(viewer, filename.c_str());
359: if (ierr) {
360: ostringstream txt;
361: txt << "Could not open PCICE coordinate file: " << filename;
362: throw ALE::Exception(txt.str().c_str());
363: }
364: PetscViewerASCIIGetPointer(viewer, &f);
365: numVerts = atoi(fgets(buf, 2048, f));
366: PetscMalloc(numVerts*dim * sizeof(PetscReal), &coords);
367: while(fgets(buf, 2048, f) != NULL) {
368: const char *x = strtok(buf, " ");
370: /* Ignore vertex number */
371: x = strtok(NULL, " ");
372: for(c = 0; c < dim; c++) {
373: coords[vertexCount*dim+c] = atof(x);
374: x = strtok(NULL, " ");
375: }
376: vertexCount++;
377: }
378: PetscViewerDestroy(&viewer);
379: numVertices = numVerts;
380: *coordinates = coords;
381: };
382: 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) {
383: return readMesh(comm, dim, basename+".nodes", basename+".lcon", useZeroBase, interpolate, debug);
384: };
385: #ifdef PETSC_OPT_SIEVE
386: 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) {
387: typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
388: Obj<Mesh> mesh = new Mesh(comm, dim, debug);
389: Obj<sieve_type> sieve = new sieve_type(comm, debug);
390: const Obj<FlexMesh> m = new FlexMesh(comm, dim, debug);
391: const Obj<FlexMesh::sieve_type> s = new FlexMesh::sieve_type(comm, debug);
392: int *cells = NULL;
393: PetscReal *coordinates = NULL;
394: int numCells = 0, numVertices = 0, numCorners = dim+1;
397: ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
398: ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
399: ALE::SieveBuilder<FlexMesh>::buildTopology(s, dim, numCells, cells, numVertices, interpolate, numCorners, -1, m->getArrowSection("orientation"));
400: m->setSieve(s);
401: m->stratify();
402: mesh->setSieve(sieve);
403: std::map<Mesh::point_type,Mesh::point_type> renumbering;
404: ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering, false);
405: mesh->stratify();
406: ALE::ISieveConverter::convertOrientation(*s, *sieve, renumbering, m->getArrowSection("orientation").ptr());
407: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
408: if (cells) {PetscFree(cells);CHKERRXX(ierr);}
409: if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
410: return mesh;
411: };
412: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename) {
413: throw ALE::Exception("Not implemented for optimized sieves");
414: };
415: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, PetscReal *coordinates[], const bool columnMajor) {
416: const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
417: if (!coordSec->size()) {
418: *numVertices = 0;
419: *dim = 0;
420: *coordinates = NULL;
421: return;
422: }
423: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
424: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
425: int size = vertices->size();
426: int embedDim = coordSec->getFiberDimension(*vertices->begin());
427: PetscReal *coords;
430: PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
431: for(Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
432: const Mesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);
433: const int row = vNumbering->getIndex(*v_iter);
435: if (columnMajor) {
436: for(int d = 0; d < embedDim; d++) {
437: coords[d*size + row] = array[d];
438: }
439: } else {
440: for(int d = 0; d < embedDim; d++) {
441: coords[row*embedDim + d] = array[d];
442: }
443: }
444: }
445: *numVertices = size;
446: *dim = embedDim;
447: *coordinates = coords;
448: };
449: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor) {
450: if (!mesh->heightStratum(0)->size()) {
451: *numElements = 0;
452: *numCorners = 0;
453: *vertices = NULL;
454: return;
455: }
456: const Obj<Mesh::sieve_type>& sieve = mesh->getSieve();
457: const Obj<Mesh::label_sequence>& elements = mesh->heightStratum(0);
458: const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
459: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
460: int size = elements->size();
461: //int corners = sieve->nCone(*elements->begin(), topology->depth())->size();
462: int corners = sieve->getConeSize(*elements->begin());
463: int *v;
466: PetscMalloc(size*corners * sizeof(int), &v);CHKERRXX(ierr);
467: for(Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
468: const Obj<Mesh::sieve_type::coneSequence> cone = sieve->cone(*e_iter);
469: Mesh::sieve_type::coneSequence::const_iterator begin = cone->begin();
470: Mesh::sieve_type::coneSequence::const_iterator end = cone->end();
472: const int row = eNumbering->getIndex(*e_iter);
473: int c = -1;
474: if (columnMajor) {
475: for(Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
476: v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
477: }
478: } else {
479: for(Mesh::sieve_type::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
480: v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
481: }
482: }
483: }
484: *numElements = size;
485: *numCorners = corners;
486: *vertices = v;
487: };
488: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer) {
489: throw ALE::Exception("Not implemented for optimized sieves");
490: };
491: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer) {
492: throw ALE::Exception("Not implemented for optimized sieves");
493: };
494: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer) {
495: throw ALE::Exception("Not implemented for optimized sieves");
496: };
497: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer) {
498: throw ALE::Exception("Not implemented for optimized sieves");
499: };
500: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh) {
501: throw ALE::Exception("Not implemented for optimized sieves");
502: };
503: #else
504: 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) {
505: Obj<Mesh> mesh = new DMMesh(comm, dim, debug);
506: Obj<sieve_type> sieve = new sieve_type(comm, debug);
507: int *cells = NULL;
508: double *coordinates = NULL;
509: int numCells = 0, numVertices = 0, numCorners = dim+1;
512: ALE::PCICE::Builder::readConnectivity(comm, adjFilename, numCorners, useZeroBase, numCells, &cells);
513: ALE::PCICE::Builder::readCoordinates(comm, coordFilename, dim, numVertices, &coordinates);
514: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildTopology(sieve, dim, numCells, cells, numVertices, interpolate, numCorners, -1, mesh->getArrowSection("orientation"));
515: mesh->setSieve(sieve);
516: mesh->stratify();
517: ALE::SieveBuilder<PETSC_MESH_TYPE>::buildCoordinates(mesh, dim, coordinates);
518: if (cells) {PetscFree(cells);CHKERRXX(ierr);}
519: if (coordinates) {PetscFree(coordinates);CHKERRXX(ierr);}
520: return mesh;
521: };
522: // Creates boundary sections:
523: // IBC[NBFS,2]: ALL
524: // BL[NBFS,1]:
525: // BNVEC[NBFS,2]:
526: // BCFUNC[NBCF,NV]: ALL
527: // IBNDFS[NBN,2]: STILL NEED 4-5
528: // BNNV[NBN,2]
529: void Builder::readBoundary(const Obj<Mesh>& mesh, const std::string& bcFilename) {
530: PetscViewer viewer;
531: FILE *f;
532: char buf[2048];
535: const Obj<Mesh::int_section_type>& ibc = mesh->getIntSection("IBC");
536: const Obj<Mesh::int_section_type>& ibndfs = mesh->getIntSection("IBNDFS");
537: const Obj<Mesh::int_section_type>& ibcnum = mesh->getIntSection("IBCNUM");
538: const Obj<Mesh::int_section_type>& ibfcon = mesh->getIntSection("IBFCON");
539: const Obj<Mesh::real_section_type>& bl = mesh->getRealSection("BL");
540: const Obj<Mesh::real_section_type>& bnvec = mesh->getRealSection("BNVEC");
541: const Obj<Mesh::real_section_type>& bnnv = mesh->getRealSection("BNNV");
542: if (mesh->commRank() != 0) {
543: #if 0
544: mesh->distributeBCValues();
545: #endif
546: return;
547: }
548: PetscViewerCreate(PETSC_COMM_SELF, &viewer);CHKERRXX(ierr);
549: PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRXX(ierr);
550: PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRXX(ierr);
551: PetscViewerFileSetName(viewer, bcFilename.c_str());CHKERRXX(ierr);
552: PetscViewerASCIIGetPointer(viewer, &f);CHKERRXX(ierr);
553: // Create IBC section
554: int numBdFaces = atoi(strtok(fgets(buf, 2048, f), " "));
555: int *tmpIBC = new int[numBdFaces*4];
556: std::map<int,std::set<int> > elem2Idx;
557: std::map<int,int> bfReorder;
558: for(int bf = 0; bf < numBdFaces; bf++) {
559: const char *x = strtok(fgets(buf, 2048, f), " ");
561: // Ignore boundary face number
562: x = strtok(NULL, " ");
563: tmpIBC[bf*4+0] = atoi(x);
564: x = strtok(NULL, " ");
565: tmpIBC[bf*4+1] = atoi(x);
566: x = strtok(NULL, " ");
567: tmpIBC[bf*4+2] = atoi(x);
568: x = strtok(NULL, " ");
569: tmpIBC[bf*4+3] = atoi(x);
570: const int elem = tmpIBC[bf*4+0]-1;
572: ibc->addFiberDimension(elem, 4);
573: ibcnum->addFiberDimension(elem, 1);
574: ibfcon->addFiberDimension(elem, 2);
575: bl->addFiberDimension(elem, 1);
576: bnvec->addFiberDimension(elem, 2);
577: elem2Idx[elem].insert(bf);
578: }
579: mesh->allocate(ibc);
580: mesh->allocate(ibcnum);
581: mesh->allocate(ibfcon);
582: mesh->allocate(bl);
583: mesh->allocate(bnvec);
584: const DMMesh::int_section_type::chart_type& chart = ibc->getChart();
585: int num = 1;
587: for(Mesh::int_section_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
588: const int elem = *p_iter;
589: int bfNum[2];
590: int k = 0;
592: for(std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
593: bfReorder[(*i_iter)+1] = num;
594: bfNum[k++] = num;
595: num++;
596: }
597: ibcnum->updatePoint(elem, bfNum);
598: }
599: for(int bf = 0; bf < numBdFaces; bf++) {
600: const int elem = tmpIBC[bf*4]-1;
602: if (elem2Idx[elem].size() > 1) {
603: if (*elem2Idx[elem].begin() == bf) {
604: int values[8];
605: int k = 0;
607: for(std::set<int>::const_iterator i_iter = elem2Idx[elem].begin(); i_iter != elem2Idx[elem].end(); ++i_iter) {
608: for(int v = 0; v < 4; ++v) {
609: values[k*4+v] = tmpIBC[*i_iter*4+v];
610: }
611: k++;
612: }
613: ibc->updatePoint(elem, values);
614: }
615: } else {
616: ibc->updatePoint(elem, &tmpIBC[bf*4]);
617: }
618: }
619: delete [] tmpIBC;
620: // Create BCFUNC section
621: int numBcFunc = atoi(strtok(fgets(buf, 2048, f), " "));
622: if (numBcFunc != 0) {throw ALE::Exception("Cannot handle BCFUNCS after rewrite");}
623: for(int bc = 0; bc < numBcFunc; bc++) {
624: #if 0
625: const char *x = strtok(fgets(buf, 2048, f), " ");
626: DMMesh::bc_value_type value;
628: // Ignore function number
629: x = strtok(NULL, " ");
630: value.rho = atof(x);
631: x = strtok(NULL, " ");
632: value.u = atof(x);
633: x = strtok(NULL, " ");
634: value.v = atof(x);
635: x = strtok(NULL, " ");
636: value.p = atof(x);
637: mesh->setBCValue(bc+1, value);
638: #endif
639: }
640: #if 0
641: mesh->distributeBCValues();
642: #endif
643: // Create IBNDFS section
644: int numBdVertices = atoi(strtok(fgets(buf, 2048, f), " "));
645: const int numElements = mesh->heightStratum(0)->size();
646: int *tmpIBNDFS = new int[numBdVertices*3];
648: for(int bv = 0; bv < numBdVertices; bv++) {
649: const char *x = strtok(fgets(buf, 2048, f), " ");
651: // Ignore boundary node number
652: x = strtok(NULL, " ");
653: tmpIBNDFS[bv*3+0] = atoi(x);
654: x = strtok(NULL, " ");
655: tmpIBNDFS[bv*3+1] = atoi(x);
656: x = strtok(NULL, " ");
657: tmpIBNDFS[bv*3+2] = atoi(x);
658: ibndfs->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, 6);
659: }
660: mesh->allocate(ibndfs);
661: for(int bv = 0; bv < numBdVertices; bv++) {
662: int values[5];
664: values[0] = tmpIBNDFS[bv*3+0];
665: // Covert to new boundary face numbers
666: values[1] = bfReorder[tmpIBNDFS[bv*3+1]];
667: values[2] = bfReorder[tmpIBNDFS[bv*3+2]];
668: values[3] = 0;
669: values[4] = 0;
670: ibndfs->updatePoint(values[0]-1+numElements, values);
671: }
672: PetscViewerDestroy(&viewer);CHKERRXX(ierr);
673: // Create BNNV[NBN,2]
674: const int dim = mesh->getDimension();
676: for(int bv = 0; bv < numBdVertices; bv++) {
677: bnnv->setFiberDimension(tmpIBNDFS[bv*3+0]-1+numElements, dim);
678: }
679: mesh->allocate(bnnv);
680: delete [] tmpIBNDFS;
681: };
682: void Builder::outputVerticesLocal(const Obj<Mesh>& mesh, int *numVertices, int *dim, double *coordinates[], const bool columnMajor) {
683: const Obj<Mesh::real_section_type>& coordSec = mesh->getRealSection("coordinates");
684: if (!coordSec->size()) {
685: *numVertices = 0;
686: *dim = 0;
687: *coordinates = NULL;
688: return;
689: }
690: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
691: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
692: int size = vertices->size();
693: int embedDim = coordSec->getFiberDimension(*vertices->begin());
694: double *coords;
697: PetscMalloc(vertices->size()*embedDim * sizeof(double), &coords);CHKERRXX(ierr);
698: for(Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
699: const DMMesh::real_section_type::value_type *array = coordSec->restrictPoint(*v_iter);
700: const int row = vNumbering->getIndex(*v_iter);
702: if (columnMajor) {
703: for(int d = 0; d < embedDim; d++) {
704: coords[d*size + row] = array[d];
705: }
706: } else {
707: for(int d = 0; d < embedDim; d++) {
708: coords[row*embedDim + d] = array[d];
709: }
710: }
711: }
712: *numVertices = size;
713: *dim = embedDim;
714: *coordinates = coords;
715: };
716: void Builder::outputElementsLocal(const Obj<Mesh>& mesh, int *numElements, int *numCorners, int *vertices[], const bool columnMajor) {
717: if (!mesh->heightStratum(0)->size()) {
718: *numElements = 0;
719: *numCorners = 0;
720: *vertices = NULL;
721: return;
722: }
723: const Obj<Mesh::sieve_type>& sieve = mesh->getSieve();
724: const Obj<Mesh::label_sequence>& elements = mesh->heightStratum(0);
725: const Obj<Mesh::numbering_type>& eNumbering = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
726: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
727: int size = elements->size();
728: //int corners = sieve->nCone(*elements->begin(), topology->depth())->size();
729: int corners = sieve->getConeSize(*elements->begin());
730: int *v;
733: PetscMalloc(elements->size()*corners * sizeof(int), &v);CHKERRXX(ierr);
734: for(Mesh::label_sequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
735: const Obj<Mesh::sieve_type::traits::coneSequence> cone = sieve->cone(*e_iter);
736: DMMesh::sieve_type::traits::coneSequence::iterator begin = cone->begin();
737: DMMesh::sieve_type::traits::coneSequence::iterator end = cone->end();
739: const int row = eNumbering->getIndex(*e_iter);
740: int c = -1;
741: if (columnMajor) {
742: for(Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
743: v[(++c)*size + row] = vNumbering->getIndex(*c_iter)+1;
744: }
745: } else {
746: for(Mesh::sieve_type::traits::coneSequence::iterator c_iter = begin; c_iter != end; ++c_iter) {
747: v[row*corners + ++c] = vNumbering->getIndex(*c_iter)+1;
748: }
749: }
750: }
751: *numElements = size;
752: *numCorners = corners;
753: *vertices = v;
754: };
757: PetscErrorCode Viewer::writeVertices(const ALE::Obj<Mesh>& mesh, PetscViewer viewer) {
758: ALE::Obj<Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
759: #if 0
760: DMMesh::field_type::patch_type patch;
761: const double *array = coordinates->restrict(patch);
762: int numVertices;
766: //FIX:
767: if (vertexBundle->getGlobalOffsets()) {
768: numVertices = vertexBundle->getGlobalOffsets()[mesh->commSize()];
769: } else {
770: numVertices = mesh->getTopology()->depthStratum(0)->size();
771: }
772: PetscViewerASCIIPrintf(viewer, "%D\n", numVertices);
773: if (mesh->commRank() == 0) {
774: int numLocalVertices = mesh->getTopology()->depthStratum(0)->size();
775: int embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
776: int vertexCount = 1;
778: for(int v = 0; v < numLocalVertices; v++) {
779: PetscViewerASCIIPrintf(viewer, "%7D ", vertexCount++);
780: for(int d = 0; d < embedDim; d++) {
781: if (d > 0) {
782: PetscViewerASCIIPrintf(viewer, " ");
783: }
784: PetscViewerASCIIPrintf(viewer, "% 12.5E", array[v*embedDim+d]);
785: }
786: PetscViewerASCIIPrintf(viewer, "\n");
787: }
788: for(int p = 1; p < mesh->commSize(); p++) {
789: double *remoteCoords;
790: MPI_Status status;
792: MPI_Recv(&numLocalVertices, 1, MPI_INT, p, 1, mesh->comm(), &status);
793: PetscMalloc(numLocalVertices*embedDim * sizeof(double), &remoteCoords);
794: MPI_Recv(remoteCoords, numLocalVertices*embedDim, MPI_DOUBLE, p, 1, mesh->comm(), &status);
795: for(int v = 0; v < numLocalVertices; v++) {
796: PetscViewerASCIIPrintf(viewer,"%7D ", vertexCount++);
797: for(int d = 0; d < embedDim; d++) {
798: if (d > 0) {
799: PetscViewerASCIIPrintf(viewer, " ");
800: }
801: PetscViewerASCIIPrintf(viewer, "% 12.5E", remoteCoords[v*embedDim+d]);
802: }
803: PetscViewerASCIIPrintf(viewer, "\n");
804: }
805: }
806: } else {
807: ALE::Obj<Mesh::bundle_type> globalOrder = coordinates->getGlobalOrder();
808: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = globalOrder->getPatch(patch);
809: const int *offsets = coordinates->getGlobalOffsets();
810: int embedDim = coordinates->getFiberDimension(patch, *mesh->getTopology()->depthStratum(0)->begin());
811: int numLocalVertices = (offsets[mesh->commRank()+1] - offsets[mesh->commRank()])/embedDim;
812: double *localCoords;
813: int k = 0;
815: PetscMalloc(numLocalVertices*embedDim * sizeof(double), &localCoords);
816: for(Mesh::bundle_type::order_type::coneSequence::iterator p_iter = cone->begin(); p_iter != cone->end(); ++p_iter) {
817: int dim = globalOrder->getFiberDimension(patch, *p_iter);
819: if (dim > 0) {
820: int offset = coordinates->getFiberOffset(patch, *p_iter);
822: for(int i = offset; i < offset+dim; ++i) {
823: localCoords[k++] = array[i];
824: }
825: }
826: }
827: if (k != numLocalVertices*embedDim) {
828: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of coordinates to send %d should be %d", k, numLocalVertices*embedDim);
829: }
830: MPI_Send(&numLocalVertices, 1, MPI_INT, 0, 1, mesh->comm());
831: MPI_Send(localCoords, numLocalVertices*embedDim, MPI_DOUBLE, 0, 1, mesh->comm());
832: PetscFree(localCoords);
833: }
834: #endif
835: return(0);
836: };
839: PetscErrorCode Viewer::writeElements(const ALE::Obj<Mesh>& mesh, PetscViewer viewer) {
840: #if 0
841: ALE::Obj<Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
842: ALE::Obj<Mesh::bundle_type> elementBundle = mesh->getBundle(topology->depth());
843: ALE::Obj<Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
844: ALE::Obj<Mesh::bundle_type> globalVertex = vertexBundle->getGlobalOrder();
845: ALE::Obj<Mesh::bundle_type> globalElement = elementBundle->getGlobalOrder();
846: DMMesh::bundle_type::patch_type patch;
847: std::string orderName("element");
848: int dim = mesh->getDimension();
849: int corners = topology->nCone(*elements->begin(), topology->depth())->size();
850: int numElements;
854: if (corners != dim+1) {
855: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "PCICE only supports simplicies");
856: }
857: if (!globalVertex) {
858: globalVertex = vertexBundle;
859: }
860: if (elementBundle->getGlobalOffsets()) {
861: numElements = elementBundle->getGlobalOffsets()[mesh->commSize()];
862: } else {
863: numElements = mesh->getTopology()->heightStratum(0)->size();
864: }
865: if (mesh->commRank() == 0) {
866: int elementCount = 1;
868: PetscViewerASCIIPrintf(viewer, "%d\n", numElements);
869: for(Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
870: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);
872: PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
873: for(Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
874: PetscViewerASCIIPrintf(viewer, " %7d", globalVertex->getIndex(patch, *c_itor).prefix);
875: }
876: PetscViewerASCIIPrintf(viewer, "\n");
877: }
878: for(int p = 1; p < mesh->commSize(); p++) {
879: int numLocalElements;
880: int *remoteVertices;
881: MPI_Status status;
883: MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
884: PetscMalloc(numLocalElements*corners * sizeof(int), &remoteVertices);
885: MPI_Recv(remoteVertices, numLocalElements*corners, MPI_INT, p, 1, mesh->comm(), &status);
886: for(int e = 0; e < numLocalElements; e++) {
887: PetscViewerASCIIPrintf(viewer, "%7d", elementCount++);
888: for(int c = 0; c < corners; c++) {
889: PetscViewerASCIIPrintf(viewer, " %7d", remoteVertices[e*corners+c]);
890: }
891: PetscViewerASCIIPrintf(viewer, "\n");
892: }
893: PetscFree(remoteVertices);
894: }
895: } else {
896: const int *offsets = elementBundle->getGlobalOffsets();
897: int numLocalElements = offsets[mesh->commRank()+1] - offsets[mesh->commRank()];
898: int *localVertices;
899: int k = 0;
901: PetscMalloc(numLocalElements*corners * sizeof(int), &localVertices);
902: for(Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
903: ALE::Obj<Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_itor);
905: if (globalElement->getFiberDimension(patch, *e_itor) > 0) {
906: for(Mesh::bundle_type::order_type::coneSequence::iterator c_itor = cone->begin(); c_itor != cone->end(); ++c_itor) {
907: localVertices[k++] = globalVertex->getIndex(patch, *c_itor).prefix;
908: }
909: }
910: }
911: if (k != numLocalElements*corners) {
912: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of vertices to send %d should be %d", k, numLocalElements*corners);
913: }
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: };
923: PetscErrorCode Viewer::writeVerticesLocal(const Obj<Mesh>& mesh, PetscViewer viewer) {
924: Obj<Mesh::real_section_type> coordinates = mesh->getRealSection("coordinates");
925: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
926: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getLocalNumbering(mesh, 0);
927: int embedDim = coordinates->getFiberDimension(*vertices->begin());
931: PetscViewerASCIIPrintf(viewer, "%D\n", vertices->size());
932: for(Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
933: const DMMesh::real_section_type::value_type *array = coordinates->restrictPoint(*v_iter);
935: PetscViewerASCIIPrintf(viewer, "%7D ", vNumbering->getIndex(*v_iter)+1);
936: for(int d = 0; d < embedDim; d++) {
937: if (d > 0) {
938: PetscViewerASCIIPrintf(viewer, " ");
939: }
940: PetscViewerASCIIPrintf(viewer, "% 12.5E", array[d]);
941: }
942: PetscViewerASCIIPrintf(viewer, "\n");
943: }
944: return(0);
945: };
948: PetscErrorCode Viewer::writeRestart(const Obj<Mesh>& mesh, PetscViewer viewer) {
949: const Obj<Mesh::real_section_type>& velocity = mesh->getRealSection("VELN");
950: const Obj<Mesh::real_section_type>& pressure = mesh->getRealSection("PN");
951: const Obj<Mesh::real_section_type>& temperature = mesh->getRealSection("TN");
952: const Obj<Mesh::numbering_type>& cNumbering = mesh->getFactory()->getNumbering(mesh, mesh->depth());
953: const Obj<Mesh::numbering_type>& vNumbering = mesh->getFactory()->getNumbering(mesh, 0);
954: const int numCells = cNumbering->getGlobalSize();
958: int blen[2];
959: MPI_Aint indices[2];
960: MPI_Datatype oldtypes[2], newtype;
961: blen[0] = 1; indices[0] = 0; oldtypes[0] = MPI_INT;
962: blen[1] = 4; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
963: MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
964: MPI_Type_commit(&newtype);
966: if (mesh->commRank() == 0) {
967: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
969: for(Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
970: if (vNumbering->isLocal(*v_iter)) {
971: const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
972: const DMMesh::real_section_type::value_type *pn = pressure->restrictPoint(*v_iter);
973: const DMMesh::real_section_type::value_type *tn = temperature->restrictPoint(*v_iter);
975: PetscViewerASCIIPrintf(viewer, "%6d% 16.8E% 16.8E% 16.8E% 16.8E\n", *v_iter-numCells+1, veln[0], veln[1], pn[0], tn[0]);
976: }
977: }
978: for(int p = 1; p < mesh->commSize(); p++) {
979: RestartType *remoteValues;
980: int numLocalElements;
981: MPI_Status status;
983: MPI_Recv(&numLocalElements, 1, MPI_INT, p, 1, mesh->comm(), &status);
984: PetscMalloc(numLocalElements * sizeof(RestartType), &remoteValues);
985: MPI_Recv(remoteValues, numLocalElements, newtype, p, 1, mesh->comm(), &status);
986: for(int e = 0; e < numLocalElements; e++) {
987: 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);
988: }
989: }
990: } else {
991: const Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
992: RestartType *localValues;
993: int numLocalElements = vNumbering->getLocalSize();
994: int k = 0;
996: PetscMalloc(numLocalElements * sizeof(RestartType), &localValues);
997: for(Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
998: if (vNumbering->isLocal(*v_iter)) {
999: const DMMesh::real_section_type::value_type *veln = velocity->restrictPoint(*v_iter);
1000: const DMMesh::real_section_type::value_type *pn = pressure->restrictPoint(*v_iter);
1001: const DMMesh::real_section_type::value_type *tn = temperature->restrictPoint(*v_iter);
1003: localValues[k].vertex = *v_iter;
1004: localValues[k].veln_x = veln[0];
1005: localValues[k].veln_y = veln[1];
1006: localValues[k].pn = pn[0];
1007: localValues[k].tn = tn[0];
1008: k++;
1009: }
1010: }
1011: if (k != numLocalElements) {
1012: SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "Invalid number of values to send for field, %d should be %d", k, numLocalElements);
1013: }
1014: MPI_Send(&numLocalElements, 1, MPI_INT, 0, 1, mesh->comm());
1015: MPI_Send(localValues, numLocalElements, newtype, 0, 1, mesh->comm());
1016: PetscFree(localValues);
1017: }
1018: MPI_Type_free(&newtype);
1019: return(0);
1020: };
1022: // This class reconstructs the local pieces of the boundary that distributed PCICE needs.
1023: // The boundary along with the boundary conditions is encoded in a collection of sections
1024: // over the PCICE mesh. These sections contain a description of the boundary topology
1025: // using elements' global names. This is unacceptable for PCICE, since it interprets
1026: // elements of the connectivity data arrays as local offsets into (some other of) these arrays.
1027: // This subroutine performs the renumbering based on the local numbering on the distributed
1028: // mesh. Essentially, we replace each global element id by its local number.
1029: //
1030: // Note: Any vertex or element number from PCICE is 1-based, but in Sieve we are 0-based. Thus
1031: // we add and subtract 1 during conversion. Also, Sieve vertices are numbered after cells.
1032: void fuseBoundary(const ALE::Obj<PETSC_MESH_TYPE>& mesh) {
1033: // Extract PCICE boundary sections
1034: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCsec = mesh->getIntSection("IBC");
1035: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBNDFSsec = mesh->getIntSection("IBNDFS");
1036: ALE::Obj<PETSC_MESH_TYPE::int_section_type> IBCNUMsec = mesh->getIntSection("IBCNUM");
1038: // Look at the sections, if debugging
1039: if (mesh->debug()) {
1040: IBCsec->view("IBCsec", mesh->comm());
1041: IBNDFSsec->view("IBNDFSsec", mesh->comm());
1042: }
1044: // Extract the local numberings
1045: ALE::Obj<PETSC_MESH_TYPE::numbering_type> vertexNum = mesh->getFactory()->getLocalNumbering(mesh, 0);
1046: ALE::Obj<PETSC_MESH_TYPE::numbering_type> elementNum = mesh->getFactory()->getLocalNumbering(mesh, mesh->depth());
1047: const int numElements = mesh->getFactory()->getNumbering(mesh, mesh->depth())->getGlobalSize();
1048: std::map<int,int> bfMap;
1049: // Declare points to the extracted numbering data
1050: const PETSC_MESH_TYPE::numbering_type::value_type *vNum, *eNum;
1052: // Create map from serial bdFace numbers to local bdFace numbers
1053: {
1054: const PETSC_MESH_TYPE::int_section_type::chart_type chart = IBCNUMsec->getChart();
1055: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator begin = chart.begin();
1056: PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator end = chart.end();
1057: int num = 0;
1059: for(PETSC_MESH_TYPE::int_section_type::chart_type::const_iterator p_iter = begin; p_iter != end; ++p_iter) {
1060: const int fiberDim = IBCNUMsec->getFiberDimension(*p_iter);
1061: const int *globalNum = IBCNUMsec->restrictPoint(*p_iter);
1063: for(int n = 0; n < fiberDim; ++n) {
1064: bfMap[globalNum[n]] = ++num;
1065: }
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: };