XdmfPartitioner.cpp 16.5 KB
Newer Older
1
2
3
4
#ifndef BUILD_EXE

extern "C"
{
5
	#include <metis.h>
6
7
}

8
#include <iostream>
9
#include <sstream>
10
#include "XdmfAttribute.hpp"
11
12
#include "XdmfAttributeCenter.hpp"
#include "XdmfAttributeType.hpp"
13
#include "XdmfGeometry.hpp"
14
#include "XdmfGeometryType.hpp"
15
#include "XdmfGridUnstructured.hpp"
16
#include "XdmfGridCollection.hpp"
17
#include "XdmfGridCollectionType.hpp"
18
#include "XdmfHDF5Writer.hpp"
19
#include "XdmfMap.hpp"
20
#include "XdmfPartitioner.hpp"
21
#include "XdmfSet.hpp"
22
#include "XdmfSetType.hpp"
23
#include "XdmfTopology.hpp"
24
#include "XdmfTopologyType.hpp"
25

26
27
28
29
30
31
boost::shared_ptr<XdmfPartitioner> XdmfPartitioner::New()
{
	boost::shared_ptr<XdmfPartitioner> p(new XdmfPartitioner());
	return p;
}

32
33
34
35
36
37
38
39
XdmfPartitioner::XdmfPartitioner()
{
}

XdmfPartitioner::~XdmfPartitioner()
{
}

40
boost::shared_ptr<XdmfGridCollection> XdmfPartitioner::partition(const boost::shared_ptr<XdmfGridUnstructured> gridToPartition, const unsigned int numberOfPartitions,
41
		const boost::shared_ptr<XdmfHDF5Writer> heavyDataWriter) const
42
{
43
44
45
	// Make sure geometry and topology are non null
	assert(gridToPartition->getGeometry() && gridToPartition->getTopology());

46
	int metisElementType;
47
	unsigned int nodesPerElement;
48

49
	boost::shared_ptr<const XdmfTopologyType> topologyType = gridToPartition->getTopology()->getType();
50
51
52
	if(topologyType == XdmfTopologyType::Triangle() || topologyType == XdmfTopologyType::Triangle_6())
	{
		metisElementType = 1;
53
		nodesPerElement = 3;
54
55
56
57
	}
	else if(topologyType == XdmfTopologyType::Quadrilateral() || topologyType == XdmfTopologyType::Quadrilateral_8())
	{
		metisElementType = 4;
58
		nodesPerElement = 4;
59
60
61
62
	}
	else if(topologyType == XdmfTopologyType::Tetrahedron() || topologyType == XdmfTopologyType::Tetrahedron_10())
	{
		metisElementType = 2;
63
		nodesPerElement = 4;
64
65
66
67
68
69
	}
	else if(topologyType == XdmfTopologyType::Hexahedron() || topologyType == XdmfTopologyType::Hexahedron_20() ||
			topologyType == XdmfTopologyType::Hexahedron_24() || topologyType == XdmfTopologyType::Hexahedron_27() ||
			topologyType == XdmfTopologyType::Hexahedron_64())
	{
		metisElementType = 3;
70
		nodesPerElement = 8;
71
72
73
74
75
76
	}
	else
	{
		assert(false);
	}

77
	bool releaseTopology = false;
78
	if(!gridToPartition->getTopology()->isInitialized())
79
	{
80
		gridToPartition->getTopology()->read();
81
		releaseTopology = true;
82
83
84
	}

	int numElements = gridToPartition->getTopology()->getNumberElements();
85
	idxtype * metisConnectivity = new idxtype[nodesPerElement * numElements];
86
	for(int i=0; i<numElements; ++i)
87
	{
88
		gridToPartition->getTopology()->getValues(i*topologyType->getNodesPerElement(), metisConnectivity + i * nodesPerElement, nodesPerElement);
89
90
91
92
93
94
	}

	int numNodes = gridToPartition->getGeometry()->getNumberPoints();

	// Need to remap connectivity for nonlinear elements so that metis handles it properly.
	std::map<idxtype, idxtype> xdmfIdToMetisId;
95
	if(nodesPerElement != topologyType->getNodesPerElement())
96
	{
97
		unsigned int index = 0;
98
		for (unsigned int i=0; i<numElements * nodesPerElement; ++i)
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
		{
			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];

	std::cout << "Entered METIS" << std::endl;
	METIS_PartMeshDual(&numElements, &numNodes, metisConnectivity, &metisElementType, &startIndex, (int*)&numberOfPartitions, &numCutEdges, elementsPartition, nodesPartition);
	//METIS_PartMeshNodal(&numElements, &numNodes, metisConnectivity, &metisElementType, &startIndex, &numPartitions, &numCutEdges, elementsPartition, nodesPartition);
	std::cout << "Exited METIS" << std::endl;

	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;
144
	for (int i=0; i<numElements; ++i)
145
146
	{
		unsigned int partitionId = elementsPartition[i];
147
		for (unsigned int j=0; j<topologyType->getNodesPerElement(); ++j)
148
		{
149
			unsigned int globalNodeId = gridToPartition->getTopology()->getValue<unsigned int>(totalIndex);
150
151
152
			if (globalToLocalNodeIdMap[partitionId].count(globalNodeId) == 0)
			{
				// Have not seen this node, need to add to map
153
				unsigned int size = globalToLocalNodeIdMap[partitionId].size();
154
155
156
157
158
159
				globalToLocalNodeIdMap[partitionId][globalNodeId] = size;
			}
			totalIndex++;
		}
		globalElementIds[partitionId].push_back(i);
	}
160
	delete [] elementsPartition;
161

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

	boost::shared_ptr<XdmfGridCollection> partitionedGrids = XdmfGridCollection::New();
165
	partitionedGrids->setType(XdmfGridCollectionType::Spatial());
166

167
	bool releaseGeometry = false;
168
	if(!gridToPartition->getGeometry()->isInitialized())
169
	{
170
		gridToPartition->getGeometry()->read();
171
		releaseGeometry = true;
172
173
	}

174
	// Split geometry and topology into proper partitions
175
176
	for(unsigned int i=0; i<numberOfPartitions; ++i)
	{
177
178
		std::map<unsigned int, unsigned int> & currNodeMap = globalToLocalNodeIdMap[i];
		std::vector<unsigned int> & currElemIds = globalElementIds[i];
179
180
181
182
183
184

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

185
			boost::shared_ptr<XdmfGridUnstructured> partitioned = XdmfGridUnstructured::New();
186
187
188
189
			partitioned->setName(name.str());
			partitionedGrids->insert(partitioned);

			// Fill in geometry for this partition
190
191
			partitioned->getGeometry()->setType(gridToPartition->getGeometry()->getType());
			unsigned int numDimensions = partitioned->getGeometry()->getType()->getDimensions();
192
			partitioned->getGeometry()->initialize(gridToPartition->getGeometry()->getArrayType(), currNodeMap.size() * numDimensions);
193
194
195

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

			if(heavyDataWriter)
			{
201
202
				partitioned->getGeometry()->accept(heavyDataWriter);
				partitioned->getGeometry()->release();
203
204
205
			}

			// Fill in topology for this partition
206
			partitioned->getTopology()->setType(gridToPartition->getTopology()->getType());
207
			partitioned->getTopology()->initialize(gridToPartition->getTopology()->getArrayType(), currElemIds.size() * topologyType->getNodesPerElement());
208
209
210
211
			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
212
				for(unsigned int j=0; j<topologyType->getNodesPerElement(); ++j)
213
				{
214
215
					unsigned int globalNodeId = currNodeMap[gridToPartition->getTopology()->getValue<unsigned int>(*iter * topologyType->getNodesPerElement() + j)];
					partitioned->getTopology()->insert(index, &globalNodeId, 1);
216
217
218
219
220
221
					index++;
				}
			}

			if(heavyDataWriter)
			{
222
223
				partitioned->getTopology()->accept(heavyDataWriter);
				partitioned->getTopology()->release();
224
			}
225

226
227
			if (generateGlobalNodeIds)
			{
228
				boost::shared_ptr<XdmfAttribute> globalNodeIds = XdmfAttribute::New();
229
				globalNodeIds->setName("GlobalNodeId");
230
231
				globalNodeIds->setType(XdmfAttributeType::GlobalId());
				globalNodeIds->setCenter(XdmfAttributeCenter::Node());
232
233
				globalNodeIds->initialize<unsigned int>(currNodeMap.size());
				for(std::map<unsigned int, unsigned int>::const_iterator iter = currNodeMap.begin(); iter != currNodeMap.end(); ++iter)
234
				{
235
					globalNodeIds->insert(iter->second, &iter->first, 1);
236
237
				}
				partitioned->insert(globalNodeIds);
238

239
240
				if(heavyDataWriter)
				{
241
242
					globalNodeIds->accept(heavyDataWriter);
					globalNodeIds->release();
243
244
245
246
				}
			}
		}
	}
247

248
249
250
251
252
253
254
255
	if(releaseGeometry)
	{
		gridToPartition->getGeometry()->release();
	}
	if(releaseTopology)
	{
		gridToPartition->getTopology()->release();
	}
256

257
	// Split attributes into proper partitions
258
	for(unsigned int i=0; i<gridToPartition->getNumberAttributes(); ++i)
259
260
	{
		boost::shared_ptr<XdmfAttribute> currAttribute = gridToPartition->getAttribute(i);
261
		bool releaseAttribute = false;
262
		if(!currAttribute->isInitialized())
263
		{
264
			currAttribute->read();
265
			releaseAttribute = true;
266
267
268
269
270
271
272
273
		}
		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)
			{
274
				boost::shared_ptr<XdmfGridUnstructured> partitioned = partitionedGrids->getGridUnstructured(partitionId);
275
				partitionId++;
276
				boost::shared_ptr<XdmfAttribute> createdAttribute = boost::shared_ptr<XdmfAttribute>();
277
				if(currAttribute->getCenter() == XdmfAttributeCenter::Grid())
278
				{
279
280
					// Insert into each partition
					createdAttribute = currAttribute;
281
				}
282
283
284
				else if(currAttribute->getCenter() == XdmfAttributeCenter::Cell() ||
						currAttribute->getCenter() == XdmfAttributeCenter::Face() ||
						currAttribute->getCenter() == XdmfAttributeCenter::Edge())
285
286
287
				{
					createdAttribute = XdmfAttribute::New();
					createdAttribute->setName(currAttribute->getName());
288
289
					createdAttribute->setCenter(currAttribute->getCenter());
					createdAttribute->setType(currAttribute->getType());
290
					unsigned int index = 0;
291
					unsigned int numValsPerComponent = currAttribute->getSize() / gridToPartition->getTopology()->getNumberElements();
292
					createdAttribute->initialize(currAttribute->getArrayType(), currElemIds.size() * numValsPerComponent);
293
294
					for(std::vector<unsigned int>::const_iterator iter = currElemIds.begin(); iter != currElemIds.end(); ++iter)
					{
295
						createdAttribute->insert(index, currAttribute, *iter * numValsPerComponent, numValsPerComponent);
296
297
298
						index += numValsPerComponent;
					}
				}
299
				else if(currAttribute->getCenter() == XdmfAttributeCenter::Node())
300
301
302
				{
					createdAttribute = XdmfAttribute::New();
					createdAttribute->setName(currAttribute->getName());
303
304
					createdAttribute->setCenter(currAttribute->getCenter());
					createdAttribute->setType(currAttribute->getType());
305
					createdAttribute->initialize(currAttribute->getArrayType(), currNodeMap.size());
306
307
					for(std::map<unsigned int, unsigned int>::const_iterator iter = currNodeMap.begin(); iter != currNodeMap.end(); ++iter)
					{
308
						createdAttribute->insert(iter->second, currAttribute, iter->first, 1);
309
310
					}
				}
311
				if(createdAttribute)
312
313
314
315
				{
					partitioned->insert(createdAttribute);
					if(heavyDataWriter)
					{
316
317
318
319
						if(!createdAttribute->isInitialized())
						{
							createdAttribute->read();
						}
320
321
						createdAttribute->accept(heavyDataWriter);
						createdAttribute->release();
322
323
324
325
					}
				}
			}
		}
326
327
328
329
		if(releaseAttribute)
		{
			currAttribute->release();
		}
330
	}
331

332
	// Split sets into proper partitions
333
	for(unsigned int i=0; i<gridToPartition->getNumberSets(); ++i)
334
335
	{
		boost::shared_ptr<XdmfSet> currSet = gridToPartition->getSet(i);
336
		bool releaseSet = false;
337
		if(!currSet->isInitialized())
338
		{
339
			currSet->read();
340
			releaseSet = true;
341
342
343
344
345
346
347
348
		}
		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)
			{
349
				boost::shared_ptr<XdmfGridUnstructured> partitioned = partitionedGrids->getGridUnstructured(partitionId);
350
				partitionId++;
351
352
353
354

				boost::shared_ptr<XdmfSet> partitionedSet = XdmfSet::New();

				if(currSet->getType() == XdmfSetType::Cell() || currSet->getType() == XdmfSetType::Face() || currSet->getType() == XdmfSetType::Edge())
355
				{
356
					for(unsigned int k=0; k<currSet->getSize(); ++k)
357
					{
358
						std::vector<unsigned int>::const_iterator val = std::find(currElemIds.begin(), currElemIds.end(), currSet->getValue<unsigned int>(k));
359
360
361
						if(val != currElemIds.end())
						{
							unsigned int valToPush = val - currElemIds.begin();
362
							partitionedSet->pushBack(valToPush);
363
364
365
						}
					}
				}
366
				else if(currSet->getType() == XdmfSetType::Node())
367
				{
368
					for(unsigned int k=0; k<currSet->getSize(); ++k)
369
					{
370
						std::map<unsigned int, unsigned int>::const_iterator val = currNodeMap.find(currSet->getValue<unsigned int>(k));
371
372
						if(val != currNodeMap.end())
						{
373
							partitionedSet->pushBack(val->second);
374
375
376
						}
					}
				}
377
				if(partitionedSet->getSize() > 0)
378
				{
379
					partitioned->insert(partitionedSet);
380
381
					partitionedSet->setName(currSet->getName());
					partitionedSet->setType(currSet->getType());
382
383
					if(heavyDataWriter)
					{
384
385
						partitionedSet->accept(heavyDataWriter);
						partitionedSet->release();
386
387
388
389
					}
				}
			}
		}
390
391
392
393
		if(releaseSet)
		{
			currSet->release();
		}
394
	}
395
396
397

	// Add XdmfMap to map boundary nodes between partitions
	std::vector<boost::shared_ptr<XdmfAttribute> > globalNodeIds;
398
	for(unsigned int i=0; i<partitionedGrids->getNumberGridUnstructureds(); ++i)
399
	{
400
		boost::shared_ptr<XdmfAttribute> globalNodeId = partitionedGrids->getGridUnstructured(i)->getAttribute("GlobalNodeId");
401
402
403
404
		if(!globalNodeId->isInitialized())
		{
			globalNodeId->read();
		}
405
406
407
408
		globalNodeIds.push_back(globalNodeId);
	}

	std::vector<boost::shared_ptr<XdmfMap> > maps = XdmfMap::New(globalNodeIds);
409
	for(unsigned int i=0; i<partitionedGrids->getNumberGridUnstructureds(); ++i)
410
411
	{
		boost::shared_ptr<XdmfMap> map = maps[i];
412
		partitionedGrids->getGridUnstructured(i)->setMap(map);
413
414
415
416
417
418
419
420
		if(heavyDataWriter)
		{
			globalNodeIds[i]->release();
			map->accept(heavyDataWriter);
			map->release();
		}
	}

421
	return partitionedGrids;
422
423
}

424
425
#else

426
#include <iostream>
427
428
429
430
#include <sstream>
#include "XdmfDomain.hpp"
#include "XdmfGridCollection.hpp"
#include "XdmfHDF5Writer.hpp"
431
#include "XdmfPartitioner.hpp"
432
433
#include "XdmfReader.hpp"
#include "XdmfWriter.hpp"
434
435
436
437

/**
 * XdmfPartitioner is a command line utility for partitioning Xdmf grids.
 * The XdmfPartitioner uses the metis library to partition Triangular, Quadrilateral, Tetrahedral,
438
 * and Hexahedral XdmfGridUnstructureds.
439
440
441
442
443
 *
 * Usage:
 *     XdmfPartitioner <path-of-file-to-partition> <num-partitions> (Optional: <path-to-output-file>)
 *
 */
444
int main(int argc, char* argv[])
445
{
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
	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";
	}

	boost::shared_ptr<XdmfReader> reader = XdmfReader::New();
	boost::shared_ptr<XdmfDomain> domain = boost::shared_dynamic_cast<XdmfDomain>(reader->read(argv[1]));

493
	if(domain->getNumberGridUnstructureds() <= 0)
494
495
496
497
498
499
500
501
502
503
	{
		std::cout << "No grids to partition" << std::endl;
		return 1;
	}

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

	boost::shared_ptr<XdmfPartitioner> partitioner = XdmfPartitioner::New();
504
	boost::shared_ptr<XdmfGridCollection> toWrite = partitioner->partition(domain->getGridUnstructured(0), numPartitions, heavyDataWriter);
505
506
507
508
509
510
511
512
513

	boost::shared_ptr<XdmfDomain> newDomain = XdmfDomain::New();
	newDomain->insert(toWrite);

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

514
	std::cout << "Wrote: " << xmlFileName.str() << std::endl;
515
}
516
517

#endif //BUILD_EXE