1 /*******************************************************************/
2 /* XDMF */
3 /* eXtensible Data Model and Format */
4 /* */
5 /* Id : Id */
6 /* Date : $Date$ */
7 /* Version : $Revision$ */
8 /* */
9 /* Author: */
10 /* Kenneth Leiter */
11 /* kenneth.leiter@arl.army.mil */
12 /* US Army Research Laboratory */
13 /* Aberdeen Proving Ground, MD */
14 /* */
15 /* Copyright @ 2010 US Army Research Laboratory */
16 /* All Rights Reserved */
17 /* See Copyright.txt or http://www.arl.hpc.mil/ice for details */
18 /* */
19 /* This software is distributed WITHOUT ANY WARRANTY; without */
20 /* even the implied warranty of MERCHANTABILITY or FITNESS */
21 /* FOR A PARTICULAR PURPOSE. See the above copyright notice */
22 /* for more information. */
23 /* */
24 /*******************************************************************/
25
26 #include "XdmfPartitioner.h"
27
28 #include <sstream>
29
30 #ifndef BUILD_EXE
31
32 extern "C"
33 {
34 #include <metis.h>
35 }
36
37 #include <map>
38 #include <vector>
39
XdmfPartitioner()40 XdmfPartitioner::XdmfPartitioner()
41 {
42 }
43
~XdmfPartitioner()44 XdmfPartitioner::~XdmfPartitioner()
45 {
46 }
47
Partition(XdmfGrid * grid,int numPartitions,XdmfElement * parentElement)48 XdmfGrid * XdmfPartitioner::Partition(XdmfGrid * grid, int numPartitions, XdmfElement * parentElement)
49 {
50 int metisElementType;
51 int numElementsPerNode;
52 switch(grid->GetTopology()->GetTopologyType())
53 {
54 case(XDMF_TRI):
55 case(XDMF_TRI_6):
56 metisElementType = 1;
57 numElementsPerNode = 3;
58 break;
59 case(XDMF_QUAD):
60 case(XDMF_QUAD_8):
61 case(XDMF_QUAD_9):
62 metisElementType = 4;
63 numElementsPerNode = 4;
64 break;
65 case(XDMF_TET):
66 case(XDMF_TET_10):
67 metisElementType = 2;
68 numElementsPerNode = 4;
69 break;
70 case(XDMF_HEX):
71 case(XDMF_HEX_20):
72 case(XDMF_HEX_24):
73 case(XDMF_HEX_27):
74 metisElementType = 3;
75 numElementsPerNode = 8;
76 break;
77 default:
78 std::cout << "Cannot partition grid with element type: " << grid->GetTopology()->GetTopologyTypeAsString() << std::endl;
79 return NULL;
80 }
81
82 int numNodes = grid->GetGeometry()->GetNumberOfPoints() / (grid->GetTopology()->GetNodesPerElement() / numElementsPerNode);
83 int numElements = grid->GetTopology()->GetNumberOfElements();
84
85 idxtype * metisConnectivity = new idxtype[numElementsPerNode * numElements];
86 for(int i=0; i<numElements; ++i)
87 {
88 grid->GetTopology()->GetConnectivity()->GetValues(i*grid->GetTopology()->GetNodesPerElement(), &metisConnectivity[i*numElementsPerNode], numElementsPerNode);
89 }
90
91 // Need to remap connectivity for quadratic elements so that metis handles it properly
92 std::map<idxtype, idxtype> xdmfConnIdxToMetisConnIdx;
93 if(numNodes != grid->GetGeometry()->GetNumberOfPoints())
94 {
95 int index = 0;
96 for(int i=0; i<numElements * numElementsPerNode; ++i)
97 {
98 std::map<idxtype, idxtype>::const_iterator val = xdmfConnIdxToMetisConnIdx.find(metisConnectivity[i]);
99 if(val != xdmfConnIdxToMetisConnIdx.end())
100 {
101 metisConnectivity[i] = val->second;
102 }
103 else
104 {
105 // Haven't seen this id before, map to index and set to new id
106 xdmfConnIdxToMetisConnIdx[metisConnectivity[i]] = index;
107 metisConnectivity[i] = index;
108 index++;
109 }
110 }
111 }
112
113 int startIndex = 0;
114 int numCutEdges = 0;
115
116 idxtype * elementsPartition = new idxtype[numElements];
117 idxtype * nodesPartition = new idxtype[numNodes];
118
119 //METIS_PartMeshDual(&numElements, &numNodes, metisConnectivity, &metisElementType, &startIndex, &numPartitions, &numCutEdges, elementsPartition, nodesPartition);
120 std::cout << "Entered METIS" << std::endl;
121 METIS_PartMeshNodal(&numElements, &numNodes, metisConnectivity, &metisElementType, &startIndex, &numPartitions, &numCutEdges, elementsPartition, nodesPartition);
122 std::cout << "Exited METIS" << std::endl;
123
124 delete [] metisConnectivity;
125 delete [] nodesPartition;
126
127 // Initialize sets to hold globalNodeIds for each partition.
128 std::vector<std::map<XdmfInt32, XdmfInt32> > globalToLocalNodeIdMap;
129 std::vector<std::map<XdmfInt32, XdmfInt32> > globalToLocalElementIdMap;
130 for(int i=0; i<numPartitions; ++i)
131 {
132 std::map<XdmfInt32, XdmfInt32> nodeMap;
133 globalToLocalNodeIdMap.push_back(nodeMap);
134 std::map<XdmfInt32, XdmfInt32> elemMap;
135 globalToLocalElementIdMap.push_back(elemMap);
136 }
137
138 // Fill in globalNodeId for each partition
139 XdmfInt32 * conn = new XdmfInt32[numElements * grid->GetTopology()->GetNodesPerElement()];
140 grid->GetTopology()->GetConnectivity()->GetValues(0, conn, numElements * grid->GetTopology()->GetNodesPerElement());
141 int totalIndex = 0;
142 for(int i=0; i<numElements; ++i)
143 {
144 for(int j=0; j<grid->GetTopology()->GetNodesPerElement(); ++j)
145 {
146 if(globalToLocalNodeIdMap[elementsPartition[i]].count(conn[totalIndex]) == 0)
147 {
148 // Have not seen this node, need to add to map
149 globalToLocalNodeIdMap[elementsPartition[i]][conn[totalIndex]] = globalToLocalNodeIdMap[elementsPartition[i]].size() - 1;
150 }
151 totalIndex++;
152 }
153 globalToLocalElementIdMap[elementsPartition[i]][i] = globalToLocalElementIdMap[elementsPartition[i]].size() - 1;
154 }
155
156 delete [] elementsPartition;
157
158 XdmfGrid * collection = new XdmfGrid();
159 collection->SetName("Collection");
160 collection->SetGridType(XDMF_GRID_COLLECTION);
161 collection->SetCollectionType(XDMF_GRID_COLLECTION_SPATIAL);
162 collection->SetDeleteOnGridDelete(true);
163
164 parentElement->Insert(collection);
165
166 for(int i=0; i<numPartitions; ++i)
167 {
168 std::map<XdmfInt32, XdmfInt32> currNodeMap = globalToLocalNodeIdMap[i];
169 std::map<XdmfInt32, XdmfInt32> currElemMap = globalToLocalElementIdMap[i];
170
171 if(currNodeMap.size() > 0)
172 {
173 std::stringstream name;
174 name << grid->GetName() << "_" << i;
175
176 XdmfGrid * partition = new XdmfGrid();
177 partition->SetName(name.str().c_str());
178
179 int numDimensions = 3;
180 if(grid->GetGeometry()->GetGeometryType() == XDMF_GEOMETRY_XY || grid->GetGeometry()->GetGeometryType() == XDMF_GEOMETRY_X_Y)
181 {
182 numDimensions = 2;
183 }
184
185 XdmfGeometry * geom = partition->GetGeometry();
186 geom->SetGeometryType(grid->GetGeometry()->GetGeometryType());
187 geom->SetNumberOfPoints(currNodeMap.size());
188 geom->SetDeleteOnGridDelete(true);
189
190 XdmfArray * points = geom->GetPoints();
191 points->SetNumberType(grid->GetGeometry()->GetPoints()->GetNumberType());
192 points->SetNumberOfElements(currNodeMap.size() * numDimensions);
193 for(std::map<XdmfInt32, XdmfInt32>::const_iterator iter = currNodeMap.begin(); iter != currNodeMap.end(); ++iter)
194 {
195 points->SetValues(iter->second * numDimensions, grid->GetGeometry()->GetPoints(), numDimensions, iter->first * 3);
196 }
197
198 XdmfTopology * top = partition->GetTopology();
199 top->SetTopologyType(grid->GetTopology()->GetTopologyType());
200 top->SetNumberOfElements(currElemMap.size());
201 top->SetDeleteOnGridDelete(true);
202
203 XdmfArray * connections = top->GetConnectivity();
204 connections->SetNumberType(grid->GetTopology()->GetConnectivity()->GetNumberType());
205 connections->SetNumberOfElements(currElemMap.size() * grid->GetTopology()->GetNodesPerElement());
206 XdmfInt32 * tmpConn = new XdmfInt32[grid->GetTopology()->GetNodesPerElement()];
207 for(std::map<XdmfInt32, XdmfInt32>::const_iterator iter = currElemMap.begin(); iter != currElemMap.end(); ++iter)
208 {
209 // Get global connectivity values for this element
210 grid->GetTopology()->GetConnectivity()->GetValues(iter->first * grid->GetTopology()->GetNodesPerElement(), tmpConn, grid->GetTopology()->GetNodesPerElement());
211 // Translate these global points to local node ids
212 for(int j=0; j<grid->GetTopology()->GetNodesPerElement(); ++j)
213 {
214 tmpConn[j] = currNodeMap[tmpConn[j]];
215 }
216 connections->SetValues(iter->second * grid->GetTopology()->GetNodesPerElement(), tmpConn, grid->GetTopology()->GetNodesPerElement());
217 }
218 delete [] tmpConn;
219 collection->Insert(partition);
220
221 // Add GlobalNodeId Attribute
222 XdmfAttribute * globalIds = new XdmfAttribute();
223 globalIds->SetName("GlobalNodeId");
224 globalIds->SetAttributeType(XDMF_ATTRIBUTE_TYPE_SCALAR);
225 globalIds->SetAttributeCenter(XDMF_ATTRIBUTE_CENTER_NODE);
226 globalIds->SetDeleteOnGridDelete(true);
227
228 XdmfArray * globalNodeIdVals = globalIds->GetValues();
229 globalNodeIdVals->SetNumberType(XDMF_INT32_TYPE);
230 globalNodeIdVals->SetNumberOfElements(currNodeMap.size());
231 for(std::map<XdmfInt32, XdmfInt32>::const_iterator iter = currNodeMap.begin(); iter != currNodeMap.end(); ++iter)
232 {
233 globalNodeIdVals->SetValues(iter->second, (XdmfInt32*)&iter->first, 1);
234 }
235 partition->Insert(globalIds);
236
237 // Split attributes and add to grid
238 for(int j=0; j<grid->GetNumberOfAttributes(); ++j)
239 {
240 XdmfAttribute * currAttribute = grid->GetAttribute(j);
241 currAttribute->Update();
242 switch(currAttribute->GetAttributeCenter())
243 {
244 case(XDMF_ATTRIBUTE_CENTER_GRID):
245 {
246 // Will continue to be true for entire collection - so insert at top level
247 collection->Insert(currAttribute);
248 break;
249 }
250 case(XDMF_ATTRIBUTE_CENTER_CELL):
251 case(XDMF_ATTRIBUTE_CENTER_FACE):
252 case(XDMF_ATTRIBUTE_CENTER_EDGE):
253 {
254 XdmfAttribute * attribute = new XdmfAttribute();
255 attribute->SetName(currAttribute->GetName());
256 attribute->SetAttributeType(currAttribute->GetAttributeType());
257 attribute->SetAttributeCenter(currAttribute->GetAttributeCenter());
258 attribute->SetDeleteOnGridDelete(true);
259 XdmfArray * attributeVals = attribute->GetValues();
260 attributeVals->SetNumberType(currAttribute->GetValues()->GetNumberType());
261 int numValsPerComponent = currAttribute->GetValues()->GetNumberOfElements() / grid->GetTopology()->GetNumberOfElements();
262 attributeVals->SetNumberOfElements(currElemMap.size() * numValsPerComponent);
263 for(std::map<XdmfInt32, XdmfInt32>::const_iterator iter = currElemMap.begin(); iter != currElemMap.end(); ++iter)
264 {
265 attributeVals->SetValues(iter->second * numValsPerComponent, currAttribute->GetValues(), numValsPerComponent, iter->first * numValsPerComponent);
266 }
267 partition->Insert(attribute);
268 break;
269 }
270 case(XDMF_ATTRIBUTE_CENTER_NODE):
271 {
272 XdmfAttribute * attribute = new XdmfAttribute();
273 attribute->SetName(currAttribute->GetName());
274 attribute->SetAttributeType(currAttribute->GetAttributeType());
275 attribute->SetAttributeCenter(currAttribute->GetAttributeCenter());
276 attribute->SetDeleteOnGridDelete(true);
277 XdmfArray * attributeVals = attribute->GetValues();
278 attributeVals->SetNumberType(currAttribute->GetValues()->GetNumberType());
279 attributeVals->SetNumberOfElements(currNodeMap.size());
280 for(std::map<XdmfInt32, XdmfInt32>::const_iterator iter = currNodeMap.begin(); iter != currNodeMap.end(); ++iter)
281 {
282 attributeVals->SetValues(iter->second, currAttribute->GetValues(), 1, iter->first);
283 }
284 partition->Insert(attribute);
285 break;
286 }
287 default:
288 {
289 std::cout << "Unknown attribute center encountered: " << currAttribute->GetAttributeCenterAsString() << std::endl;
290 break;
291 }
292 }
293 }
294
295 // Split sets and add to grid
296 for(int j=0; j<grid->GetNumberOfSets(); ++j)
297 {
298 XdmfSet * currSet = grid->GetSets(j);
299 currSet->Update();
300 switch(currSet->GetSetType())
301 {
302 case(XDMF_SET_TYPE_CELL):
303 case(XDMF_SET_TYPE_FACE):
304 case(XDMF_SET_TYPE_EDGE):
305 {
306 std::vector<XdmfInt32> myIds;
307 for(int j=0; j<currSet->GetIds()->GetNumberOfElements(); j++)
308 {
309 std::map<XdmfInt32, XdmfInt32>::const_iterator val = currElemMap.find(currSet->GetIds()->GetValueAsInt32(j));
310 if(val != currNodeMap.end())
311 {
312 myIds.push_back(val->second);
313 }
314 }
315 if(myIds.size() != 0)
316 {
317 XdmfSet * set = new XdmfSet();
318 set->SetName(currSet->GetName());
319 set->SetSetType(currSet->GetSetType());
320 set->SetSize(myIds.size());
321 set->SetDeleteOnGridDelete(true);
322 XdmfArray * ids = set->GetIds();
323 ids->SetNumberType(XDMF_INT32_TYPE);
324 ids->SetNumberOfElements(myIds.size());
325 ids->SetValues(0, &myIds[0], myIds.size());
326 partition->Insert(set);
327 }
328 break;
329 }
330 case(XDMF_SET_TYPE_NODE):
331 {
332 std::vector<XdmfInt32> myIds;
333 for(int j=0; j<currSet->GetIds()->GetNumberOfElements(); j++)
334 {
335 std::map<XdmfInt32, XdmfInt32>::const_iterator val = currNodeMap.find(currSet->GetIds()->GetValueAsInt32(j));
336 if(val != currNodeMap.end())
337 {
338 myIds.push_back(val->second);
339 }
340 }
341 if(myIds.size() != 0)
342 {
343 XdmfSet * set = new XdmfSet();
344 set->SetName(currSet->GetName());
345 set->SetSetType(currSet->GetSetType());
346 set->SetSize(myIds.size());
347 set->SetDeleteOnGridDelete(true);
348 XdmfArray * ids = set->GetIds();
349 ids->SetNumberType(XDMF_INT32_TYPE);
350 ids->SetNumberOfElements(myIds.size());
351 ids->SetValues(0, &myIds[0], myIds.size());
352 partition->Insert(set);
353 }
354 break;
355 }
356 default:
357 {
358 std::cout << "Unknown set type encountered: " << currSet->GetSetTypeAsString() << std::endl;
359 break;
360 }
361 }
362 }
363 std::cout << "HERE" << std::endl;
364 }
365 }
366
367 return collection;
368 }
369
370 #else
371
372 /**
373 * XdmfPartitioner is a command line utility for partitioning Xdmf grids.
374 * The XdmfPartitioner uses the metis library to partition Triangular, Quadrilateral, Tetrahedral,
375 * and Hexahedral XdmfGrids.
376 *
377 * Usage:
378 * XdmfPartitioner <path-of-file-to-partition> <num-partitions> (Optional: <path-to-output-file>)
379 *
380 */
main(int argc,char * argv[])381 int main(int argc, char* argv[])
382 {
383 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>)";
384 std::string meshName = "";
385
386 if (argc < 3)
387 {
388 cout << usage << endl;
389 return 1;
390 }
391
392 FILE * refFile = fopen(argv[1], "r");
393 if (refFile)
394 {
395 // Success
396 meshName = argv[1];
397 fclose(refFile);
398 }
399 else
400 {
401 cout << "Cannot open file: " << argv[1] << endl;
402 return 1;
403 }
404
405 int numPartitions = atoi(argv[2]);
406
407 if (argc >= 4)
408 {
409 meshName = argv[3];
410 }
411
412 if(meshName.find_last_of("/\\") != std::string::npos)
413 {
414 meshName = meshName.substr(meshName.find_last_of("/\\")+1, meshName.length());
415 }
416
417 if (meshName.rfind(".") != std::string::npos)
418 {
419 meshName = meshName.substr(0, meshName.rfind("."));
420 }
421
422 if(argc < 4)
423 {
424 meshName = meshName + "-partitioned";
425 }
426
427 XdmfDOM dom;
428 XdmfInt32 error = dom.Parse(argv[1]);
429 if(error == XDMF_FAIL)
430 {
431 std::cout << "File does not appear to be a valid Xdmf file" << std::endl;
432 return 1;
433 }
434 XdmfXmlNode gridElement = dom.FindElementByPath("/Xdmf/Domain/Grid");
435 if(gridElement == NULL)
436 {
437 std::cout << "Cannot parse Xdmf file!" << std::endl;
438 return 1;
439 }
440
441 XdmfGrid * grid = new XdmfGrid();
442 grid->SetDOM(&dom);
443 grid->SetElement(gridElement);
444 grid->Update();
445
446 XdmfDOM newDOM;
447 XdmfRoot newRoot;
448 XdmfDomain newDomain;
449
450 newRoot.SetDOM(&newDOM);
451 newRoot.Build();
452 newRoot.Insert(&newDomain);
453
454 XdmfPartitioner partitioner;
455 XdmfGrid * partitionedGrid = partitioner.Partition(grid, numPartitions, &newDomain);
456
457 for(int i=0; i<partitionedGrid->GetNumberOfChildren(); ++i)
458 {
459 XdmfGrid * child = partitionedGrid->GetChild(i);
460
461 // Set heavy data set names for geometry and topology
462 std::stringstream heavyPointName;
463 heavyPointName << meshName << ".h5:/" << child->GetName() << "/XYZ";
464 child->GetGeometry()->GetPoints()->SetHeavyDataSetName(heavyPointName.str().c_str());
465
466 std::stringstream heavyConnName;
467 heavyConnName << meshName << ".h5:/" << child->GetName() << "/Connections";
468 child->GetTopology()->GetConnectivity()->SetHeavyDataSetName(heavyConnName.str().c_str());
469
470 // Set heavy data set names for mesh attributes and sets
471 for(int i=0; i<child->GetNumberOfAttributes(); i++)
472 {
473 std::stringstream heavyAttrName;
474 heavyAttrName << meshName << ".h5:/" << child->GetName() << "/Attribute/" << child->GetAttribute(i)->GetAttributeCenterAsString() << "/" << child->GetAttribute(i)->GetName();
475 child->GetAttribute(i)->GetValues()->SetHeavyDataSetName(heavyAttrName.str().c_str());
476 }
477
478 for(int i=0; i<child->GetNumberOfSets(); i++)
479 {
480 std::stringstream heavySetName;
481 heavySetName << meshName << ".h5:/" << child->GetName() << "/Set/" << child->GetSets(i)->GetSetTypeAsString() << "/" << child->GetSets(i)->GetName();
482 child->GetSets(i)->GetIds()->SetHeavyDataSetName(heavySetName.str().c_str());
483 }
484 }
485
486 partitionedGrid->Build();
487
488 std::stringstream outputFileName;
489 outputFileName << meshName << ".xmf";
490
491 newDOM.Write(outputFileName.str().c_str());
492
493 std::cout << "Wrote: " << outputFileName.str().c_str() << std::endl;
494 }
495
496 #endif //BUILD_EXE
497