← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/fem/DofMap.cpp

  • Committer: Niclas Jansson
  • Date: 2010-10-17 10:21:52 UTC
  • Revision ID: njansson@csc.kth.se-20101017102152-e6u3c9uwxa4z19c8
Minor fixes in configure.ac

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/fem/DofMap.cpp
2
2
// Licensed under the GNU LGPL Version 2.1.
3
3
 
4
4
// Modified by Martin Alnes, 2008
 
5
// Modified by Niclas Jansson, 2009
5
6
 
6
7
// First added:  2007-03-01
7
 
// Last changed: 2008-04-10
 
8
// Last changed: 2009-11-01
8
9
 
 
10
#include <dolfin/config/dolfin_config.h>
9
11
#include <dolfin/common/types.h>
10
12
#include <dolfin/mesh/Cell.h>
 
13
#include <dolfin/mesh/Vertex.h>
11
14
#include "UFCCell.h"
12
15
#include "DofMap.h"
13
16
#include "SubSystem.h"
16
19
#include "UFC.h"
17
20
#include <dolfin/main/MPI.h>
18
21
 
 
22
#include <dolfin/mesh/BoundaryMesh.h>
 
23
#include <dolfin/mesh/Facet.h>
 
24
#include <dolfin/mesh/MeshData.h>
 
25
#include <dolfin/mesh/MeshFunction.h>
 
26
#include <dolfin/mesh/GlobalFacetMap.h>
 
27
#include <cstring>
 
28
#include <cstdlib>
 
29
 
 
30
#ifdef HAVE_MPI
 
31
#include <mpi.h>
 
32
#endif
19
33
using namespace dolfin;
20
34
 
21
35
//-----------------------------------------------------------------------------
22
 
DofMap::DofMap(ufc::dof_map& dof_map, Mesh& mesh) : dof_map(0), 
 
36
DofMap::DofMap(ufc::dof_map& dof_map, Mesh& mesh, bool dof_map_local) : dof_map(0), 
23
37
               ufc_dof_map(&dof_map), ufc_dof_map_local(false), 
24
38
               dolfin_mesh(mesh), num_cells(mesh.numCells()), 
25
 
               partitions(0)
 
39
               partitions(0), _type_(-1), _local_size(0), v_map(0)
26
40
{
27
 
  init();
 
41
 
 
42
 
 
43
  // Assume responsibilty for ufc_dof_map
 
44
  if(dof_map_local) 
 
45
    ufc_dof_map_local = ufc_dof_map;
 
46
  init();  
 
47
 
 
48
  if(dolfin::MPI::numProcesses() > 1)
 
49
    build();
28
50
}
29
51
//-----------------------------------------------------------------------------
30
 
DofMap::DofMap(ufc::dof_map& dof_map, Mesh& mesh, MeshFunction<uint>& partitions)
31
 
             : dof_map(0), ufc_dof_map(&dof_map), ufc_dof_map_local(false), 
32
 
               dolfin_mesh(mesh), num_cells(mesh.numCells()), 
33
 
               partitions(&partitions)
 
52
DofMap::DofMap(ufc::dof_map& dof_map, Mesh& mesh, MeshFunction<uint>& partitions,
 
53
               bool dof_map_local) : dof_map(0), ufc_dof_map(&dof_map), 
 
54
               ufc_dof_map_local(false), dolfin_mesh(mesh), num_cells(mesh.numCells()), 
 
55
               partitions(&partitions),  _type_(-1), _local_size(0), v_map(0)
34
56
{
 
57
  // Assume responsibilty for ufc_dof_map
 
58
  if(dof_map_local) 
 
59
    ufc_dof_map_local = ufc_dof_map;
35
60
  init();
 
61
 
 
62
  if(dolfin::MPI::numProcesses() > 1)
 
63
    build();
36
64
}
37
65
//-----------------------------------------------------------------------------
38
66
DofMap::DofMap(const std::string signature, Mesh& mesh) 
39
67
  : dof_map(0), ufc_dof_map(0), ufc_dof_map_local(false),
40
 
    dolfin_mesh(mesh), num_cells(mesh.numCells()), partitions(0)
 
68
    dolfin_mesh(mesh), num_cells(mesh.numCells()), partitions(0), 
 
69
    _type_(-1), _local_size(0), v_map(0)
41
70
{
42
71
  // Create ufc dof map from signature
43
72
  ufc_dof_map = ElementLibrary::create_dof_map(signature);
48
77
  ufc_dof_map_local = ufc_dof_map;
49
78
 
50
79
  init();
 
80
 
 
81
  if(dolfin::MPI::numProcesses() > 1)
 
82
    build();
51
83
}
52
84
//-----------------------------------------------------------------------------
53
85
DofMap::DofMap(const std::string signature, Mesh& mesh, 
54
86
               MeshFunction<uint>& partitions) 
55
87
  : dof_map(0), ufc_dof_map(0), 
56
88
    ufc_dof_map_local(false), dolfin_mesh(mesh), num_cells(mesh.numCells()),
57
 
    partitions(&partitions)
 
89
    partitions(&partitions), _type_(-1), _local_size(0), v_map(0)
58
90
{
59
91
  // Create ufc dof map from signature
60
92
  ufc_dof_map = ElementLibrary::create_dof_map(signature);
65
97
  ufc_dof_map_local = ufc_dof_map;
66
98
 
67
99
  init();
 
100
 
 
101
  if(dolfin::MPI::numProcesses() > 1)
 
102
    build();
68
103
}
69
104
//-----------------------------------------------------------------------------
70
105
DofMap::~DofMap()
71
106
{
72
107
  if (dof_map)
73
 
  {
74
 
    for (uint i = 0; i < num_cells; ++i)
75
 
      delete [] dof_map[i];
76
108
    delete [] dof_map;
77
 
  }
78
109
 
79
110
  if (ufc_dof_map_local)
80
111
    delete ufc_dof_map_local;
 
112
 
 
113
  if(v_map)
 
114
    delete[] v_map;
81
115
}
82
116
//-----------------------------------------------------------------------------
83
117
DofMap* DofMap::extractDofMap(const Array<uint>& sub_system, uint& offset) const
84
118
{
85
119
  // Check that dof map has not be re-ordered
86
 
  if (dof_map)
87
 
    error("Dof map has been re-ordered. Don't yet know how to extract sub dof maps.");
 
120
  //  if (dof_map)
 
121
  //    error("Dof map has been re-ordered. Don't yet know how to extract sub dof maps.");
88
122
 
89
123
  // Reset offset
90
124
  offset = 0;
95
129
  message(2, "Offset for sub system: %d", offset);
96
130
 
97
131
  if (partitions)
98
 
    return new DofMap(*sub_dof_map, dolfin_mesh, *partitions);
 
132
    return new DofMap(*sub_dof_map, dolfin_mesh, *partitions, true);
99
133
  else
100
 
    return new DofMap(*sub_dof_map, dolfin_mesh);
 
134
    return new DofMap(*sub_dof_map, dolfin_mesh, true);
101
135
}
102
136
//-----------------------------------------------------------------------------
103
137
ufc::dof_map* DofMap::extractDofMap(const ufc::dof_map& dof_map, uint& offset, const Array<uint>& sub_system) const
168
202
    UFCCell ufc_cell(*cell);
169
203
    for (; !cell.end(); ++cell)
170
204
    {
171
 
      ufc_cell.update(*cell);
 
205
      ufc_cell.update(*cell, dolfin_mesh.distdata());
172
206
      ufc_dof_map->init_cell(ufc_mesh, ufc_cell);
173
207
    }
174
208
    ufc_dof_map->init_cell_finalize();
182
216
  // Either lookup pretabulated values (if build() has been called)
183
217
  // or ask the ufc::dof_map to tabulate the values
184
218
 
185
 
  if (dof_map)
186
 
  {
187
 
    for (uint i = 0; i < local_dimension(); i++)
188
 
      dofs[i] = dof_map[cell_index][i];
189
 
    //memcpy(dofs, dof_map[cell_index], sizeof(uint)*local_dimension()); // FIXME: Maybe memcpy() can be used to speed this up? Test this!
190
 
  }
191
 
  else
192
 
    ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);
193
 
}
194
 
//-----------------------------------------------------------------------------
195
 
void DofMap::build(UFC& ufc)
196
 
{
197
 
  printf("BUILDING\n");
198
 
 
199
 
  // Initialize new dof map
200
 
  if (dof_map)
201
 
    delete dof_map;
202
 
  dof_map = new uint*[dolfin_mesh.numCells()];
203
 
  
204
 
  // for all processes
205
 
  uint current_dof = 0;
206
 
  for (uint p = 0; p < MPI::numProcesses(); ++p)
207
 
  {
208
 
    // for all cells
209
 
    for (CellIterator c(dolfin_mesh); !c.end(); ++c)
210
 
    {
211
 
      // if cell in partition belonging to process p
212
 
      if ((*partitions)(*c) != p)
213
 
        continue;
214
 
 
215
 
      dof_map[c->index()] = new uint[local_dimension()];
216
 
      //dolfin_debug2("cpu %d building cell %d", MPI::processNumber(), c->index());
217
 
      ufc.update(*c);
218
 
      ufc_dof_map->tabulate_dofs(ufc.dofs[0], ufc.mesh, ufc.cell);
219
 
 
220
 
      for (uint i=0; i < ufc_dof_map->local_dimension(); ++i)
221
 
      {
222
 
        const uint dof = ufc.dofs[0][i];
223
 
        //dolfin_debug3("ufc.dofs[%d][%d] = %d", 0, MPI::processNumber(), ufc.dofs[0][i]);
224
 
 
225
 
        std::map<uint, uint>::iterator it = map.find(dof);
226
 
        if (it != map.end())
227
 
        {
228
 
          //dolfin_debug2("cpu %d dof %d already computed", MPI::processNumber(), it->second);
229
 
          dof_map[c->index()][i] = it->second;
230
 
        }
231
 
        else
232
 
        {
233
 
          dof_map[c->index()][i] = current_dof;
234
 
          map[dof] = current_dof++;
235
 
        }
236
 
      }
237
 
    }  
238
 
  }
239
 
}
240
 
//-----------------------------------------------------------------------------
241
 
std::map<dolfin::uint, dolfin::uint> DofMap::getMap() const
 
219
  if(_type_ == 0)
 
220
  {
 
221
    Cell c(dolfin_mesh, cell_index);
 
222
    for (uint i = 0; i < local_dimension(); i++)
 
223
      dofs[i] = dolfin_mesh.distdata().get_global(c.entities(0)[i], 0);
 
224
  }
 
225
  else if(_type_ == 1)
 
226
  {
 
227
    *dofs = dolfin_mesh.distdata().get_cell_global(cell_index);
 
228
  }
 
229
  else if(_type_ == 2 && v_map)
 
230
  {
 
231
    Cell c(dolfin_mesh, cell_index);
 
232
    uint gdim = ufc_dof_map->geometric_dimension();
 
233
    uint num_entities = c.numEntities(0);
 
234
    dolfin_assert(gdim * num_entities == local_dimension());
 
235
    for (uint k = 0; k < gdim; k++)
 
236
      for (uint i = 0;  i < num_entities; i++)
 
237
        dofs[i + k * (gdim + 1)] = v_map[c.entities(0)[i]] + k;
 
238
  }
 
239
  else if(_type_ == 3)
 
240
  {
 
241
    uint gdim = ufc_dof_map->geometric_dimension();
 
242
    dolfin_assert(local_dimension() == gdim);
 
243
    for(uint i = 0; i < gdim; i++)
 
244
      dofs[i] = (cell_index + i * dolfin_mesh.numCells()) + _offset_;
 
245
  }
 
246
  else if (dof_map)
 
247
  {
 
248
    uint offset = local_dimension() * cell_index;
 
249
    memcpy(dofs, &dof_map[offset], sizeof(uint)*local_dimension());
 
250
  }
 
251
  else
 
252
    ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);
 
253
} 
 
254
//-----------------------------------------------------------------------------
 
255
void DofMap::tabulate_dofs(uint* dofs, const ufc::cell& ufc_cell, uint cell_index) const
 
256
{
 
257
  // Either lookup pretabulated values (if build() has been called)
 
258
  // or ask the ufc::dof_map to tabulate the values
 
259
  if(_type_ == 0)
 
260
  {
 
261
    Cell c(dolfin_mesh, cell_index);
 
262
    for (uint i = 0; i < local_dimension(); i++)
 
263
      dofs[i] = dolfin_mesh.distdata().get_global(c.entities(0)[i], 0);
 
264
  }
 
265
  else if(_type_ == 1)
 
266
  {
 
267
    *dofs = dolfin_mesh.distdata().get_cell_global(cell_index);
 
268
  }
 
269
  else if(_type_ == 2 && v_map)
 
270
  {
 
271
    Cell c(dolfin_mesh, cell_index);
 
272
    uint gdim = ufc_dof_map->geometric_dimension();
 
273
    uint num_entities = c.numEntities(0);
 
274
    for (uint k = 0; k < gdim; k++)
 
275
      for (uint i = 0;  i < num_entities; i++)
 
276
        dofs[i + k * (gdim + 1)] = v_map[c.entities(0)[i]] + k;
 
277
  }
 
278
  else if(_type_ == 3)
 
279
  {
 
280
    uint gdim = ufc_dof_map->geometric_dimension();
 
281
    for(uint i = 0; i < gdim; i++)
 
282
      dofs[i] = (cell_index + i * dolfin_mesh.numCells()) + _offset_;
 
283
 
 
284
  }
 
285
  else if (dof_map)
 
286
  {
 
287
    uint offset = local_dimension() * cell_index;
 
288
    memcpy(dofs, &dof_map[offset], sizeof(uint)*local_dimension());
 
289
  }
 
290
  else
 
291
    ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);
 
292
 
 
293
}//-----------------------------------------------------------------------------
 
294
void DofMap::build()
 
295
{
 
296
 
 
297
  if( dof_map ) 
 
298
    delete [] dof_map;
 
299
 
 
300
  map.clear();
 
301
 
 
302
 
 
303
  if(MPI::numProcesses() == 1) {
 
304
    /*
 
305
     uint *dofs =  new uint[local_dimension()];
 
306
     
 
307
     dof_map = new uint*[dolfin_mesh.numCells()];
 
308
 
 
309
     uint num = 0;
 
310
 
 
311
     uint tt = local_dimension() / ufc_dof_map->geometric_dimension();
 
312
     for(CellIterator c(dolfin_mesh); !c.end(); ++c) {
 
313
       
 
314
       dof_map[c->index()] = new uint[local_dimension()];    
 
315
       ufc.update(*c, dolfin_mesh.distdata());
 
316
            
 
317
       for(uint j = 0; j < ufc.form.rank(); j++) {
 
318
         ufc_dof_map->tabulate_dofs(dofs, ufc.mesh, ufc.cell);      
 
319
 
 
320
         if(tt == 0) {
 
321
           for(uint i = 0; i < local_dimension(); i++) {
 
322
             const uint dof = dofs[i];
 
323
 
 
324
             std::map<uint, uint>::iterator it = map.find(dof);
 
325
             if(it != map.end())
 
326
               dof_map[c->index()][i] = it->second;
 
327
             else {
 
328
               dof_map[c->index()][i] = num;
 
329
               map[dof] = num++;
 
330
#ifdef BLOCKED
 
331
               dof_map[c->index()][i] /= 2;
 
332
               map[dof] /= 2;
 
333
#endif
 
334
             }
 
335
           }         
 
336
         }
 
337
         else 
 
338
           for(uint i = 0; i < tt; i++) {
 
339
             for(uint k = 0; k < ufc_dof_map->geometric_dimension(); k++) {
 
340
               const uint dof = dofs[i + k * tt];
 
341
               
 
342
               std::map<uint, uint>::iterator it = map.find(dof);
 
343
               if (it != map.end())
 
344
                 dof_map[c->index()][i + k * tt] = it->second;
 
345
               else {
 
346
                 dof_map[c->index()][i + k * tt] = num;
 
347
                 map[dof] = num++;
 
348
#ifdef BLOCKED
 
349
                 dof_map[c->index()][i + k * tt] /= 2;
 
350
                 map[dof] /= 2;
 
351
#endif
 
352
               }
 
353
             }
 
354
           }    
 
355
       }
 
356
     }
 
357
     
 
358
 
 
359
//             for(CellIterator c(dolfin_mesh); !c.end(); ++c) {
 
360
//         ufc.update(*c);
 
361
//          for(uint j = 0; j< ufc.form.rank(); j++) {
 
362
//      tabulate_dofs(dofs, ufc.cell, c->index());
 
363
//       for(uint i = 0; i < local_dimension(); i++)
 
364
//       cout<< dofs[i] << " ";
 
365
//       cout<<endl;
 
366
//       }
 
367
 
 
368
 
 
369
         //   }
 
370
         */
 
371
    //     delete[] dofs;
 
372
  }
 
373
  else { 
 
374
#ifdef HAVE_MPI
 
375
    uint *dofs =  new uint[local_dimension()];
 
376
 
 
377
    uint pe_size = MPI::numProcesses();
 
378
    uint rank = MPI::processNumber();
 
379
 
 
380
    dolfin_mesh.renumber();
 
381
 
 
382
    if (ufc_dof_map->global_dimension() == dolfin_mesh.distdata().global_numVertices()) {
 
383
      _type_ = 0;
 
384
      _local_size = dolfin_mesh.numVertices() - dolfin_mesh.distdata().num_ghost(0);
 
385
    }
 
386
    else if(ufc_dof_map->global_dimension() == dolfin_mesh.distdata().global_numCells()) {
 
387
      _type_ = 1;      
 
388
      _local_size = dolfin_mesh.numCells();
 
389
    }
 
390
    else if(ufc_dof_map->global_dimension() == 
 
391
       ufc_dof_map->geometric_dimension() * dolfin_mesh.distdata().global_numVertices()) {
 
392
      
 
393
      uint gdim = ufc_dof_map->geometric_dimension();
 
394
      uint num_local = dolfin_mesh.numVertices() - dolfin_mesh.distdata().num_ghost(0);
 
395
      
 
396
      uint num_dofs = gdim * num_local;
 
397
      uint offset = 0;
 
398
 
 
399
#if ( MPI_VERSION > 1 )
 
400
      MPI_Exscan(&num_dofs, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
401
#else 
 
402
      MPI_Scan(&num_dofs, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
403
      offset -= num_dofs;
 
404
#endif
 
405
      _map<uint, uint> v_offset;
 
406
 
 
407
      for(VertexIterator v(dolfin_mesh); !v.end(); ++v) {
 
408
        if(!dolfin_mesh.distdata().is_ghost(v->index(), 0)) {
 
409
          v_offset[dolfin_mesh.distdata().get_global(*v)] = offset; 
 
410
          offset += gdim;
 
411
        }
 
412
      }
 
413
 
 
414
      Array<uint> *ghost_buff = new Array<uint>[pe_size];
 
415
      for(MeshGhostIterator iter(dolfin_mesh.distdata(), 0); !iter.end(); ++iter)
 
416
        ghost_buff[iter.owner()].push_back(dolfin_mesh.distdata().get_global(iter.index(), 0)); 
 
417
                
 
418
      MPI_Status status;
 
419
      Array<uint> send_buff;
 
420
      uint src,dest;
 
421
      uint recv_size = dolfin_mesh.distdata().num_ghost(0); 
 
422
      int recv_count, recv_size_gh, send_size;  
 
423
      
 
424
      for(uint i = 0; i < pe_size; i++) {
 
425
        send_size = ghost_buff[i].size();
 
426
        MPI_Reduce(&send_size, &recv_size_gh, 1, 
 
427
                   MPI_INT, MPI_SUM, i, MPI::DOLFIN_COMM);
 
428
      }
 
429
      
 
430
      uint *recv_ghost = new uint[ recv_size_gh];
 
431
      uint *recv_buff = new uint[ recv_size ];
 
432
      
 
433
      for(uint j=1; j < pe_size; j++){
 
434
        src = (rank - j + pe_size) % pe_size;
 
435
        dest = (rank + j) % pe_size;
 
436
        
 
437
        MPI_Sendrecv(&ghost_buff[dest][0], ghost_buff[dest].size(),
 
438
                     MPI_UNSIGNED, dest, 1, recv_ghost, recv_size_gh, 
 
439
                     MPI_UNSIGNED, src, 1, MPI::DOLFIN_COMM, &status);
 
440
        MPI_Get_count(&status,MPI_UNSIGNED,&recv_count);
 
441
        
 
442
        for(int k=0; k < recv_count; k++)
 
443
          send_buff.push_back(v_offset[recv_ghost[k]]);
 
444
        
 
445
        MPI_Sendrecv(&send_buff[0], send_buff.size(), MPI_UNSIGNED, src, 2,
 
446
                     recv_buff, recv_size , MPI_UNSIGNED, dest, 2, 
 
447
                     MPI::DOLFIN_COMM,&status);
 
448
        MPI_Get_count(&status,MPI_UNSIGNED,&recv_count);
 
449
 
 
450
        for(int j=0; j < recv_count; j++)
 
451
          v_offset[ghost_buff[dest][j]] = recv_buff[j];
 
452
 
 
453
        send_buff.clear();
 
454
      }
 
455
 
 
456
      delete[] recv_ghost;
 
457
      delete[] recv_buff;
 
458
 
 
459
      if( v_map )
 
460
        delete[] v_map;
 
461
 
 
462
      v_map = new uint[dolfin_mesh.numVertices()];
 
463
      for(VertexIterator v(dolfin_mesh); !v.end(); ++v)
 
464
        v_map[v->index()] = v_offset[dolfin_mesh.distdata().get_global(v->index(), 0)];
 
465
      
 
466
      _type_ = 2;
 
467
      v_offset.clear();
 
468
      
 
469
      _local_size = gdim * num_local;
 
470
 
 
471
      for(uint i =0; i < pe_size; i++)
 
472
        ghost_buff[i].clear();
 
473
      delete[] ghost_buff;
 
474
      
 
475
    }
 
476
    else if(ufc_dof_map->global_dimension() == 
 
477
       ufc_dof_map->geometric_dimension() * dolfin_mesh.distdata().global_numCells()) {
 
478
 
 
479
      uint gdim = ufc_dof_map->geometric_dimension();
 
480
      uint num_dofs =  gdim * dolfin_mesh.numCells();      
 
481
      uint offset = 0;
 
482
 
 
483
#if ( MPI_VERSION > 1 )
 
484
      MPI_Exscan(&num_dofs, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
485
#else 
 
486
      MPI_Scan(&num_dofs, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
487
      offset -= num_dofs;
 
488
#endif
 
489
      _offset_ = offset;
 
490
      _type_ = 3;
 
491
      _local_size =  gdim * dolfin_mesh.numCells();      
 
492
    }
 
493
    else {
 
494
 
 
495
 
 
496
      BoundaryMesh interior_boundary;
 
497
      interior_boundary.init_interior(dolfin_mesh);
 
498
      MeshFunction<uint>* cell_map =interior_boundary.data().meshFunction("cell map");
 
499
 
 
500
      std::vector<uint> send_buffer;                                                
 
501
      _set<uint> shared_dofs, forbidden_dofs, owned_dofs;                            
 
502
      _map<uint, uint> dof_vote;                                                
 
503
      _map<uint, std::vector<uint> > dof2index;
 
504
 
 
505
      uint n = local_dimension(); 
 
506
      dof_map = new uint[n * dolfin_mesh.numCells()];
 
507
      uint *facet_dofs = new uint[num_facet_dofs()];
 
508
 
 
509
      Cell c_tmp(dolfin_mesh, 1);
 
510
      UFCCell ufc_cell(c_tmp);      
 
511
                                                                                      
 
512
      // Initialize random number generator differently on each process             
 
513
      srand((uint)time(0) + MPI::processNumber());   
 
514
 
 
515
      // Decide ownership of shared dofs                                            
 
516
      for (CellIterator bc(interior_boundary); !bc.end(); ++bc)                     
 
517
      {                                                                             
 
518
        Facet f(dolfin_mesh, cell_map->get(*bc));
 
519
        Cell c(dolfin_mesh, f.entities(dolfin_mesh.topology().dim())[0]);
 
520
 
 
521
        uint local_facet = c.index(f);                                          
 
522
 
 
523
        ufc_cell.update(c, dolfin_mesh.distdata());                                  
 
524
        ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);                        
 
525
        ufc_dof_map->tabulate_facet_dofs(facet_dofs, local_facet);
 
526
 
 
527
        for (uint i = 0; i < num_facet_dofs(); i++)    
 
528
        {                                                                         
 
529
          // Assign an ownership vote for each "shared" dof                       
 
530
          if (shared_dofs.find(dofs[i]) == shared_dofs.end())                     
 
531
          {                                                                       
 
532
            shared_dofs.insert(dofs[i]);                                          
 
533
            dof_vote[dofs[i]] = (uint) rand() + (uint) MPI::processNumber();
 
534
            send_buffer.push_back(dofs[i]);                                       
 
535
            send_buffer.push_back(dof_vote[dofs[i]]);                             
 
536
          }                                                                       
 
537
        }                                                                         
 
538
      }  
 
539
 
 
540
      // Decide ownership of "shared" dofs                                          
 
541
      MPI_Status status;
 
542
      int recv_count;
 
543
      uint src, dest, max_recv;
 
544
      uint num_proc = MPI::numProcesses();                                         
 
545
      uint proc_num = MPI::processNumber();                                        
 
546
      uint local_size = send_buffer.size();
 
547
      MPI_Allreduce(&local_size, &max_recv, 1,
 
548
                   MPI_UNSIGNED, MPI_MAX, MPI::DOLFIN_COMM);
 
549
      uint *recv_buffer = new uint[max_recv];                                       
 
550
      for(uint k = 1; k < MPI::numProcesses(); ++k)                                
 
551
      {                                                                             
 
552
        src = (proc_num - k + num_proc) % num_proc;                                 
 
553
        dest = (proc_num +k) % num_proc;                                            
 
554
        
 
555
        MPI_Sendrecv(&send_buffer[0], send_buffer.size(), MPI_UNSIGNED, dest, 1,
 
556
                     recv_buffer, max_recv, MPI_UNSIGNED, src, 1,
 
557
                     MPI::DOLFIN_COMM, &status);
 
558
        MPI_Get_count(&status, MPI_UNSIGNED, &recv_count);
 
559
        
 
560
        for (int i = 0; i < recv_count; i += 2)                                    
 
561
        {                                                                           
 
562
          if (shared_dofs.find(recv_buffer[i]) != shared_dofs.end())                
 
563
          {                                                                         
 
564
            // Move dofs with higher ownership votes from shared to forbidden       
 
565
            if (recv_buffer[i+1] < dof_vote[recv_buffer[i]] ) {                     
 
566
              forbidden_dofs.insert(recv_buffer[i]);                                
 
567
              shared_dofs.erase(recv_buffer[i]);                                    
 
568
            }                                                                       
 
569
          }                                                                         
 
570
        }                                                                           
 
571
      }                                                                             
 
572
      
 
573
      send_buffer.clear();  
 
574
      
 
575
      // Mark all non forbidden dofs as owned by the processes                      
 
576
      for (CellIterator c(dolfin_mesh); !c.end(); ++c)                                     
 
577
      {                                                                             
 
578
        ufc_cell.update(*c, dolfin_mesh.distdata());                                  
 
579
        ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);  
 
580
        for (uint i = 0; i < n; i++)                                                
 
581
        {                                                                           
 
582
          if (forbidden_dofs.find(dofs[i]) == forbidden_dofs.end())                 
 
583
          {                                                                         
 
584
            // Mark dof as owned                                                    
 
585
            owned_dofs.insert(dofs[i]);                                             
 
586
          }                                                                         
 
587
          
 
588
          // Create mapping from dof to dof_map offset                              
 
589
          dof2index[dofs[i]].push_back(c->index() * n + i);                         
 
590
        }                                                                           
 
591
      }                                                                             
 
592
      
 
593
      // Compute offset for owned and non shared dofs                               
 
594
      uint range = owned_dofs.size();                                               
 
595
      _local_size = range;
 
596
      uint offset = 0;
 
597
#if ( MPI_VERSION > 1 )
 
598
      MPI_Exscan(&range, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
599
#else 
 
600
      MPI_Scan(&range, &offset, 1, MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
601
      offset -= range;
 
602
#endif    
 
603
 
 
604
      // Compute renumbering for local and owned shared dofs                        
 
605
      for (_set<uint>::iterator it = owned_dofs.begin();                        
 
606
           it != owned_dofs.end(); ++it, offset++)                                  
 
607
      {                                                                             
 
608
        for(std::vector<uint>::iterator di = dof2index[*it].begin();                
 
609
            di != dof2index[*it].end(); ++di)                                       
 
610
          dof_map[*di] = offset;                                                   
 
611
        
 
612
        if (shared_dofs.find(*it) != shared_dofs.end())                             
 
613
        {                                                                           
 
614
          send_buffer.push_back(*it);                                               
 
615
          send_buffer.push_back(offset);                                            
 
616
        }                                                                           
 
617
      }                                                                             
 
618
      
 
619
 
 
620
      // Exchange new dof numbers for shared dofs                                   
 
621
      delete[] recv_buffer;                                                         
 
622
      local_size = send_buffer.size();
 
623
      MPI_Allreduce(&local_size, &max_recv, 1,
 
624
                   MPI_UNSIGNED, MPI_MAX, MPI::DOLFIN_COMM);     
 
625
      recv_buffer = new uint[max_recv];                                             
 
626
      for(uint k = 1; k < MPI::numProcesses(); ++k)                                
 
627
      {                                                                             
 
628
        src = (proc_num - k + num_proc) % num_proc;                                 
 
629
        dest = (proc_num +k) % num_proc;                                            
 
630
        
 
631
        MPI_Sendrecv(&send_buffer[0], send_buffer.size(), MPI_UNSIGNED, dest, 1,
 
632
                     recv_buffer, max_recv, MPI_UNSIGNED, src, 1,
 
633
                     MPI::DOLFIN_COMM, &status);
 
634
        MPI_Get_count(&status, MPI_UNSIGNED, &recv_count);
 
635
        
 
636
        for (int i = 0; i < recv_count; i += 2)                                    
 
637
        {                                                                           
 
638
          // Assign new dof number for shared dofs                                  
 
639
          if (forbidden_dofs.find(recv_buffer[i]) != forbidden_dofs.end())          
 
640
          {                                                                         
 
641
            for(std::vector<uint>::iterator di = dof2index[recv_buffer[i]].begin(); 
 
642
                di != dof2index[recv_buffer[i]].end(); ++di)                        
 
643
              dof_map[*di] = recv_buffer[i+1];                                       
 
644
          }                                                                         
 
645
        }                                                                           
 
646
      }                                                                             
 
647
      delete[] recv_buffer;                                                         
 
648
      map.clear();
 
649
      delete[] facet_dofs;
 
650
    }
 
651
    delete[] dofs;   
 
652
 
 
653
#endif
 
654
  }
 
655
}
 
656
//-----------------------------------------------------------------------------
 
657
std::map<dolfin::uint, dolfin::uint> DofMap::getMap() //FIXME: const
242
658
{
243
659
  return map;
244
660
}
353
769
    UFCCell ufc_cell(*cell);
354
770
    for (; !cell.end(); ++cell)
355
771
    {
356
 
      ufc_cell.update(*cell);
 
772
      ufc_cell.update(*cell, dolfin_mesh.distdata());
357
773
 
358
774
      ufc_dof_map->tabulate_dofs(dofs, ufc_mesh, ufc_cell);
359
775
 
386
802
    UFCCell ufc_cell(*cell);
387
803
    for (; !cell.end(); ++cell)
388
804
    {
389
 
      ufc_cell.update(*cell);
 
805
      ufc_cell.update(*cell, dolfin_mesh.distdata());
390
806
 
391
807
      ufc_dof_map->tabulate_coordinates(coordinates, ufc_cell);
392
808