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 Tabs tab0(0);
00331 Tabs tab1;
00332 int verb=3;
00333
00334
00335 int nNodalFuncs = omnipresentNodalFuncs().size();
00336 int nElemFuncs = omnipresentElemFuncs().size();
00337
00338 int nNodesetFuncs = pointScalarFields().size() - nNodalFuncs;
00339
00340 int nNodesets = funcsForNodeset.size();
00341
00342 PLAYA_ROOT_MSG1(verb, tab0 << "ExodusWriter::writeFields()");
00343 PLAYA_ROOT_MSG2(verb, tab1 << "nNodalFuncs = " << nNodalFuncs);
00344 PLAYA_ROOT_MSG2(verb, tab1 << "nElemFuncs = " << nElemFuncs);
00345 PLAYA_ROOT_MSG2(verb, tab1 << "nNodesetFuncs = " << nNodesetFuncs);
00346 PLAYA_ROOT_MSG2(verb, tab1 << "nNodesets = " << nNodesets);
00347
00348
00349 Set<int> nsFuncSet;
00350 Map<int, Array<int> > funcToNSMap;
00351
00352 for (int i=0; i<nNodesets; i++)
00353 {
00354 const Array<int>& f = *(funcsForNodeset[i]);
00355 for (int j=0; j<f.size(); j++)
00356 {
00357 nsFuncSet.put(f[j]);
00358 if (funcToNSMap.containsKey(f[j]))
00359 {
00360 funcToNSMap[f[j]].append(i);
00361 }
00362 else
00363 {
00364 funcToNSMap.put(f[j], tuple(i));
00365 }
00366 }
00367 }
00368 Array<int> nsFuncs = nsFuncSet.elements();
00369 TEUCHOS_TEST_FOR_EXCEPT(nsFuncs.size() != nNodesetFuncs);
00370
00371 Map<int, int > funcIDToNSFuncIndex;
00372 for (int i=0; i<nNodesetFuncs; i++) funcIDToNSFuncIndex.put(nsFuncs[i],i);
00373
00374 Array<Array<int> > nsFuncNodesets(nsFuncs.size());
00375 for (int i=0; i<nNodesetFuncs; i++)
00376 {
00377 nsFuncNodesets[i] = funcToNSMap.get(nsFuncs[i]);
00378 }
00379
00380 Array<int> nodesetFuncTruthTable(nNodesetFuncs * nNodesets, 0);
00381 for (int i=0; i<nNodesetFuncs; i++)
00382 {
00383 for (int j=0; j<nsFuncNodesets[i].size(); j++)
00384 {
00385 int ns = nsFuncNodesets[i][j];
00386 nodesetFuncTruthTable[ns*nNodesetFuncs + i] = 1;
00387 }
00388
00389 nsFuncNodesets[i] = funcToNSMap.get(nsFuncs[i]);
00390 }
00391
00392
00393
00394 int ierr;
00395
00396 if (nNodalFuncs > 0)
00397 {
00398 ierr = ex_put_var_param(exoid, "N", nNodalFuncs);
00399 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00400 }
00401
00402 if (nElemFuncs > 0)
00403 {
00404 ierr = ex_put_var_param(exoid, "E", nElemFuncs);
00405 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00406 }
00407
00408
00409 if (nNodesets > 0)
00410 {
00411 ierr = ex_put_var_param(exoid, "M", nNodesetFuncs);
00412 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00413
00414 ierr = ex_put_nset_var_tab(exoid, nNodesets,
00415 nNodesetFuncs, &(nodesetFuncTruthTable[0]));
00416 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00417
00418 Array<std::string> nsFuncNames(nNodesetFuncs);
00419 Array<const char*> nsNameP;
00420
00421 for (int i=0; i<nNodesetFuncs; i++)
00422 {
00423 nsFuncNames[i] = pointScalarNames()[nsFuncs[i]];
00424 }
00425 getCharpp(nsFuncNames, nsNameP);
00426
00427 ierr = ex_put_var_names(exoid, "M", nNodesetFuncs, (char**)&(nsNameP[0]));
00428 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00429 }
00430
00431
00432 Array<std::string> nodalFuncNames(nNodalFuncs);
00433 Array<std::string> elemFuncNames(nElemFuncs);
00434 Array<const char*> nNameP;
00435 Array<const char*> eNameP;
00436
00437 for (int i=0; i<nNodalFuncs; i++)
00438 {
00439 nodalFuncNames[i] = pointScalarNames()[omnipresentNodalFuncs[i]];
00440 }
00441 getCharpp(nodalFuncNames, nNameP);
00442
00443 for (int i=0; i<nElemFuncs; i++)
00444 {
00445 elemFuncNames[i] = cellScalarNames()[omnipresentElemFuncs[i]];
00446 }
00447 getCharpp(elemFuncNames, eNameP);
00448
00449 if (nNodalFuncs > 0)
00450 {
00451 ierr = ex_put_var_names(exoid, "N", nNodalFuncs, (char**)&(nNameP[0]));
00452 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00453
00454 Array<double> funcVals;
00455 Array<int> nodeID(mesh().numCells(0));
00456 for (int i=0; i<mesh().numCells(0); i++) nodeID[i]=i;
00457
00458 for (int i=0; i<nNodalFuncs; i++)
00459 {
00460 int f = omnipresentNodalFuncs[i];
00461 pointScalarFields()[f]->getDataBatch(0, nodeID, tuple(f), funcVals);
00462 int t = 1;
00463 int numNodes = funcVals.size();
00464 ierr = ex_put_nodal_var(exoid, t, i+1, numNodes, &(funcVals[0]));
00465 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00466 }
00467
00468 for (int i=0; i<nNodesetFuncs; i++)
00469 {
00470 const Array<int>& ns = nsFuncNodesets[i];
00471 int fid = nsFuncs[i];
00472
00473 for (int s=0; s<ns.size(); s++)
00474 {
00475 const Array<int>& nodes = *(nodesForNodeset[ns[s]]);
00476 pointScalarFields()[fid]->getDataBatch(0, nodes, tuple(fid), funcVals);
00477 int t = 1;
00478 int numNodes = funcVals.size();
00479 int id = nsID[ns[s]];
00480 ierr = ex_put_nset_var(exoid, t, i+1, id, numNodes, &(funcVals[0]));
00481 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00482 }
00483 }
00484 }
00485
00486 if (nElemFuncs > 0)
00487 {
00488 ierr = ex_put_var_names(exoid, "E", nElemFuncs, (char**)&(eNameP[0]));
00489 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00490
00491 Array<double> funcVals;
00492 int dim = mesh().spatialDim();
00493 Array<int> elemID(mesh().numCells(dim));
00494 for (int i=0; i<mesh().numCells(dim); i++) elemID[i]=i;
00495
00496 for (int i=0; i<nElemFuncs; i++)
00497 {
00498 int f = omnipresentElemFuncs[i];
00499 cellScalarFields()[f]->getDataBatch(dim, elemID, tuple(f), funcVals);
00500 int t = 1;
00501 int numElems = funcVals.size();
00502 ierr = ex_put_elem_var(exoid, t, i+1, 1, numElems, &(funcVals[0]));
00503 TEUCHOS_TEST_FOR_EXCEPT(ierr < 0);
00504 }
00505 }
00506
00507
00508 #else
00509 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Exodus not enabled");
00510 #endif
00511
00512 }
00513
00514
00515
00516 std::string ExodusWriter::elemType(const CellType& type) const
00517 {
00518 switch(type)
00519 {
00520 case TriangleCell:
00521 return "TRIANGLE";
00522 case TetCell:
00523 return "TETRA";
00524 default:
00525 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "cell type=" << type << " cannot be used as a "
00526 "maximal-dimension cell in exodus");
00527 }
00528 return "NULL";
00529 }
00530
00531
00532
00533 void ExodusWriter::writeParallelInfo(const std::string& parfile) const
00534 {
00535 std::ofstream os(parfile.c_str());
00536
00537 int dim = mesh().spatialDim();
00538 int nCells = mesh().numCells(dim);
00539 int nPts = mesh().numCells(0);
00540
00541 os << myRank() << " " << nProc() << std::endl;
00542
00543 os << nPts << std::endl;
00544 for (int i=0; i<nPts; i++)
00545 {
00546 os << i << " " << mesh().mapLIDToGID(0,i)
00547 << " " << mesh().ownerProcID(0,i) << std::endl;
00548 }
00549
00550 os << nCells << std::endl;
00551 for (int c=0; c<nCells; c++)
00552 {
00553 os << c << " " << mesh().mapLIDToGID(dim,c)
00554 << " " << mesh().ownerProcID(dim,c) << std::endl;
00555 }
00556
00557 for (int i=0; i<comments().length(); i++)
00558 {
00559 os << "# " << comments()[i] << std::endl;
00560 }
00561 }
00562
00563
00564
00565
00566
00567 void ExodusWriter::findNodeSets(
00568 Array<CellFilter>& nodesetFilters,
00569 Array<int>& omnipresentFuncs,
00570 Array<RCP<Array<int> > >& funcsForNodeset,
00571 Array<RCP<Array<int> > >& nodesForNodeset,
00572 Array<int>& nsID,
00573 Array<int>& nNodesPerSet,
00574 Array<int>& nsNodePtr,
00575 RCP<Array<int> > allNodes
00576 ) const
00577 {
00578
00579 int verb = 2;
00580
00581 const Array<RCP<FieldBase> >& f = pointScalarFields();
00582 CellFilter maximal = new MaximalCellFilter();
00583
00584 nNodesPerSet.resize(0);
00585 nsNodePtr.resize(0);
00586 nsID.resize(0);
00587
00588 Map<CellFilter, RCP<Array<int> > > tmp;
00589
00590 for (int i=0; i<f.size(); i++)
00591 {
00592 const CellFilter& cf = f[i]->domain();
00593 if (cf==maximal)
00594 {
00595 SUNDANCE_MSG2(verb, "function #" << i << " is defined on all nodes");
00596 omnipresentFuncs.append(i);
00597 continue;
00598 }
00599 if (!tmp.containsKey(cf))
00600 {
00601 RCP<Array<int> > a = rcp(new Array<int>());
00602 tmp.put(cf, a);
00603 }
00604 SUNDANCE_MSG2(verb, "function #" << i << " is defined on CF " << cf);
00605 tmp[cf]->append(i);
00606 }
00607
00608 int nodesetID=1;
00609 int nodePtr=0;
00610 nodesetFilters.resize(0);
00611 funcsForNodeset.resize(0);
00612 nodesForNodeset.resize(0);
00613
00614 for (Map<CellFilter, RCP<Array<int> > >::const_iterator
00615 i=tmp.begin(); i!=tmp.end(); i++)
00616 {
00617 const CellFilter& cf = i->first;
00618 nodesetFilters.append(cf);
00619 funcsForNodeset.append(i->second);
00620 RCP<Array<int> > cells
00621 = cellSetToLIDArray(connectedNodeSet(cf, mesh()));
00622 nodesForNodeset.append(cells);
00623 int nn = cells->size();
00624 nNodesPerSet.append(nn);
00625 nsID.append(nodesetID++);
00626 nsNodePtr.append(nodePtr);
00627 nodePtr += nn;
00628 }
00629
00630 SUNDANCE_MSG2(verb, "node set IDs = " << nsID);
00631 SUNDANCE_MSG2(verb, "num nodes = " << nNodesPerSet);
00632 SUNDANCE_MSG2(verb, "node set pointers = " << nsNodePtr);
00633
00634
00635 int numNodes = nodePtr;
00636 allNodes->resize(numNodes);
00637
00638 int k=0;
00639 for (int i=0; i<nsID.size(); i++)
00640 {
00641 SUNDANCE_MSG4(verb, "node set " << i << " funcs = "
00642 << *funcsForNodeset[i]);
00643 SUNDANCE_MSG4(verb, "node set " << i
00644 << " nodes = " << *nodesForNodeset[i]);
00645 const Array<int>& myCells = *(nodesForNodeset[i]);
00646 for (int c=0; c<myCells.size(); c++)
00647 {
00648 (*allNodes)[k++] = myCells[c];
00649 }
00650 }
00651
00652 SUNDANCE_MSG4(verb, "all nodes = " << *allNodes);
00653 }
00654
00655
00656
00657 void ExodusWriter::findBlocks(
00658 Array<CellFilter>& blockFilters,
00659 Array<int>& omnipresentFuncs,
00660 Array<RCP<Array<int> > >& funcsForBlock,
00661 Array<RCP<Array<int> > >& elemsForBlock,
00662 Array<int>& blockIDs,
00663 Array<int>& nElemsPerBlock,
00664 Array<int>& elemBlockPtr,
00665 RCP<Array<int> > allElems
00666 ) const
00667 {
00668
00669 int verb=0;
00670 const Array<RCP<FieldBase> >& f = cellScalarFields();
00671 CellFilter maximal = new MaximalCellFilter();
00672
00673 nElemsPerBlock.resize(0);
00674 elemBlockPtr.resize(0);
00675 blockIDs.resize(0);
00676
00677 Map<CellFilter, RCP<Array<int> > > tmp;
00678
00679 for (int i=0; i<f.size(); i++)
00680 {
00681 const CellFilter& cf = f[i]->domain();
00682 if (cf==maximal)
00683 {
00684 SUNDANCE_MSG2(verb, "function #" << i << " is defined on all nodes");
00685 omnipresentFuncs.append(i);
00686 continue;
00687 }
00688 if (!tmp.containsKey(cf))
00689 {
00690 RCP<Array<int> > a = rcp(new Array<int>());
00691 tmp.put(cf, a);
00692 }
00693 SUNDANCE_MSG2(verb, "function #" << i << " is defined on CF " << cf);
00694 tmp[cf]->append(i);
00695 }
00696
00697 int blockID=1;
00698 int blockPtr=0;
00699 blockFilters.resize(0);
00700 funcsForBlock.resize(0);
00701 elemsForBlock.resize(0);
00702
00703 for (Map<CellFilter, RCP<Array<int> > >::const_iterator
00704 i=tmp.begin(); i!=tmp.end(); i++)
00705 {
00706 const CellFilter& cf = i->first;
00707 blockFilters.append(cf);
00708 funcsForBlock.append(i->second);
00709 RCP<Array<int> > cells
00710 = cellSetToLIDArray(cf.getCells(mesh()));
00711 elemsForBlock.append(cells);
00712 int nn = cells->size();
00713 nElemsPerBlock.append(nn);
00714 blockIDs.append(blockID++);
00715 elemBlockPtr.append(blockPtr);
00716 blockPtr += nn;
00717 }
00718
00719 SUNDANCE_MSG2(verb, "block IDs = " << blockIDs);
00720 SUNDANCE_MSG2(verb, "num elems = " << nElemsPerBlock);
00721 SUNDANCE_MSG2(verb, "block pointers = " << elemBlockPtr);
00722
00723
00724 int numElems = blockPtr;
00725 allElems->resize(numElems);
00726
00727 int k=0;
00728 for (int i=0; i<blockIDs.size(); i++)
00729 {
00730 SUNDANCE_MSG2(verb, "block " << i << " funcs = "
00731 << *funcsForBlock[i]);
00732 SUNDANCE_MSG2(verb, "block " << i
00733 << " elems = " << *elemsForBlock[i]);
00734 const Array<int>& myCells = *(elemsForBlock[i]);
00735 for (int c=0; c<myCells.size(); c++)
00736 {
00737 (*allElems)[k++] = myCells[c];
00738 }
00739 }
00740
00741 SUNDANCE_MSG2(verb, "all elems = " << *allElems);
00742 }
00743
00744 void ExodusWriter::getCharpp(const Array<std::string>& s, Array<const char*>& p) const
00745 {
00746 p.resize(s.size());
00747 for (int i=0; i<p.size(); i++) p[i] = s[i].c_str();
00748 }
00749
00750
