1
// Copyright (C) 2006-2007 Anders Logg.
2
// Licensed under the GNU LGPL Version 2.1.
4
// First added: 2006-06-02
5
// Last changed: 2007-03-13
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>
18
using namespace dolfin;
20
//-----------------------------------------------------------------------------
21
dolfin::uint TopologyComputation::computeEntities(Mesh& mesh, uint dim)
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).
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:
32
// 1. Iterate over cells and count new entities
34
// 2. Allocate memory / prepare data structures
36
// 3. Iterate over cells and add new entities
38
// Get mesh topology and connectivity
39
MeshTopology& topology = mesh.topology();
41
MeshConnectivity& ce = topology(topology.dim(), dim);
42
MeshConnectivity& ev = topology(dim, 0);
44
// Check if entities have already been computed
45
if ( topology.size(dim) > 0 )
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);
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);
59
//message("Computing mesh entities of topological dimension %d.", dim);
61
// Compute connectivity dim - dim if not already computed
62
computeConnectivity(mesh, mesh.topology().dim(), mesh.topology().dim());
65
const CellType& cell_type = mesh.type();
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++)
73
entities[i] = new uint[n];
74
for (uint j = 0; j < n; j++)
78
// Count the number of entities
79
uint num_entities = 0;
80
for (MeshEntityIterator c(mesh, mesh.topology().dim()); !c.end(); ++c)
82
// Get vertices from cell
83
const uint* vertices = c->entities(0);
84
dolfin_assert(vertices);
87
cell_type.createEntities(entities, dim, vertices);
90
num_entities += countEntities(mesh, *c, entities, m, n, dim);
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);
99
uint current_entity = 0;
100
for (MeshEntityIterator c(mesh, mesh.topology().dim()); !c.end(); ++c)
102
// Get vertices from cell
103
const uint* vertices = c->entities(0);
104
dolfin_assert(vertices);
107
cell_type.createEntities(entities, dim, vertices);
109
// Add new entities to the mesh
110
addEntities(mesh, *c, entities, m, n, dim, ce, ev, current_entity);
113
// Delete temporary data
114
for (uint i = 0; i < m; i++)
115
delete [] entities[i];
118
//message("Created %d new entities.", num_entities);
122
//-----------------------------------------------------------------------------
123
void TopologyComputation::computeConnectivity(Mesh& mesh, uint d0, uint d1)
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:
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
134
// Each of these functions assume a set of preconditions that we
137
// Get mesh topology and connectivity
138
MeshTopology& topology = mesh.topology();
139
MeshConnectivity& connectivity = topology(d0, d1);
141
// Check if connectivity has already been computed
142
if ( connectivity.size() > 0 )
145
//message("Computing mesh connectivity %d - %d.", d0, d1);
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);
153
// Check if connectivity still needs to be computed
154
if ( connectivity.size() > 0 )
157
// Decide how to compute the connectivity
160
// Compute connectivity d1 - d0 and take transpose
161
computeConnectivity(mesh, d1, d0);
162
computeFromTranspose(mesh, d0, d1);
166
// These connections should already exist
167
dolfin_assert(!(d0 > 0 && d1 == 0));
169
// Choose how to take intersection
171
if ( d0 == 0 && d1 == 0 )
172
d = mesh.topology().dim();
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);
180
//----------------------------------------------------------------------------
181
void TopologyComputation::computeFromTranspose(Mesh& mesh, uint d0, uint d1)
183
// The transpose is computed in three steps:
185
// 1. Iterate over entities of dimension d1 and count the number
186
// of connections for each entity of dimension d0
188
// 2. Allocate memory / prepare data structures
190
// 3. Iterate again over entities of dimension d1 and add connections
191
// for each entity of dimension d0
193
//message("Computing mesh connectivity %d - %d from transpose.", d0, d1);
195
// Get mesh topology and connectivity
196
MeshTopology& topology = mesh.topology();
197
MeshConnectivity& connectivity = topology(d0, d1);
199
// Need connectivity d1 - d0
200
dolfin_assert(topology(d1, d0).size() > 0);
203
Array<uint> tmp(topology.size(d0));
205
// Reset size for each entity
206
for (uint i = 0; i < tmp.size(); i++)
209
// Count the number of connections
210
for (MeshEntityIterator e1(mesh, d1); !e1.end(); ++e1)
211
for (MeshEntityIterator e0(*e1, d0); !e0.end(); ++e0)
214
// Initialize the number of connections
215
connectivity.init(tmp);
217
// Reset current position for each entity
218
for (uint i = 0; i < tmp.size(); i++)
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()]++);
226
//----------------------------------------------------------------------------
227
void TopologyComputation::computeFromIntersection(Mesh& mesh,
228
uint d0, uint d1, uint d)
230
// The intersection is computed in three steps:
232
// 1. Nested iteration over mesh - d0 - d - d1 and count the connections
234
// 2. Allocate memory / prepare data structures
236
// 3. Nested iteration over mesh - d0 - d - d1 and add the connections
238
//message("Computing mesh connectivity %d - %d from intersection %d - %d - %d.",
239
// d0, d1, d0, d, d1);
241
// Get mesh topology and connectivity
242
MeshTopology& topology = mesh.topology();
243
MeshConnectivity& connectivity = topology(d0, d1);
246
dolfin_assert(d0 >= d1);
248
// Need connectivity d0 - d and d - d1
249
dolfin_assert(topology(d0, d).size() > 0);
250
dolfin_assert(topology(d, d1).size() > 0);
253
Array<uint> tmp(topology.size(d0));
255
// Reset size for each entity
256
for (uint i = 0; i < tmp.size(); i++)
259
// FIXME: Check efficiency of std::set, compare with std::vector
261
// A set with connected entities
262
std::set<uint> entities;
264
// Iterate over all entities of dimension d0
265
for (MeshEntityIterator e0(mesh, d0); !e0.end(); ++e0)
267
// Clear set of connected entities
270
// Iterate over all connected entities of dimension d
271
for (MeshEntityIterator e(*e0, d); !e.end(); ++e)
273
// Iterate over all connected entities of dimension d1
274
for (MeshEntityIterator e1(*e, d1); !e1.end(); ++e1)
278
// An entity is not a neighbor to itself
279
if ( e0->index() != e1->index() )
280
entities.insert(e1->index());
284
// Entity e1 must be completely contained in e0
285
if ( contains(*e0, *e1) )
286
entities.insert(e1->index());
291
// Count the number of connected entities
292
tmp[e0->index()] = entities.size();
295
// Initialize the number of connections
296
connectivity.init(tmp);
298
// Iterate over all entities of dimension d
299
for (MeshEntityIterator e0(mesh, d0); !e0.end(); ++e0)
301
// Clear set of connected entities
304
// Iterate over all connected entities of dimension d
305
for (MeshEntityIterator e(*e0, d); !e.end(); ++e)
307
// Iterate over all connected entities of dimension d1
308
for (MeshEntityIterator e1(*e, d1); !e1.end(); ++e1)
312
// An entity is not a neighbor to itself
313
if ( e0->index() != e1->index() )
314
entities.insert(e1->index());
318
// Entity e1 must be completely contained in e0
319
if ( contains(*e0, *e1) )
320
entities.insert(e1->index());
325
// Add the connected entities
327
for (std::set<uint>::iterator it = entities.begin(); it != entities.end(); ++it)
328
connectivity.set(e0->index(), *it, pos++);
331
//----------------------------------------------------------------------------
332
dolfin::uint TopologyComputation::countEntities(Mesh& mesh, MeshEntity& cell,
333
uint** entities, uint m, uint n,
336
// For each entity, we iterate over connected and previously visited
337
// cells to see if the entity has already been counted.
339
// Needs to be a cell
340
dolfin_assert(cell.dim() == mesh.topology().dim());
342
// Iterate over the given list of entities
343
uint num_entities = 0;
344
for (uint i = 0; i < m; i++)
346
// Iterate over connected cells and look for entity
347
for (MeshEntityIterator c(cell, mesh.topology().dim()); !c.end(); ++c)
349
// Check only previously visited cells
350
if ( c->index() >= cell.index() )
353
// Check for vertices
354
if ( contains(c->entities(0), c->numEntities(0), entities[i], n) )
361
// Entity found, don't need to count
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)
375
// We repeat the same algorithm as in countEntities() but this time
376
// we add any entities that are new.
378
// Needs to be a cell
379
dolfin_assert(cell.dim() == mesh.topology().dim());
381
// Iterate over the given list of entities
382
for (uint i = 0; i < m; i++)
384
// Iterate over connected cells and look for entity
385
for (MeshEntityIterator c(cell, mesh.topology().dim()); !c.end(); ++c)
387
// Check only previously visited cells
388
if ( c->index() >= cell.index() )
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++)
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) )
400
// Entity already exists, so pick the index
401
ce.set(cell.index(), e.index(), i);
407
// Entity does not exist, so create it
408
ce.set(cell.index(), current_entity, i);
409
ev.set(current_entity, entities[i]);
414
// Entity found, don't need to create
419
//----------------------------------------------------------------------------
420
bool TopologyComputation::contains(MeshEntity& e0, MeshEntity& e1)
423
return contains(e0.entities(0), e0.numEntities(0),
424
e1.entities(0), e1.numEntities(0));
426
//----------------------------------------------------------------------------
427
bool TopologyComputation::contains(uint* v0, uint n0, uint* v1, uint n1)
432
for (uint i1 = 0; i1 < n1; i1++)
435
for (uint i0 = 0; i0 < n0; i0++)
437
if ( v0[i0] == v1[i1] )
449
//----------------------------------------------------------------------------