← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/mesh/TopologyComputation.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/mesh/TopologyComputation.cpp
1
 
// Copyright (C) 2006-2007 Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// First added:  2006-06-02
5
 
// Last changed: 2007-03-13
6
 
 
7
 
#include <set>
8
 
#include <dolfin/dolfin_log.h>
9
 
#include <dolfin/Array.h>
10
 
#include <dolfin/CellType.h>
11
 
#include <dolfin/Mesh.h>
12
 
#include <dolfin/MeshTopology.h>
13
 
#include <dolfin/MeshConnectivity.h>
14
 
#include <dolfin/MeshEntity.h>
15
 
#include <dolfin/MeshEntityIterator.h>
16
 
#include <dolfin/TopologyComputation.h>
17
 
 
18
 
using namespace dolfin;
19
 
 
20
 
//-----------------------------------------------------------------------------
21
 
dolfin::uint TopologyComputation::computeEntities(Mesh& mesh, uint dim)
22
 
{
23
 
  // Generating an entity of topological dimension dim is equivalent
24
 
  // to generating the connectivity dim - 0 (connections to vertices)
25
 
  // and the connectivity mesh.topology().dim() - dim (connections from cells).
26
 
  //
27
 
  // We generate entities by iterating over all cells and generating a
28
 
  // new entity only on its first occurence. Entities also contained
29
 
  // in a previously visited cell are not generated. The new entities
30
 
  // are computed in three steps:
31
 
  //
32
 
  //   1. Iterate over cells and count new entities
33
 
  //
34
 
  //   2. Allocate memory / prepare data structures
35
 
  //
36
 
  //   3. Iterate over cells and add new entities
37
 
  
38
 
  // Get mesh topology and connectivity
39
 
  MeshTopology& topology = mesh.topology();
40
 
 
41
 
  MeshConnectivity& ce = topology(topology.dim(), dim);
42
 
  MeshConnectivity& ev = topology(dim, 0);
43
 
 
44
 
  // Check if entities have already been computed
45
 
  if ( topology.size(dim) > 0 )
46
 
  {
47
 
    // Make sure we really have the connectivity
48
 
    if ( (ce.size() == 0 && dim != topology.dim()) || (ev.size() == 0 && dim != 0) )
49
 
      error("Entities of topological dimension %d exist but connectivity is missing.", dim);
50
 
    return topology.size(dim);
51
 
  }
52
 
  else
53
 
  {
54
 
    // Make sure connectivity does not already exist
55
 
    if ( ce.size() > 0 || ev.size() > 0 )
56
 
      error("Connectivity for topological dimension %d exists but entities are missing.", dim);
57
 
  }
58
 
 
59
 
  //message("Computing mesh entities of topological dimension %d.", dim);
60
 
 
61
 
  // Compute connectivity dim - dim if not already computed
62
 
  computeConnectivity(mesh, mesh.topology().dim(), mesh.topology().dim());
63
 
 
64
 
  // Get cell type
65
 
  const CellType& cell_type = mesh.type();
66
 
 
67
 
  // Initialize local array of entities
68
 
  const uint m = cell_type.numEntities(dim);
69
 
  const uint n = cell_type.numVertices(dim);
70
 
  uint** entities = new uint*[m];
71
 
  for (uint i = 0; i < m; i++)
72
 
  {
73
 
    entities[i] = new uint[n];
74
 
    for (uint j = 0; j < n; j++)
75
 
      entities[i][j] = 0;
76
 
  }
77
 
 
78
 
  // Count the number of entities
79
 
  uint num_entities = 0;
80
 
  for (MeshEntityIterator c(mesh, mesh.topology().dim()); !c.end(); ++c)
81
 
  {
82
 
    // Get vertices from cell
83
 
    const uint* vertices = c->entities(0);
84
 
    dolfin_assert(vertices);
85
 
    
86
 
    // Create entities
87
 
    cell_type.createEntities(entities, dim, vertices);
88
 
    
89
 
    // Count new entities
90
 
    num_entities += countEntities(mesh, *c, entities, m, n, dim);
91
 
  }
92
 
 
93
 
  // Initialize the number of entities and connections
94
 
  topology.init(dim, num_entities);
95
 
  ce.init(mesh.numCells(), m);
96
 
  ev.init(num_entities, n);
97
 
 
98
 
  // Add new entities
99
 
  uint current_entity = 0;
100
 
  for (MeshEntityIterator c(mesh, mesh.topology().dim()); !c.end(); ++c)
101
 
  {
102
 
    // Get vertices from cell
103
 
    const uint* vertices = c->entities(0);
104
 
    dolfin_assert(vertices);
105
 
    
106
 
    // Create entities
107
 
    cell_type.createEntities(entities, dim, vertices);
108
 
    
109
 
    // Add new entities to the mesh
110
 
    addEntities(mesh, *c, entities, m, n, dim, ce, ev, current_entity);
111
 
  }
112
 
 
113
 
  // Delete temporary data
114
 
  for (uint i = 0; i < m; i++)
115
 
    delete [] entities[i];
116
 
  delete [] entities;
117
 
 
118
 
  //message("Created %d new entities.", num_entities);
119
 
 
120
 
  return num_entities;
121
 
}
122
 
//-----------------------------------------------------------------------------
123
 
void TopologyComputation::computeConnectivity(Mesh& mesh, uint d0, uint d1)
124
 
{
125
 
  // This is where all the logic takes place to find a stragety for
126
 
  // the connectivity computation. For any given pair (d0, d1), the
127
 
  // connectivity is computed by suitably combining the following
128
 
  // basic building blocks:
129
 
  //
130
 
  //   1. computeEntities():     d  - 0  from dim - 0
131
 
  //   2. computeTranspose():    d0 - d1 from d1 - d0
132
 
  //   3. computeIntersection(): d0 - d1 from d0 - d' - d1
133
 
  //
134
 
  // Each of these functions assume a set of preconditions that we
135
 
  // need to satisfy.
136
 
 
137
 
  // Get mesh topology and connectivity
138
 
  MeshTopology& topology = mesh.topology();
139
 
  MeshConnectivity& connectivity = topology(d0, d1);
140
 
 
141
 
  // Check if connectivity has already been computed
142
 
  if ( connectivity.size() > 0 )
143
 
    return;
144
 
 
145
 
  //message("Computing mesh connectivity %d - %d.", d0, d1);
146
 
 
147
 
  // Compute entities if they don't exist
148
 
  if ( topology.size(d0) == 0 )
149
 
    computeEntities(mesh, d0);
150
 
  if ( topology.size(d1) == 0 )
151
 
    computeEntities(mesh, d1);
152
 
 
153
 
  // Check if connectivity still needs to be computed
154
 
  if ( connectivity.size() > 0 )
155
 
    return;
156
 
 
157
 
  // Decide how to compute the connectivity
158
 
  if ( d0 < d1 )
159
 
  {
160
 
    // Compute connectivity d1 - d0 and take transpose
161
 
    computeConnectivity(mesh, d1, d0);
162
 
    computeFromTranspose(mesh, d0, d1);
163
 
  }
164
 
  else
165
 
  {
166
 
    // These connections should already exist
167
 
    dolfin_assert(!(d0 > 0 && d1 == 0));
168
 
 
169
 
    // Choose how to take intersection
170
 
    uint d = 0;
171
 
    if ( d0 == 0 && d1 == 0 )
172
 
      d = mesh.topology().dim();
173
 
 
174
 
    // Compute connectivity d0 - d - d1 and take intersection
175
 
    computeConnectivity(mesh, d0, d);
176
 
    computeConnectivity(mesh, d, d1);
177
 
    computeFromIntersection(mesh, d0, d1, d);
178
 
  }
179
 
}
180
 
//----------------------------------------------------------------------------
181
 
void TopologyComputation::computeFromTranspose(Mesh& mesh, uint d0, uint d1)
182
 
{
183
 
  // The transpose is computed in three steps:
184
 
  //
185
 
  //   1. Iterate over entities of dimension d1 and count the number
186
 
  //      of connections for each entity of dimension d0
187
 
  //
188
 
  //   2. Allocate memory / prepare data structures
189
 
  //
190
 
  //   3. Iterate again over entities of dimension d1 and add connections
191
 
  //      for each entity of dimension d0
192
 
 
193
 
  //message("Computing mesh connectivity %d - %d from transpose.", d0, d1);
194
 
  
195
 
  // Get mesh topology and connectivity
196
 
  MeshTopology& topology = mesh.topology();
197
 
  MeshConnectivity& connectivity = topology(d0, d1);
198
 
 
199
 
  // Need connectivity d1 - d0
200
 
  dolfin_assert(topology(d1, d0).size() > 0);
201
 
 
202
 
  // Temporary array
203
 
  Array<uint> tmp(topology.size(d0));
204
 
 
205
 
  // Reset size for each entity
206
 
  for (uint i = 0; i < tmp.size(); i++)
207
 
    tmp[i] = 0;
208
 
 
209
 
  // Count the number of connections
210
 
  for (MeshEntityIterator e1(mesh, d1); !e1.end(); ++e1)
211
 
    for (MeshEntityIterator e0(*e1, d0); !e0.end(); ++e0)
212
 
      tmp[e0->index()]++;
213
 
 
214
 
  // Initialize the number of connections
215
 
  connectivity.init(tmp);
216
 
 
217
 
  // Reset current position for each entity
218
 
  for (uint i = 0; i < tmp.size(); i++)
219
 
    tmp[i] = 0;
220
 
  
221
 
  // Add the connections
222
 
  for (MeshEntityIterator e1(mesh, d1); !e1.end(); ++e1)
223
 
    for (MeshEntityIterator e0(*e1, d0); !e0.end(); ++e0)
224
 
      connectivity.set(e0->index(), e1->index(), tmp[e0->index()]++);
225
 
}
226
 
//----------------------------------------------------------------------------
227
 
void TopologyComputation::computeFromIntersection(Mesh& mesh,
228
 
                                             uint d0, uint d1, uint d)
229
 
{
230
 
  // The intersection is computed in three steps:
231
 
  //
232
 
  //   1. Nested iteration over mesh - d0 - d - d1 and count the connections
233
 
  //
234
 
  //   2. Allocate memory / prepare data structures
235
 
  //
236
 
  //   3. Nested iteration over mesh - d0 - d - d1 and add the connections
237
 
 
238
 
  //message("Computing mesh connectivity %d - %d from intersection %d - %d - %d.",
239
 
  //            d0, d1, d0, d, d1);
240
 
 
241
 
  // Get mesh topology and connectivity
242
 
  MeshTopology& topology = mesh.topology();
243
 
  MeshConnectivity& connectivity = topology(d0, d1);
244
 
 
245
 
  // Need d0 >= d1
246
 
  dolfin_assert(d0 >= d1);
247
 
 
248
 
  // Need connectivity d0 - d and d - d1
249
 
  dolfin_assert(topology(d0, d).size() > 0);
250
 
  dolfin_assert(topology(d, d1).size() > 0);
251
 
 
252
 
  // Temporary array
253
 
  Array<uint> tmp(topology.size(d0));
254
 
 
255
 
  // Reset size for each entity
256
 
  for (uint i = 0; i < tmp.size(); i++)
257
 
    tmp[i] = 0;
258
 
 
259
 
  // FIXME: Check efficiency of std::set, compare with std::vector
260
 
 
261
 
  // A set with connected entities
262
 
  std::set<uint> entities;
263
 
 
264
 
  // Iterate over all entities of dimension d0
265
 
  for (MeshEntityIterator e0(mesh, d0); !e0.end(); ++e0)
266
 
  {
267
 
    // Clear set of connected entities
268
 
    entities.clear();
269
 
 
270
 
    // Iterate over all connected entities of dimension d
271
 
    for (MeshEntityIterator e(*e0, d); !e.end(); ++e)
272
 
    {
273
 
      // Iterate over all connected entities of dimension d1
274
 
      for (MeshEntityIterator e1(*e, d1); !e1.end(); ++e1)
275
 
      {
276
 
        if ( d0 == d1 )
277
 
        {
278
 
          // An entity is not a neighbor to itself
279
 
          if ( e0->index() != e1->index() )
280
 
            entities.insert(e1->index());
281
 
        }
282
 
        else
283
 
        {
284
 
          // Entity e1 must be completely contained in e0
285
 
          if ( contains(*e0, *e1) )
286
 
            entities.insert(e1->index());
287
 
        }
288
 
      }
289
 
    }
290
 
 
291
 
    // Count the number of connected entities
292
 
    tmp[e0->index()] = entities.size();
293
 
  }
294
 
 
295
 
  // Initialize the number of connections
296
 
  connectivity.init(tmp);
297
 
 
298
 
  // Iterate over all entities of dimension d
299
 
  for (MeshEntityIterator e0(mesh, d0); !e0.end(); ++e0)
300
 
  {
301
 
    // Clear set of connected entities
302
 
    entities.clear();
303
 
 
304
 
    // Iterate over all connected entities of dimension d
305
 
    for (MeshEntityIterator e(*e0, d); !e.end(); ++e)
306
 
    {
307
 
      // Iterate over all connected entities of dimension d1
308
 
      for (MeshEntityIterator e1(*e, d1); !e1.end(); ++e1)
309
 
      {
310
 
        if ( d0 == d1 )
311
 
        {
312
 
          // An entity is not a neighbor to itself
313
 
          if ( e0->index() != e1->index() )
314
 
            entities.insert(e1->index());
315
 
        }
316
 
        else
317
 
        {
318
 
          // Entity e1 must be completely contained in e0
319
 
          if ( contains(*e0, *e1) )
320
 
            entities.insert(e1->index());
321
 
        }
322
 
      }
323
 
    }
324
 
 
325
 
    // Add the connected entities
326
 
    uint pos = 0;
327
 
    for (std::set<uint>::iterator it = entities.begin(); it != entities.end(); ++it)
328
 
      connectivity.set(e0->index(), *it, pos++);
329
 
  }
330
 
}
331
 
//----------------------------------------------------------------------------
332
 
dolfin::uint TopologyComputation::countEntities(Mesh& mesh, MeshEntity& cell,
333
 
                                           uint** entities, uint m, uint n,
334
 
                                           uint dim)
335
 
{
336
 
  // For each entity, we iterate over connected and previously visited
337
 
  // cells to see if the entity has already been counted.
338
 
  
339
 
  // Needs to be a cell
340
 
  dolfin_assert(cell.dim() == mesh.topology().dim());
341
 
 
342
 
  // Iterate over the given list of entities
343
 
  uint num_entities = 0;
344
 
  for (uint i = 0; i < m; i++)
345
 
  {
346
 
    // Iterate over connected cells and look for entity
347
 
    for (MeshEntityIterator c(cell, mesh.topology().dim()); !c.end(); ++c)
348
 
    {
349
 
      // Check only previously visited cells
350
 
      if ( c->index() >= cell.index() )
351
 
        continue;
352
 
 
353
 
      // Check for vertices
354
 
      if ( contains(c->entities(0), c->numEntities(0), entities[i], n) )
355
 
        goto found;
356
 
    }
357
 
    
358
 
    // Increase counter
359
 
    num_entities++;
360
 
    
361
 
    // Entity found, don't need to count
362
 
  found:
363
 
    ;
364
 
  }
365
 
 
366
 
  return num_entities;
367
 
}
368
 
//----------------------------------------------------------------------------
369
 
void TopologyComputation::addEntities(Mesh& mesh, MeshEntity& cell,
370
 
                                 uint** entities, uint m, uint n, uint dim,
371
 
                                 MeshConnectivity& ce,
372
 
                                 MeshConnectivity& ev,
373
 
                                 uint& current_entity)
374
 
{
375
 
  // We repeat the same algorithm as in countEntities() but this time
376
 
  // we add any entities that are new.
377
 
  
378
 
  // Needs to be a cell
379
 
  dolfin_assert(cell.dim() == mesh.topology().dim());
380
 
 
381
 
  // Iterate over the given list of entities
382
 
  for (uint i = 0; i < m; i++)
383
 
  {
384
 
    // Iterate over connected cells and look for entity
385
 
    for (MeshEntityIterator c(cell, mesh.topology().dim()); !c.end(); ++c)
386
 
    {
387
 
      // Check only previously visited cells
388
 
      if ( c->index() >= cell.index() )
389
 
        continue;
390
 
      
391
 
      // Check all entities of dimension dim in connected cell
392
 
      uint num_other_entities = c->numEntities(dim);
393
 
      uint* other_entities = c->entities(dim);
394
 
      for (uint j = 0; j < num_other_entities; j++)
395
 
      {
396
 
        // Can't use iterators since connectivity has not been computed
397
 
        MeshEntity e(mesh, dim, other_entities[j]);
398
 
        if ( contains(e.entities(0), e.numEntities(0), entities[i], n) )
399
 
        {
400
 
          // Entity already exists, so pick the index
401
 
          ce.set(cell.index(), e.index(), i);
402
 
          goto found;
403
 
        }
404
 
      }
405
 
    }
406
 
    
407
 
    // Entity does not exist, so create it
408
 
    ce.set(cell.index(), current_entity, i);
409
 
    ev.set(current_entity, entities[i]);
410
 
    
411
 
    // Increase counter
412
 
    current_entity++;
413
 
    
414
 
    // Entity found, don't need to create
415
 
  found:
416
 
    ;
417
 
  }
418
 
}
419
 
//----------------------------------------------------------------------------
420
 
bool TopologyComputation::contains(MeshEntity& e0, MeshEntity& e1)
421
 
{
422
 
  // Check vertices
423
 
  return contains(e0.entities(0), e0.numEntities(0),
424
 
                  e1.entities(0), e1.numEntities(0));
425
 
}
426
 
//----------------------------------------------------------------------------
427
 
bool TopologyComputation::contains(uint* v0, uint n0, uint* v1, uint n1)
428
 
{
429
 
  dolfin_assert(v0);
430
 
  dolfin_assert(v1);
431
 
 
432
 
  for (uint i1 = 0; i1 < n1; i1++)
433
 
  {
434
 
    bool found = false;
435
 
    for (uint i0 = 0; i0 < n0; i0++)
436
 
    {
437
 
      if ( v0[i0] == v1[i1] )
438
 
      {
439
 
        found = true;
440
 
        break;
441
 
      }
442
 
    }
443
 
    if ( !found )
444
 
      return false;
445
 
  }
446
 
 
447
 
  return true;
448
 
}
449
 
//----------------------------------------------------------------------------