XdmfPartitioner.cpp 23 KB
Newer Older
Kenneth Leiter's avatar
Kenneth Leiter committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
/*****************************************************************************/
/*                                    XDMF                                   */
/*                       eXtensible Data Model and Format                    */
/*                                                                           */
/*  Id : XdmfPartitioner.cpp                                                 */
/*                                                                           */
/*  Author:                                                                  */
/*     Kenneth Leiter                                                        */
/*     kenneth.leiter@arl.army.mil                                           */
/*     US Army Research Laboratory                                           */
/*     Aberdeen Proving Ground, MD                                           */
/*                                                                           */
/*     Copyright @ 2011 US Army Research Laboratory                          */
/*     All Rights Reserved                                                   */
/*     See Copyright.txt for details                                         */
/*                                                                           */
/*     This software is distributed WITHOUT ANY WARRANTY; without            */
/*     even the implied warranty of MERCHANTABILITY or FITNESS               */
/*     FOR A PARTICULAR PURPOSE.  See the above copyright notice             */
/*     for more information.                                                 */
/*                                                                           */
/*****************************************************************************/

24 25 26 27
#ifndef BUILD_EXE

extern "C"
{
28
#include <metis.h>
29 30
}

31
#include <iostream>
32
#include <sstream>
33
#include "XdmfAttribute.hpp"
34 35
#include "XdmfAttributeCenter.hpp"
#include "XdmfAttributeType.hpp"
36
#include "XdmfGeometry.hpp"
37
#include "XdmfGeometryType.hpp"
38
#include "XdmfGridCollection.hpp"
39
#include "XdmfGridCollectionType.hpp"
40
#include "XdmfHeavyDataWriter.hpp"
41
#include "XdmfMap.hpp"
42
#include "XdmfPartitioner.hpp"
43
#include "XdmfSet.hpp"
44
#include "XdmfSetType.hpp"
45
#include "XdmfTopology.hpp"
46
#include "XdmfTopologyType.hpp"
47
#include "XdmfUnstructuredGrid.hpp"
48
#include "XdmfError.hpp"
49

50 51
#include "XdmfArrayType.hpp"

52
shared_ptr<XdmfPartitioner>
53
XdmfPartitioner::New()
54
{
55
  shared_ptr<XdmfPartitioner> p(new XdmfPartitioner());
56
  return p;
57 58
}

59 60 61 62 63 64 65 66
XdmfPartitioner::XdmfPartitioner()
{
}

XdmfPartitioner::~XdmfPartitioner()
{
}

67 68
void
XdmfPartitioner::ignore(const shared_ptr<const XdmfSet> set)
69 70 71 72
{
  mIgnoredSets.insert(set);
}

73 74
shared_ptr<XdmfGridCollection>
XdmfPartitioner::partition(const shared_ptr<XdmfUnstructuredGrid> gridToPartition,
75
                           const unsigned int numberOfPartitions,
76
                           const shared_ptr<XdmfHeavyDataWriter> heavyDataWriter) const
77
{
78
  // Make sure geometry and topology are non null
79
  if(!(gridToPartition->getGeometry() && gridToPartition->getTopology()))
80 81 82
    XdmfError::message(XdmfError::FATAL, 
                       "Current grid's geometry or topology is null in "
                       "XdmfPartitioner::partition");
83

84
  const shared_ptr<XdmfGeometry> geometry =
85
    gridToPartition->getGeometry();
86
  const shared_ptr<const XdmfGeometryType> geometryType =
87
    geometry->getType();
88
  const shared_ptr<XdmfTopology> topology =
89
    gridToPartition->getTopology();
90
  const shared_ptr<const XdmfTopologyType> topologyType =
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
    topology->getType();

  int metisElementType;
  unsigned int nodesPerElement;

  if(topologyType == XdmfTopologyType::Triangle() ||
     topologyType == XdmfTopologyType::Triangle_6()) {
    metisElementType = 1;
    nodesPerElement = 3;
  }
  else if(topologyType == XdmfTopologyType::Quadrilateral() ||
          topologyType == XdmfTopologyType::Quadrilateral_8()) {
    metisElementType = 4;
    nodesPerElement = 4;
  }
  else if(topologyType == XdmfTopologyType::Tetrahedron() ||
          topologyType == XdmfTopologyType::Tetrahedron_10()) {
    metisElementType = 2;
    nodesPerElement = 4;
  }
  else if(topologyType == XdmfTopologyType::Hexahedron() ||
          topologyType == XdmfTopologyType::Hexahedron_20() ||
          topologyType == XdmfTopologyType::Hexahedron_24() ||
          topologyType == XdmfTopologyType::Hexahedron_27() ||
115
          topologyType == XdmfTopologyType::Hexahedron_64() ||
116 117 118 119 120 121 122
          topologyType == XdmfTopologyType::Hexahedron_125() ||
          topologyType == XdmfTopologyType::Hexahedron_216() ||
          topologyType == XdmfTopologyType::Hexahedron_343() ||
          topologyType == XdmfTopologyType::Hexahedron_512() ||
          topologyType == XdmfTopologyType::Hexahedron_729() ||
          topologyType == XdmfTopologyType::Hexahedron_1000() ||
          topologyType == XdmfTopologyType::Hexahedron_1331()) {
123 124 125 126
    metisElementType = 3;
    nodesPerElement = 8;
  }
  else {
127 128 129 130
    XdmfError::message(XdmfError::FATAL, 
                       "Topology type is not 'Triangle', 'Quadrilateral', "
                       "'Tetrahedron', 'Hexahedron' in "
                       "XdmfPartition::partition");
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
  }

  bool releaseTopology = false;
  if(!topology->isInitialized()) {
    topology->read();
    releaseTopology = true;
  }

  int numElements = topology->getNumberElements();
  idxtype * metisConnectivity = new idxtype[nodesPerElement * numElements];
  unsigned int connectivityOffset = 0;
  idxtype * metisConnectivityPtr = metisConnectivity;
  for(int i=0; i<numElements; ++i) {
    topology->getValues(connectivityOffset,
                        metisConnectivityPtr,
                        nodesPerElement);
    connectivityOffset += topologyType->getNodesPerElement();
    metisConnectivityPtr += nodesPerElement;
  }

  int numNodes = geometry->getNumberPoints();

  // Need to remap connectivity for nonlinear elements so that metis handles
  // it properly.
  std::map<idxtype, idxtype> xdmfIdToMetisId;
  if(nodesPerElement != topologyType->getNodesPerElement()) {
    unsigned int index = 0;
    for (unsigned int i=0; i<numElements * nodesPerElement; ++i) {
      std::map<idxtype, idxtype>::const_iterator val =
        xdmfIdToMetisId.find(metisConnectivity[i]);
      if (val != xdmfIdToMetisId.end()) {
        metisConnectivity[i] = val->second;
      }
      else {
        // Haven't seen this id before, map to index and set to new id
        xdmfIdToMetisId[metisConnectivity[i]] = index;
        metisConnectivity[i] = index;
        index++;
      }
    }
    numNodes = index;
  }

  int startIndex = 0;
  int numCutEdges = 0;

  idxtype * elementsPartition = new idxtype[numElements];
  idxtype * nodesPartition = new idxtype[numNodes];

180
  /*
181 182 183 184 185 186 187 188 189
    METIS_PartMeshDual(&numElements,
    &numNodes,
    metisConnectivity,
    &metisElementType,
    &startIndex,
    (int*)&numberOfPartitions,
    &numCutEdges,
    elementsPartition,
    nodesPartition);
190
  */
191

192 193 194 195 196
  METIS_PartMeshNodal(&numElements,
                      &numNodes,
                      metisConnectivity,
                      &metisElementType,
                      &startIndex,
197
                      (int*)&numberOfPartitions,
198 199 200
                      &numCutEdges,
                      elementsPartition,
                      nodesPartition);
201

202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
  delete [] metisConnectivity;
  delete [] nodesPartition;

  // For each partition, map global to local node id
  std::vector<std::map<unsigned int, unsigned int> > globalToLocalNodeIdMap;
  // For each partition, list global element id.
  std::vector<std::vector<unsigned int> > globalElementIds;
  for(unsigned int i=0; i<numberOfPartitions; ++i) {
    std::map<unsigned int, unsigned int> nodeMap;
    globalToLocalNodeIdMap.push_back(nodeMap);
    std::vector<unsigned int> elementIds;
    globalElementIds.push_back(elementIds);
  }

  // Fill in globalNodeId for each partition
  unsigned int totalIndex = 0;
  for (int i=0; i<numElements; ++i) {
    unsigned int partitionId = elementsPartition[i];
    for (unsigned int j=0; j<topologyType->getNodesPerElement(); ++j) {
      unsigned int globalNodeId = topology->getValue<unsigned int>(totalIndex);
      if (globalToLocalNodeIdMap[partitionId].count(globalNodeId) == 0) {
        // Have not seen this node, need to add to map
        unsigned int size = globalToLocalNodeIdMap[partitionId].size();
        globalToLocalNodeIdMap[partitionId][globalNodeId] = size;
      }
      totalIndex++;
    }
    globalElementIds[partitionId].push_back(i);
  }
  delete [] elementsPartition;

  bool generateGlobalNodeIds = !gridToPartition->getAttribute("GlobalNodeId");

235
  shared_ptr<XdmfGridCollection> partitionedGrids =
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
    XdmfGridCollection::New();
  partitionedGrids->setType(XdmfGridCollectionType::Spatial());

  bool releaseGeometry = false;
  if(!geometry->isInitialized()) {
    geometry->read();
    releaseGeometry = true;
  }

  // Split geometry and topology into proper partitions
  for(unsigned int i=0; i<numberOfPartitions; ++i) {
    std::map<unsigned int, unsigned int> & currNodeMap =
      globalToLocalNodeIdMap[i];
    std::vector<unsigned int> & currElemIds = globalElementIds[i];

    if(currElemIds.size() > 0) {
      std::stringstream name;
      name << gridToPartition->getName() << "_" << i;

255
      shared_ptr<XdmfUnstructuredGrid> partitioned =
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303
        XdmfUnstructuredGrid::New();
      partitioned->setName(name.str());
      partitionedGrids->insert(partitioned);

      // Fill in geometry for this partition
      partitioned->getGeometry()->setType(geometryType);
      unsigned int numDimensions = geometryType->getDimensions();
      partitioned->getGeometry()->initialize(geometry->getArrayType(),
                                             currNodeMap.size() * numDimensions);

      for(std::map<unsigned int, unsigned int>::const_iterator iter =
            currNodeMap.begin();
          iter != currNodeMap.end();
          ++iter) {
        partitioned->getGeometry()->insert(iter->second * numDimensions,
                                           geometry,
                                           iter->first * numDimensions,
                                           numDimensions);
      }

      if(heavyDataWriter) {
        partitioned->getGeometry()->accept(heavyDataWriter);
        partitioned->getGeometry()->release();
      }

      // Fill in topology for this partition
      partitioned->getTopology()->setType(topologyType);
      partitioned->getTopology()->initialize(topology->getArrayType(),
                                             currElemIds.size() * topologyType->getNodesPerElement());
      unsigned int index = 0;
      for(std::vector<unsigned int>::const_iterator iter = currElemIds.begin();
          iter != currElemIds.end();
          ++iter) {
        // Translate these global node ids to local node ids
        for(unsigned int j=0; j<topologyType->getNodesPerElement(); ++j) {
          unsigned int globalNodeId =
            currNodeMap[topology->getValue<unsigned int>(*iter * topologyType->getNodesPerElement() + j)];
          partitioned->getTopology()->insert(index, &globalNodeId, 1);
          index++;
        }
      }

      if(heavyDataWriter) {
        partitioned->getTopology()->accept(heavyDataWriter);
        partitioned->getTopology()->release();
      }

      if (generateGlobalNodeIds) {
304
        shared_ptr<XdmfAttribute> globalNodeIds = XdmfAttribute::New();
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
        globalNodeIds->setName("GlobalNodeId");
        globalNodeIds->setType(XdmfAttributeType::GlobalId());
        globalNodeIds->setCenter(XdmfAttributeCenter::Node());
        globalNodeIds->initialize<unsigned int>(currNodeMap.size());
        for(std::map<unsigned int, unsigned int>::const_iterator iter =
              currNodeMap.begin();
            iter != currNodeMap.end();
            ++iter) {
          globalNodeIds->insert(iter->second, &iter->first, 1);
        }
        partitioned->insert(globalNodeIds);

        if(heavyDataWriter) {
          globalNodeIds->accept(heavyDataWriter);
          globalNodeIds->release();
        }
      }
    }
  }

  if(releaseGeometry) {
    gridToPartition->getGeometry()->release();
  }
  if(releaseTopology) {
    gridToPartition->getTopology()->release();
  }

  // Split attributes into proper partitions
  for(unsigned int i=0; i<gridToPartition->getNumberAttributes(); ++i) {
334
    shared_ptr<XdmfAttribute> currAttribute =
335 336 337 338 339 340 341 342 343 344 345 346
      gridToPartition->getAttribute(i);
    bool releaseAttribute = false;
    if(!currAttribute->isInitialized()) {
      currAttribute->read();
      releaseAttribute = true;
    }
    unsigned int partitionId = 0;
    for(unsigned int j=0; j<numberOfPartitions; ++j) {
      std::map<unsigned int, unsigned int> & currNodeMap =
        globalToLocalNodeIdMap[j];
      std::vector<unsigned int> & currElemIds = globalElementIds[j];
      if(currElemIds.size() > 0) {
347
        shared_ptr<XdmfUnstructuredGrid> partitioned =
348 349
          partitionedGrids->getUnstructuredGrid(partitionId);
        partitionId++;
350 351
        shared_ptr<XdmfAttribute> createdAttribute =
          shared_ptr<XdmfAttribute>();
352 353 354 355 356 357 358 359 360 361 362 363
        if(currAttribute->getCenter() == XdmfAttributeCenter::Grid()) {
          // Insert into each partition
          createdAttribute = currAttribute;
        }
        else if(currAttribute->getCenter() == XdmfAttributeCenter::Cell() ||
                currAttribute->getCenter() == XdmfAttributeCenter::Face() ||
                currAttribute->getCenter() == XdmfAttributeCenter::Edge()) {
          createdAttribute = XdmfAttribute::New();
          createdAttribute->setName(currAttribute->getName());
          createdAttribute->setCenter(currAttribute->getCenter());
          createdAttribute->setType(currAttribute->getType());
          unsigned int index = 0;
364 365
          const unsigned int numberComponents = 
	    currAttribute->getSize() / topology->getNumberElements();
366
          createdAttribute->initialize(currAttribute->getArrayType(),
367
                                       currElemIds.size() * numberComponents);
368 369 370 371 372 373
          for(std::vector<unsigned int>::const_iterator iter =
                currElemIds.begin();
              iter != currElemIds.end();
              ++iter) {
            createdAttribute->insert(index,
                                     currAttribute,
374 375 376
                                     *iter * numberComponents,
                                     numberComponents);
            index += numberComponents;
377 378 379 380 381 382 383
          }
        }
        else if(currAttribute->getCenter() == XdmfAttributeCenter::Node()) {
          createdAttribute = XdmfAttribute::New();
          createdAttribute->setName(currAttribute->getName());
          createdAttribute->setCenter(currAttribute->getCenter());
          createdAttribute->setType(currAttribute->getType());
384 385 386 387
          createdAttribute->initialize(currAttribute->getArrayType(), 
				       currNodeMap.size());
	  const unsigned int numberComponents = 
	    currAttribute->getSize() / geometry->getNumberPoints();
388 389 390 391
          for(std::map<unsigned int, unsigned int>::const_iterator iter =
                currNodeMap.begin();
              iter != currNodeMap.end();
              ++iter) {
392
            createdAttribute->insert(iter->second * numberComponents,
393
                                     currAttribute,
394 395
                                     iter->first * numberComponents,
                                     numberComponents);
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416
          }
        }
        if(createdAttribute) {
          partitioned->insert(createdAttribute);
          if(heavyDataWriter) {
            if(!createdAttribute->isInitialized()) {
              createdAttribute->read();
            }
            createdAttribute->accept(heavyDataWriter);
            createdAttribute->release();
          }
        }
      }
    }
    if(releaseAttribute) {
      currAttribute->release();
    }
  }

  // Split sets into proper partitions
  for(unsigned int i=0; i<gridToPartition->getNumberSets(); ++i) {
417
    shared_ptr<XdmfSet> currSet = gridToPartition->getSet(i);
418 419 420 421 422 423 424 425 426 427 428 429
    if(mIgnoredSets.find(currSet) == mIgnoredSets.end()) {
      bool releaseSet = false;
      if(!currSet->isInitialized()) {
        currSet->read();
        releaseSet = true;
      }
      unsigned int partitionId = 0;
      for(unsigned int j=0; j<numberOfPartitions; ++j) {
        std::map<unsigned int, unsigned int> & currNodeMap =
          globalToLocalNodeIdMap[j];
        std::vector<unsigned int> & currElemIds = globalElementIds[j];
        if(currElemIds.size() > 0) {
430
          shared_ptr<XdmfUnstructuredGrid> partitioned =
431 432
            partitionedGrids->getUnstructuredGrid(partitionId);
          partitionId++;
433 434

          shared_ptr<XdmfSet> partitionedSet = XdmfSet::New();
435
          std::vector<unsigned int> partitionedSetIndex;
436

437 438 439 440 441 442 443 444 445 446 447 448 449
          if(currSet->getType() == XdmfSetType::Cell() ||
             currSet->getType() == XdmfSetType::Face() ||
             currSet->getType() == XdmfSetType::Edge()) {
            for(unsigned int k=0; k<currSet->getSize(); ++k) {
              std::vector<unsigned int>::const_iterator val =
                std::find(currElemIds.begin(),
                          currElemIds.end(),
                          currSet->getValue<unsigned int>(k));
              if(val != currElemIds.end()) {
                unsigned int valToPush = val - currElemIds.begin();
                partitionedSet->pushBack(valToPush);
                partitionedSetIndex.push_back(k);
              }
450 451
            }
          }
452 453 454 455 456 457 458 459
          else if(currSet->getType() == XdmfSetType::Node()) {
            for(unsigned int k=0; k<currSet->getSize(); ++k) {
              std::map<unsigned int, unsigned int>::const_iterator val =
                currNodeMap.find(currSet->getValue<unsigned int>(k));
              if(val != currNodeMap.end()) {
                partitionedSet->pushBack(val->second);
                partitionedSetIndex.push_back(k);
              }
460 461
            }
          }
462

463 464 465
          if(partitionedSet->getSize() > 0) {

            for(unsigned int k=0; k<currSet->getNumberAttributes(); ++k) {
466
              const shared_ptr<XdmfAttribute> currAttribute =
467 468 469 470 471 472 473
                currSet->getAttribute(k);
              if(currAttribute->getCenter() == XdmfAttributeCenter::Node()) {
                bool releaseAttribute = false;
                if(!currAttribute->isInitialized()) {
                  currAttribute->read();
                  releaseAttribute = true;
                }
474
                const shared_ptr<XdmfAttribute> partitionedAttribute =
475 476 477 478 479 480 481
                  XdmfAttribute::New();
                partitionedAttribute->setCenter(currAttribute->getCenter());
                partitionedAttribute->setName(currAttribute->getName());
                partitionedAttribute->setType(currAttribute->getType());
                partitionedAttribute->initialize(currAttribute->getArrayType(),
                                                 partitionedSetIndex.size());
                for(unsigned int l=0; l<partitionedSetIndex.size(); ++l) {
482
                  partitionedAttribute->insert(l,
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500
                                               currAttribute,
                                               partitionedSetIndex[l]);
                }
                partitionedSet->insert(partitionedAttribute);
                if(heavyDataWriter) {
                  partitionedAttribute->accept(heavyDataWriter);
                  partitionedAttribute->release();
                }
              }
            }

            partitioned->insert(partitionedSet);
            partitionedSet->setName(currSet->getName());
            partitionedSet->setType(currSet->getType());
            if(heavyDataWriter) {
              partitionedSet->accept(heavyDataWriter);
              partitionedSet->release();
            }
501 502 503
          }
        }
      }
504 505 506
      if(releaseSet) {
        currSet->release();
      }
507 508 509 510
    }
  }

  // Add XdmfMap to map boundary nodes between partitions
511
  std::vector<shared_ptr<XdmfAttribute> > globalNodeIds;
512 513 514
  for(unsigned int i=0;
      i<partitionedGrids->getNumberUnstructuredGrids();
      ++i) {
515
    shared_ptr<XdmfAttribute> globalNodeId =
516
      partitionedGrids->getUnstructuredGrid(i)->getAttribute("GlobalNodeId");
517 518 519 520
    if(!globalNodeId->isInitialized()) {
      globalNodeId->read();
    }
    globalNodeIds.push_back(globalNodeId);
521 522
  }

523
  std::vector<shared_ptr<XdmfMap> > maps = XdmfMap::New(globalNodeIds);
524 525 526
  for(unsigned int i=0;
      i<partitionedGrids->getNumberUnstructuredGrids();
      ++i) {
527
    shared_ptr<XdmfMap> map = maps[i];
528 529
    map->setName("Subdomain Boundary");
    partitionedGrids->getUnstructuredGrid(i)->insert(map);
530 531 532 533 534 535 536 537
    if(heavyDataWriter) {
      globalNodeIds[i]->release();
      map->accept(heavyDataWriter);
      map->release();
    }
  }

  return partitionedGrids;
538 539
}

540 541
#else

542
#include <iostream>
543 544 545 546
#include <sstream>
#include "XdmfDomain.hpp"
#include "XdmfGridCollection.hpp"
#include "XdmfHDF5Writer.hpp"
547
#include "XdmfPartitioner.hpp"
548 549
#include "XdmfReader.hpp"
#include "XdmfWriter.hpp"
550 551 552

/**
 * XdmfPartitioner is a command line utility for partitioning Xdmf grids.
553 554
 * The XdmfPartitioner uses the metis library to partition Triangular,
 * Quadrilateral, Tetrahedral, and Hexahedral XdmfGridUnstructureds.
555 556
 *
 * Usage:
557 558
 *     XdmfPartitioner <path-of-file-to-partition> <num-partitions>
 *                     (Optional: <path-to-output-file>)
559 560
 *
 */
561
int main(int argc, char* argv[])
562
{
563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602
  std::string usage = "Partitions an Xdmf grid using the metis library: \n \n "
    "Usage: \n \n   XdmfPartitioner <path-of-file-to-partition> "
    "<num-partitions> (Optional: <path-to-output-file>)";
  std::string meshName = "";

  if (argc < 3) {
    std::cout << usage << std::endl;
    return 1;
  }

  FILE * refFile = fopen(argv[1], "r");
  if (refFile) {
    // Success
    meshName = argv[1];
    fclose(refFile);
  }
  else {
    std::cout << "Cannot open file: " << argv[1] << std::endl;
    return 1;
  }

  unsigned int numPartitions = atoi(argv[2]);

  if (argc >= 4) {
    meshName = argv[3];
  }

  if(meshName.find_last_of("/\\") != std::string::npos) {
    meshName = meshName.substr(meshName.find_last_of("/\\")+1,
                               meshName.length());
  }

  if (meshName.rfind(".") != std::string::npos) {
    meshName = meshName.substr(0, meshName.rfind("."));
  }

  if(argc < 4) {
    meshName = meshName + "-partitioned";
  }

603 604 605
  shared_ptr<XdmfReader> reader = XdmfReader::New();
  shared_ptr<XdmfDomain> domain =
    shared_dynamic_cast<XdmfDomain>(reader->read(argv[1]));
606 607 608 609 610 611 612 613

  if(domain->getNumberUnstructuredGrids() <= 0) {
    std::cout << "No grids to partition" << std::endl;
    return 1;
  }

  std::stringstream heavyFileName;
  heavyFileName << meshName << ".h5";
614
  shared_ptr<XdmfHDF5Writer> heavyDataWriter =
615 616
    XdmfHDF5Writer::New(heavyFileName.str());

617 618
  shared_ptr<XdmfPartitioner> partitioner = XdmfPartitioner::New();
  shared_ptr<XdmfGridCollection> toWrite =
619 620 621 622
    partitioner->partition(domain->getUnstructuredGrid(0),
                           numPartitions,
                           heavyDataWriter);

623
  shared_ptr<XdmfDomain> newDomain = XdmfDomain::New();
624 625 626 627
  newDomain->insert(toWrite);

  std::stringstream xmlFileName;
  xmlFileName << meshName << ".xmf";
628 629
  shared_ptr<XdmfWriter> writer = XdmfWriter::New(xmlFileName.str(),
                                                  heavyDataWriter);
630 631 632
  newDomain->accept(writer);

  std::cout << "Wrote: " << xmlFileName.str() << std::endl;
633
}
634

635
#endif // BUILD_EXE