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