00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceExodusWriter.hpp"
00032 #include "PlayaExceptions.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "SundanceMaximalCellFilter.hpp"
00035 #include "SundanceCellFilter.hpp"
00036 #include "SundanceVertexSort.hpp"
00037 #include "PlayaTabs.hpp"
00038 #include "Teuchos_XMLObject.hpp"
00039
00040 #ifdef HAVE_SUNDANCE_EXODUS
00041 #include "exodusII.h"
00042 #endif
00043
00044
00045 using namespace Sundance;
00046 using namespace Teuchos;
00047 using namespace std;
00048
00049
00050 void ExodusWriter::write() const
00051 {
00052
00053 std::string exoFile = filename();
00054 std::string parFile = filename();
00055 if (nProc() > 1)
00056 {
00057 exoFile = exoFile + "-" + Teuchos::toString(nProc()) + "-" + Teuchos::toString(myRank());
00058 parFile = parFile + "-" + Teuchos::toString(nProc()) + "-" + Teuchos::toString(myRank());
00059 }
00060 exoFile = exoFile + ".exo";
00061 parFile = parFile + ".pxo";
00062
00063 if (nProc() > 1) writeParallelInfo(parFile);
00064 #ifdef HAVE_SUNDANCE_EXODUS
00065 int ws = 8;
00066 int exoid = ex_create(exoFile.c_str(), EX_CLOBBER, &ws, &ws);
00067
00068 TEUCHOS_TEST_FOR_EXCEPTION(exoid < 0, std::runtime_error, "failure to create file "
00069 << filename());
00070
00071
00072 Array<CellFilter> nsFilters;
00073 Array<int> omniNodalFuncs;
00074 Array<RCP<Array<int> > > funcsForNodeset;
00075 Array<RCP<Array<int> > > nodesForNodeset;
00076 Array<int> nsID;
00077 Array<int> nsNodesPerSet;
00078 Array<int> nsNodePtr;
00079 RCP<Array<int> > allNodes=rcp(new Array<int>());
00080
00081 Array<CellFilter> blockFilters;
00082 Array<int> omniElemFuncs;
00083 Array<RCP<Array<int> > > funcsForBlock;
00084 Array<RCP<Array<int> > > elemsForBlock;
00085 Array<int> blockID;
00086 Array<int> nElemsPerBlock;
00087 Array<int> blockElemPtr;
00088 RCP<Array<int> > allElems=rcp(new Array<int>());
00089
00090
00091 findNodeSets(nsFilters, omniNodalFuncs, funcsForNodeset,
00092 nodesForNodeset, nsID, nsNodesPerSet, nsNodePtr, allNodes);
00093
00094 findBlocks(blockFilters, omniElemFuncs, funcsForBlock,
00095 elemsForBlock, blockID, nElemsPerBlock, blockElemPtr, allElems);
00096
00097
00098
00099 writeMesh(exoid, nsFilters, nsID, nsNodesPerSet, nsNodePtr, allNodes );
00100
00101 writeFields(exoid, nsFilters, omniNodalFuncs, omniElemFuncs,
00102 funcsForNodeset,
00103 nodesForNodeset, nsID);
00104
00105
00106 ex_close(exoid);
00107 #else
00108 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Exodus not enabled");
00109 #endif
00110 }
00111
00112
00113 void ExodusWriter::offset(Array<int>& x) const
00114 {
00115 for (int i=0; i<x.size(); i++) x[i]++;
00116 }
00117
00118
00119 void ExodusWriter::writeMesh(int exoid,
00120 const Array<CellFilter>& nodesetFilters,
00121 const Array<int>& nsID,
00122 const Array<int>& nNodesPerSet,
00123 const Array<int>& nsNodePtr,
00124 const RCP<Array<int> >& allNodes) const
00125 {
00126 #ifdef HAVE_SUNDANCE_EXODUS
00127
00128 int ierr = 0;
00129
00130 int dim = mesh().spatialDim();
00131 int nElems = mesh().numCells(dim);
00132
00133 Array<int> ssLabels = mesh().getAllLabelsForDimension(dim-1).elements();
00134 int numSS = 0;
00135
00136 for (int ss=0; ss<ssLabels.size(); ss++)
00137 {
00138 if (ssLabels[ss] != 0) numSS++;
00139 }
00140
00141 int numNS = nsID.size();
00142
00143
00144
00145 int nElemBlocks = mesh().numLabels(dim);
00146
00147 ierr = ex_put_init(
00148 exoid,
00149 filename().c_str(),
00150 dim,
00151 mesh().numCells(0),
00152 nElems,
00153 nElemBlocks,
00154 numNS,
00155 numSS
00156 );
00157 TEUCHOS_TEST_FOR_EXCEPT(ierr<0);
00158
00159
00160
00161 int nPts = mesh().numCells(0);
00162
00163 Array<double> x(nPts);
00164 Array<double> y(nPts);
00165 Array<double> z(nPts * (dim > 2));
00166
00167 for (int n=0; n<nPts; n++)
00168 {
00169 const Point& P = mesh().nodePosition(n);
00170 x[n] = P[0];
00171 y[n] = P[1];
00172 if (dim==3) z[n] = P[2];
00173 }
00174
00175 if (dim==2)
00176 {
00177 ierr = ex_put_coord(exoid, &(x[0]), &(y[0]), (void*) 0);
00178 }
00179 else
00180 {
00181 ierr = ex_put_coord(exoid, &(x[0]), &(y[0]), &(z[0]));
00182 }
00183
00184 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00185
00186 if (dim==2)
00187 {
00188 Array<std::string> cn;
00189 cn.append("x");
00190 cn.append("y");
00191 Array<const char*> pp;
00192 getCharpp(cn, pp);
00193 ierr = ex_put_coord_names(exoid,(char**) &(pp[0]));
00194 }
00195 else
00196 {
00197 Array<std::string> cn;
00198 cn.append("x");
00199 cn.append("y");
00200 cn.append("z");
00201 Array<const char*> pp;
00202 getCharpp(cn, pp);
00203 ierr = ex_put_coord_names(exoid, (char**)&(pp[0]));
00204 }
00205
00206 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00207
00208
00209
00210 Array<int> blockLabels = mesh().getAllLabelsForDimension(dim).elements();
00211 int nodesPerElem = dim+1;
00212 std::string eType = elemType(mesh().cellType(dim));
00213
00214 bool minBlockIsZero = false;
00215 for (int b=0; b<blockLabels.size(); b++)
00216 {
00217 if (blockLabels[b] < 1) {minBlockIsZero = true; break;}
00218 }
00219
00220
00221 for (int b=0; b<blockLabels.size(); b++)
00222 {
00223
00224 int numBlockAttr = 0;
00225 Array<int> blockElemLIDs;
00226 Array<int> nodeLIDs;
00227 Array<int> orient;
00228 mesh().getLIDsForLabel(dim, blockLabels[b], blockElemLIDs);
00229 int numElemsThisBlock = blockElemLIDs.size();
00230
00231
00232 mesh().getFacetLIDs(dim, blockElemLIDs, 0, nodeLIDs, orient);
00233 offset(nodeLIDs);
00234 ierr = ex_put_elem_block(
00235 exoid, blockLabels[b]+minBlockIsZero, eType.c_str(),
00236 numElemsThisBlock, nodesPerElem, numBlockAttr
00237 );
00238
00239 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00240
00241 ierr = ex_put_elem_conn(exoid, blockLabels[b]+minBlockIsZero,
00242 &(nodeLIDs[0]));
00243 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00244 }
00245
00246
00247 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00248
00249
00250
00251
00252 for (int ss=0; ss<ssLabels.size(); ss++)
00253 {
00254
00255 if (ssLabels[ss]==0) continue;
00256 Array<int> sideLIDs;
00257 RCP<Array<int> > elemLIDs = rcp(new Array<int>());
00258 RCP<Array<int> > facets = rcp(new Array<int>());
00259 MaximalCofacetBatch maxCofacetBatch;
00260
00261
00262 mesh().getLIDsForLabel(dim-1, ssLabels[ss], sideLIDs);
00263 sort(&(sideLIDs[0]), &(sideLIDs[sideLIDs.size()-1]));
00264
00265
00266 mesh().getMaxCofacetLIDs(sideLIDs, maxCofacetBatch);
00267
00268 maxCofacetBatch.getSpecifiedCofacets(0, elemLIDs, facets);
00269
00270 int numSides = sideLIDs.size();
00271 int numDists = 0;
00272
00273 for (int i=0; i<elemLIDs->size(); i++)
00274 {
00275 (*facets)[i] = ufcFacetIndexToExFacetIndex(dim, (*facets)[i]);
00276 }
00277
00278 offset(sideLIDs);
00279 offset(*elemLIDs);
00280 offset(*facets);
00281
00282
00283 ierr = ex_put_side_set_param(exoid, ssLabels[ss], numSides, numDists);
00284 ierr = ex_put_side_set(exoid, ssLabels[ss], &((*elemLIDs)[0]), &((*facets)[0]));
00285 }
00286
00287
00288 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00289
00290
00291
00292 if (nsID.size() > 0)
00293 {
00294
00295
00296 Array<int> nsDistPerSet(nsID.size(), 0);
00297 Array<int> nsDistPtr(nsID.size(), 0);
00298 Array<int> emptyDist(1, 0);
00299
00300 offset(*allNodes);
00301
00302 ierr = ex_put_concat_node_sets( exoid,
00303 (int*) &(nsID[0]),
00304 (int*) &(nNodesPerSet[0]),
00305 &(nsDistPerSet[0]),
00306 (int*) &(nsNodePtr[0]),
00307 &(nsDistPtr[0]),
00308 &((*allNodes)[0]),
00309 &(emptyDist[0]));
00310
00311 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00312 }
00313
00314 #else
00315 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Exodus not enabled");
00316 #endif
00317 }
00318
00319
00320 void ExodusWriter::writeFields(int exoid,
00321 const Array<CellFilter>& nodesetFilters,
00322 const Array<int>& omnipresentNodalFuncs,
00323 const Array<int>& omnipresentElemFuncs,
00324 const Array<RCP<Array<int> > >& funcsForNodeset,
00325 const Array<RCP<Array<int> > >& nodesForNodeset,
00326 const Array<int>& nsID) const
00327 {
00328
00329 #ifdef HAVE_SUNDANCE_EXODUS
00330 int nNodalFuncs = omnipresentNodalFuncs().size();
00331 int nElemFuncs = omnipresentElemFuncs().size();
00332
00333 int nNodesetFuncs = pointScalarFields().size() - nNodalFuncs;
00334
00335 int nNodesets = funcsForNodeset.size();
00336
00337
00338
00339 Set<int> nsFuncSet;
00340 Map<int, Array<int> > funcToNSMap;
00341
00342 for (int i=0; i<nNodesets; i++)
00343 {
00344 const Array<int>& f = *(funcsForNodeset[i]);
00345 for (int j=0; j<f.size(); j++)
00346 {
00347 nsFuncSet.put(f[j]);
00348 if (funcToNSMap.containsKey(f[j]))
00349 {
00350 funcToNSMap[f[j]].append(i);
00351 }
00352 else
00353 {
00354 funcToNSMap.put(f[j], tuple(i));
00355 }
00356 }
00357 }
00358 Array<int> nsFuncs = nsFuncSet.elements();
00359 TEUCHOS_TEST_FOR_EXCEPT(nsFuncs.size() != nNodesetFuncs);
00360
00361 Map<int, int > funcIDToNSFuncIndex;
00362 for (int i=0; i<nNodesetFuncs; i++) funcIDToNSFuncIndex.put(nsFuncs[i],i);
00363
00364 Array<Array<int> > nsFuncNodesets(nsFuncs.size());
00365 for (int i=0; i<nNodesetFuncs; i++)
00366 {
00367 nsFuncNodesets[i] = funcToNSMap.get(nsFuncs[i]);
00368 }
00369
00370 Array<int> nodesetFuncTruthTable(nNodesetFuncs * nNodesets, 0);
00371 for (int i=0; i<nNodesetFuncs; i++)
00372 {
00373 for (int j=0; j<nsFuncNodesets[i].size(); j++)
00374 {
00375 int ns = nsFuncNodesets[i][j];
00376 nodesetFuncTruthTable[ns*nNodesetFuncs + i] = 1;
00377 }
00378
00379 nsFuncNodesets[i] = funcToNSMap.get(nsFuncs[i]);
00380 }
00381
00382
00383
00384 int ierr;
00385
00386 if (nNodalFuncs > 0)
00387 {
00388 ierr = ex_put_var_param(exoid, "N", nNodalFuncs);
00389 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00390 }
00391
00392 if (nElemFuncs > 0)
00393 {
00394 ierr = ex_put_var_param(exoid, "E", nElemFuncs);
00395 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00396 }
00397
00398
00399 if (nNodesets > 0)
00400 {
00401 ierr = ex_put_var_param(exoid, "M", nNodesetFuncs);
00402 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00403
00404 ierr = ex_put_nset_var_tab(exoid, nNodesets,
00405 nNodesetFuncs, &(nodesetFuncTruthTable[0]));
00406 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00407
00408 Array<std::string> nsFuncNames(nNodesetFuncs);
00409 Array<const char*> nsNameP;
00410
00411 for (int i=0; i<nNodesetFuncs; i++)
00412 {
00413 nsFuncNames[i] = pointScalarNames()[nsFuncs[i]];
00414 }
00415 getCharpp(nsFuncNames, nsNameP);
00416
00417 ierr = ex_put_var_names(exoid, "M", nNodesetFuncs, (char**)&(nsNameP[0]));
00418 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00419 }
00420
00421
00422 Array<std::string> nodalFuncNames(nNodalFuncs);
00423 Array<std::string> elemFuncNames(nElemFuncs);
00424 Array<const char*> nNameP;
00425 Array<const char*> eNameP;
00426
00427 for (int i=0; i<nNodalFuncs; i++)
00428 {
00429 nodalFuncNames[i] = pointScalarNames()[omnipresentNodalFuncs[i]];
00430 }
00431 getCharpp(nodalFuncNames, nNameP);
00432
00433 for (int i=0; i<nElemFuncs; i++)
00434 {
00435 elemFuncNames[i] = cellScalarNames()[omnipresentElemFuncs[i]];
00436 }
00437 getCharpp(elemFuncNames, eNameP);
00438
00439 if (nNodalFuncs > 0)
00440 {
00441 ierr = ex_put_var_names(exoid, "N", nNodalFuncs, (char**)&(nNameP[0]));
00442 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00443
00444 Array<double> funcVals;
00445 Array<int> nodeID(mesh().numCells(0));
00446 for (int i=0; i<mesh().numCells(0); i++) nodeID[i]=i;
00447
00448 for (int i=0; i<nNodalFuncs; i++)
00449 {
00450 int f = omnipresentNodalFuncs[i];
00451 pointScalarFields()[f]->getDataBatch(0, nodeID, tuple(f), funcVals);
00452 int t = 1;
00453 int numNodes = funcVals.size();
00454 ierr = ex_put_nodal_var(exoid, t, i+1, numNodes, &(funcVals[0]));
00455 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00456 }
00457
00458 for (int i=0; i<nNodesetFuncs; i++)
00459 {
00460 const Array<int>& ns = nsFuncNodesets[i];
00461 int fid = nsFuncs[i];
00462
00463 for (int s=0; s<ns.size(); s++)
00464 {
00465 const Array<int>& nodes = *(nodesForNodeset[ns[s]]);
00466 pointScalarFields()[fid]->getDataBatch(0, nodes, tuple(fid), funcVals);
00467 int t = 1;
00468 int numNodes = funcVals.size();
00469 int id = nsID[ns[s]];
00470 ierr = ex_put_nset_var(exoid, t, i+1, id, numNodes, &(funcVals[0]));
00471 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00472 }
00473 }
00474 }
00475
00476 if (nElemFuncs > 0)
00477 {
00478 ierr = ex_put_var_names(exoid, "E", nElemFuncs, (char**)&(eNameP[0]));
00479 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00480
00481 Array<double> funcVals;
00482 int dim = mesh().spatialDim();
00483 Array<int> elemID(mesh().numCells(dim));
00484 for (int i=0; i<mesh().numCells(dim); i++) elemID[i]=i;
00485
00486 for (int i=0; i<nElemFuncs; i++)
00487 {
00488 int f = omnipresentElemFuncs[i];
00489 cellScalarFields()[f]->getDataBatch(dim, elemID, tuple(f), funcVals);
00490 int t = 1;
00491 int numElems = funcVals.size();
00492 ierr = ex_put_elem_var(exoid, t, i+1, 1, numElems, &(funcVals[0]));
00493 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00494 }
00495 }
00496
00497
00498 #else
00499 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Exodus not enabled");
00500 #endif
00501
00502 }
00503
00504
00505
00506 std::string ExodusWriter::elemType(const CellType& type) const
00507 {
00508 switch(type)
00509 {
00510 case TriangleCell:
00511 return "TRIANGLE";
00512 case TetCell:
00513 return "TETRA";
00514 default:
00515 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "cell type=" << type << " cannot be used as a "
00516 "maximal-dimension cell in exodus");
00517 }
00518 return "NULL";
00519 }
00520
00521
00522
00523 void ExodusWriter::writeParallelInfo(const std::string& parfile) const
00524 {
00525 std::ofstream os(parfile.c_str());
00526
00527 int dim = mesh().spatialDim();
00528 int nCells = mesh().numCells(dim);
00529 int nPts = mesh().numCells(0);
00530
00531 os << myRank() << " " << nProc() << std::endl;
00532
00533 os << nPts << std::endl;
00534 for (int i=0; i<nPts; i++)
00535 {
00536 os << i << " " << mesh().mapLIDToGID(0,i)
00537 << " " << mesh().ownerProcID(0,i) << std::endl;
00538 }
00539
00540 os << nCells << std::endl;
00541 for (int c=0; c<nCells; c++)
00542 {
00543 os << c << " " << mesh().mapLIDToGID(dim,c)
00544 << " " << mesh().ownerProcID(dim,c) << std::endl;
00545 }
00546
00547 for (int i=0; i<comments().length(); i++)
00548 {
00549 os << "# " << comments()[i] << std::endl;
00550 }
00551 }
00552
00553
00554
00555
00556
00557 void ExodusWriter::findNodeSets(
00558 Array<CellFilter>& nodesetFilters,
00559 Array<int>& omnipresentFuncs,
00560 Array<RCP<Array<int> > >& funcsForNodeset,
00561 Array<RCP<Array<int> > >& nodesForNodeset,
00562 Array<int>& nsID,
00563 Array<int>& nNodesPerSet,
00564 Array<int>& nsNodePtr,
00565 RCP<Array<int> > allNodes
00566 ) const
00567 {
00568
00569 int verb = 0;
00570
00571 const Array<RCP<FieldBase> >& f = pointScalarFields();
00572 CellFilter maximal = new MaximalCellFilter();
00573
00574 nNodesPerSet.resize(0);
00575 nsNodePtr.resize(0);
00576 nsID.resize(0);
00577
00578 Map<CellFilter, RCP<Array<int> > > tmp;
00579
00580 for (int i=0; i<f.size(); i++)
00581 {
00582 const CellFilter& cf = f[i]->domain();
00583 if (cf==maximal)
00584 {
00585 SUNDANCE_MSG2(verb, "function #" << i << " is defined on all nodes");
00586 omnipresentFuncs.append(i);
00587 continue;
00588 }
00589 if (!tmp.containsKey(cf))
00590 {
00591 RCP<Array<int> > a = rcp(new Array<int>());
00592 tmp.put(cf, a);
00593 }
00594 SUNDANCE_MSG2(verb, "function #" << i << " is defined on CF " << cf);
00595 tmp[cf]->append(i);
00596 }
00597
00598 int nodesetID=1;
00599 int nodePtr=0;
00600 nodesetFilters.resize(0);
00601 funcsForNodeset.resize(0);
00602 nodesForNodeset.resize(0);
00603
00604 for (Map<CellFilter, RCP<Array<int> > >::const_iterator
00605 i=tmp.begin(); i!=tmp.end(); i++)
00606 {
00607 const CellFilter& cf = i->first;
00608 nodesetFilters.append(cf);
00609 funcsForNodeset.append(i->second);
00610 RCP<Array<int> > cells
00611 = cellSetToLIDArray(connectedNodeSet(cf, mesh()));
00612 nodesForNodeset.append(cells);
00613 int nn = cells->size();
00614 nNodesPerSet.append(nn);
00615 nsID.append(nodesetID++);
00616 nsNodePtr.append(nodePtr);
00617 nodePtr += nn;
00618 }
00619
00620 SUNDANCE_MSG2(verb, "node set IDs = " << nsID);
00621 SUNDANCE_MSG2(verb, "num nodes = " << nNodesPerSet);
00622 SUNDANCE_MSG2(verb, "node set pointers = " << nsNodePtr);
00623
00624
00625 int numNodes = nodePtr;
00626 allNodes->resize(numNodes);
00627
00628 int k=0;
00629 for (int i=0; i<nsID.size(); i++)
00630 {
00631 SUNDANCE_MSG2(verb, "node set " << i << " funcs = "
00632 << *funcsForNodeset[i]);
00633 SUNDANCE_MSG2(verb, "node set " << i
00634 << " nodes = " << *nodesForNodeset[i]);
00635 const Array<int>& myCells = *(nodesForNodeset[i]);
00636 for (int c=0; c<myCells.size(); c++)
00637 {
00638 (*allNodes)[k++] = myCells[c];
00639 }
00640 }
00641
00642 SUNDANCE_MSG2(verb, "all nodes = " << *allNodes);
00643 }
00644
00645
00646
00647 void ExodusWriter::findBlocks(
00648 Array<CellFilter>& blockFilters,
00649 Array<int>& omnipresentFuncs,
00650 Array<RCP<Array<int> > >& funcsForBlock,
00651 Array<RCP<Array<int> > >& elemsForBlock,
00652 Array<int>& blockIDs,
00653 Array<int>& nElemsPerBlock,
00654 Array<int>& elemBlockPtr,
00655 RCP<Array<int> > allElems
00656 ) const
00657 {
00658
00659 int verb=0;
00660 const Array<RCP<FieldBase> >& f = cellScalarFields();
00661 CellFilter maximal = new MaximalCellFilter();
00662
00663 nElemsPerBlock.resize(0);
00664 elemBlockPtr.resize(0);
00665 blockIDs.resize(0);
00666
00667 Map<CellFilter, RCP<Array<int> > > tmp;
00668
00669 for (int i=0; i<f.size(); i++)
00670 {
00671 const CellFilter& cf = f[i]->domain();
00672 if (cf==maximal)
00673 {
00674 SUNDANCE_MSG2(verb, "function #" << i << " is defined on all nodes");
00675 omnipresentFuncs.append(i);
00676 continue;
00677 }
00678 if (!tmp.containsKey(cf))
00679 {
00680 RCP<Array<int> > a = rcp(new Array<int>());
00681 tmp.put(cf, a);
00682 }
00683 SUNDANCE_MSG2(verb, "function #" << i << " is defined on CF " << cf);
00684 tmp[cf]->append(i);
00685 }
00686
00687 int blockID=1;
00688 int blockPtr=0;
00689 blockFilters.resize(0);
00690 funcsForBlock.resize(0);
00691 elemsForBlock.resize(0);
00692
00693 for (Map<CellFilter, RCP<Array<int> > >::const_iterator
00694 i=tmp.begin(); i!=tmp.end(); i++)
00695 {
00696 const CellFilter& cf = i->first;
00697 blockFilters.append(cf);
00698 funcsForBlock.append(i->second);
00699 RCP<Array<int> > cells
00700 = cellSetToLIDArray(cf.getCells(mesh()));
00701 elemsForBlock.append(cells);
00702 int nn = cells->size();
00703 nElemsPerBlock.append(nn);
00704 blockIDs.append(blockID++);
00705 elemBlockPtr.append(blockPtr);
00706 blockPtr += nn;
00707 }
00708
00709 SUNDANCE_MSG2(verb, "block IDs = " << blockIDs);
00710 SUNDANCE_MSG2(verb, "num elems = " << nElemsPerBlock);
00711 SUNDANCE_MSG2(verb, "block pointers = " << elemBlockPtr);
00712
00713
00714 int numElems = blockPtr;
00715 allElems->resize(numElems);
00716
00717 int k=0;
00718 for (int i=0; i<blockIDs.size(); i++)
00719 {
00720 SUNDANCE_MSG2(verb, "block " << i << " funcs = "
00721 << *funcsForBlock[i]);
00722 SUNDANCE_MSG2(verb, "block " << i
00723 << " elems = " << *elemsForBlock[i]);
00724 const Array<int>& myCells = *(elemsForBlock[i]);
00725 for (int c=0; c<myCells.size(); c++)
00726 {
00727 (*allElems)[k++] = myCells[c];
00728 }
00729 }
00730
00731 SUNDANCE_MSG2(verb, "all elems = " << *allElems);
00732 }
00733
00734 void ExodusWriter::getCharpp(const Array<std::string>& s, Array<const char*>& p) const
00735 {
00736 p.resize(s.size());
00737 for (int i=0; i<p.size(); i++) p[i] = s[i].c_str();
00738 }
00739
00740
