← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/mesh/MeshPartition.cpp

  • Committer: Niclas Jansson
  • Date: 2011-06-10 14:33:43 UTC
  • Revision ID: njansson@csc.kth.se-20110610143343-d21p4am8rghiojfm
Added rudimentary header to binary files

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mesh/MeshPartition.cpp
2
2
// Licensed under the GNU LGPL Version 2.1.
3
3
//
4
4
// Modified by Anders Logg, 2008.
 
5
// Modified by Niclas Jansson, 2008-2009.
5
6
//
6
7
// First added:  2007-04-03
7
 
// Last changed: 2008-02-11
 
8
// Last changed: 2009-06-20
8
9
 
 
10
#include <dolfin/config/dolfin_config.h>
9
11
#include <dolfin/graph/Graph.h>
10
12
#include <dolfin/graph/GraphPartition.h>
11
 
#include "MeshPartition.h"
12
 
#include "MeshFunction.h"
 
13
#include <dolfin/mesh/MeshPartition.h>
 
14
#include <dolfin/mesh/MeshFunction.h>
 
15
#include <dolfin/mesh/MeshRenumber.h>
13
16
#include <dolfin/parameter/parameters.h>
14
17
 
 
18
#include <dolfin/mesh/Vertex.h>
 
19
#include <dolfin/mesh/Cell.h>
 
20
 
 
21
#ifdef HAVE_MPI
 
22
#include <mpi.h>
 
23
#include <parmetis.h>
 
24
#endif
 
25
 
 
26
 
 
27
 
 
28
 
15
29
using namespace dolfin;
16
30
 
17
31
//-----------------------------------------------------------------------------
28
42
    GraphPartition::edgecut(graph, num_partitions, partitions.values());
29
43
}
30
44
//-----------------------------------------------------------------------------
 
45
void MeshPartition::partition(Mesh& mesh, MeshFunction<uint>& partitions)
 
46
{
 
47
  partitionCommonMetis(mesh, partitions, 0);
 
48
}
 
49
//-----------------------------------------------------------------------------
 
50
void MeshPartition::partition(Mesh& mesh, MeshFunction<uint>& partitions,
 
51
                              MeshFunction<uint>& weight)
 
52
{
 
53
  partitionCommonMetis(mesh, partitions, &weight);
 
54
}
 
55
//-----------------------------------------------------------------------------
 
56
#ifdef HAVE_MPI
 
57
//-----------------------------------------------------------------------------
 
58
void MeshPartition::partitionCommonMetis(Mesh& mesh, 
 
59
                                         MeshFunction<uint>& partitions,
 
60
                                         MeshFunction<uint>* weight)
 
61
{
 
62
 
 
63
  // Metis assumes vertices numbered from process 0 
 
64
  MeshRenumber::renumber_vertices(mesh);
 
65
 
 
66
  float ubvec = 1.05; 
 
67
  int numflag = 0;    // C-style numbering
 
68
  int edgecut = 0;    
 
69
  int wgtflag, ncon;
 
70
  if( weight ) {
 
71
    wgtflag = 2;    // Weights on vertices only
 
72
    ncon = 1;       // One weight per vertex
 
73
  }
 
74
  else {
 
75
    wgtflag = 0;    // Turn off graph weights 
 
76
    ncon = 0;       // No weights on vertices
 
77
  }
 
78
 
 
79
  // Duplicate MPI communicator
 
80
  MPI_Comm comm;
 
81
  MPI_Comm_dup(MPI::DOLFIN_COMM, &comm);
 
82
 
 
83
  // Get information about the PE
 
84
  int size, rank;  
 
85
  MPI_Comm_size(MPI::DOLFIN_COMM, &size);
 
86
  MPI_Comm_rank(MPI::DOLFIN_COMM, &rank);
 
87
 
 
88
 
 
89
  idxtype *elmdist = new idxtype[size + 1];
 
90
  int ncells = mesh.numCells();
 
91
  elmdist[rank] = ncells;
 
92
  MPI_Allgather(&elmdist[rank], 1, MPI_INT, elmdist, 
 
93
                1, MPI_INT, MPI::DOLFIN_COMM);
 
94
 
 
95
  idxtype *elmwgt = NULL;
 
96
  if( weight ) {
 
97
    elmwgt = new idxtype[ncells];
 
98
    for(CellIterator c(mesh); !c.end(); ++c) 
 
99
      elmwgt[c->index()] = static_cast<idxtype>(weight->get(*c));
 
100
  }
 
101
 
 
102
  int sum_elm = elmdist[0];  
 
103
  int tmp_elm;
 
104
  elmdist[0] = 0;
 
105
  for(int i=1;i<size+1;i++){    
 
106
    tmp_elm = elmdist[i];
 
107
    elmdist[i] = sum_elm;
 
108
    sum_elm = tmp_elm + sum_elm;
 
109
  }
 
110
 
 
111
  int nvertices = mesh.type().numVertices(mesh.topology().dim());
 
112
  int ncnodes = nvertices - 1 ;
 
113
 
 
114
  idxtype *eptr = new idxtype[ncells + 1];
 
115
  eptr[0] = 0;
 
116
  for(uint i=1;i < (mesh.numCells() + 1);i++)
 
117
    eptr[i] = eptr[i-1] + nvertices;
 
118
 
 
119
  int *eind =  new idxtype[nvertices * ncells];  
 
120
  int i = 0;
 
121
  for(CellIterator c(mesh); !c.end(); ++c)
 
122
    for(VertexIterator v(*c); !v.end(); ++v)
 
123
      eind[i++] = mesh.distdata().get_global(*v);
 
124
 
 
125
  idxtype *part = new idxtype[ncells];
 
126
 
 
127
  float *tpwgts = new float[size];
 
128
  for(i=0; i<size; i++)
 
129
    tpwgts[i] =  1.0/(float)(size); 
 
130
 
 
131
  // default options
 
132
  int options[3] = {1, 0, 15};
 
133
 
 
134
 
 
135
  ParMETIS_V3_PartMeshKway(elmdist, eptr, eind, elmwgt, &wgtflag,&numflag,
 
136
                           &ncon,&ncnodes,&size, tpwgts, &ubvec,
 
137
                           options, &edgecut, part,&comm);
 
138
 
 
139
  delete[] eind;
 
140
  delete[] elmdist;
 
141
  delete[] tpwgts;
 
142
  delete[] eptr;
 
143
  if(weight)
 
144
    delete[] elmwgt;
 
145
 
 
146
  // Create partition function
 
147
  partitions.init(mesh, mesh.topology().dim());
 
148
  partitions = size;
 
149
  for(CellIterator cell(mesh); !cell.end(); ++cell)
 
150
    partitions.set(*cell, (uint) part[ cell->index() ]);
 
151
 
 
152
  delete[] part;
 
153
  MPI_Comm_free(&comm);
 
154
}
 
155
//-----------------------------------------------------------------------------
 
156
void MeshPartition::partition_geom(Mesh& mesh, MeshFunction<uint>& partitions)
 
157
{
 
158
  // Duplicate MPI communicator
 
159
  MPI_Comm comm; 
 
160
  MPI_Comm_dup(MPI::DOLFIN_COMM, &comm);
 
161
 
 
162
  int size, rank;
 
163
  // Get information about the PE
 
164
  MPI_Comm_size(MPI::DOLFIN_COMM, &size);
 
165
  MPI_Comm_rank(MPI::DOLFIN_COMM, &rank);
 
166
  
 
167
  // Gather number of locally stored vertices for each processor
 
168
  idxtype *vtxdist = new idxtype[size+1];  
 
169
  vtxdist[rank] = static_cast<idxtype> (mesh.numVertices());
 
170
 
 
171
 
 
172
  MPI_Allgather(&vtxdist[rank], 1, MPI_INT, vtxdist, 1, 
 
173
                MPI_INT, MPI::DOLFIN_COMM);
 
174
 
 
175
  int i,tmp;
 
176
  int sum = vtxdist[0];  
 
177
  vtxdist[0] = 0;
 
178
  for(i=1;i<size+1;i++){    
 
179
    tmp = vtxdist[i];
 
180
    vtxdist[i] = sum;
 
181
    sum = tmp + sum;
 
182
  }
 
183
 
 
184
  idxtype *part = new idxtype[mesh.numVertices()];
 
185
  int gdim =  static_cast<int>( mesh.geometry().dim() );
 
186
  float *xdy = new float[gdim * mesh.numVertices()];
 
187
 
 
188
  i = 0;
 
189
  for(VertexIterator vertex(mesh); !vertex.end(); ++vertex) {
 
190
    xdy[i] = static_cast<float>(vertex->point().x());
 
191
    xdy[i+1] = static_cast<float>(vertex->point().y());
 
192
    if(gdim > 2)
 
193
      xdy[i+2] = static_cast<float>(vertex->point().z());
 
194
    i +=gdim;
 
195
  }
 
196
 
 
197
  ParMETIS_V3_PartGeom(vtxdist,&gdim,xdy,part,&comm);
 
198
 
 
199
  // Create meshfunction from partitions
 
200
  partitions.init(mesh,0);
 
201
  for(VertexIterator vertex(mesh); !vertex.end(); ++vertex)
 
202
    partitions.set(*vertex, static_cast<uint>( part[vertex->index()]) );
 
203
 
 
204
  delete[] xdy;
 
205
  delete[] part;
 
206
  delete[] vtxdist;
 
207
  MPI_Comm_free(&comm);
 
208
}
 
209
//-----------------------------------------------------------------------------
 
210
#else
 
211
//-----------------------------------------------------------------------------
 
212
void MeshPartition::partitionCommonMetis(Mesh& mesh, 
 
213
                                         MeshFunction<uint>& partitions,
 
214
                                         MeshFunction<uint>* weight)
 
215
{
 
216
  error("ParMetis needs MPI");
 
217
}
 
218
//-----------------------------------------------------------------------------
 
219
void MeshPartition::partition_geom(Mesh& mesh, MeshFunction<uint>& partitions)
 
220
{
 
221
  error("ParMetis needs MPI");
 
222
}
 
223
//-----------------------------------------------------------------------------
 
224
#endif
 
225