1
// Copyright (C) 2003-2006 Johan Hoffman and Anders Logg.
2
// Licensed under the GNU GPL Version 2.
4
// Modified by Par Ingelstrom 2004.
7
// Last changed: 2006-02-20
9
#include <dolfin/dolfin_log.h>
10
#include <dolfin/Mesh.h>
11
#include <dolfin/Cell.h>
12
#include <dolfin/Vertex.h>
13
#include <dolfin/TetMeshRefinement.h>
15
using namespace dolfin;
17
//-----------------------------------------------------------------------------
18
bool TetMeshRefinement::checkRule(Cell& cell, int no_marked_edges)
20
dolfin_assert(cell.type() == Cell::tetrahedron);
22
// Choose refinement rule
24
if ( checkRuleRegular(cell, no_marked_edges) )
27
if ( checkRuleIrregular1(cell, no_marked_edges) )
30
if ( checkRuleIrregular2(cell, no_marked_edges) )
33
if ( checkRuleIrregular3(cell, no_marked_edges) )
36
if ( checkRuleIrregular4(cell, no_marked_edges) )
39
// We didn't find a matching rule for refinement
42
//-----------------------------------------------------------------------------
43
void TetMeshRefinement::refine(Cell& cell, Mesh& mesh)
45
// Refine cell according to marker
46
switch ( cell.marker() ) {
47
case Cell::marked_for_no_ref:
48
refineNoRefine(cell, mesh);
50
case Cell::marked_for_reg_ref:
51
refineRegular(cell, mesh);
53
case Cell::marked_for_irr_ref_1:
54
refineIrregular1(cell, mesh);
56
case Cell::marked_for_irr_ref_2:
57
refineIrregular2(cell, mesh);
59
case Cell::marked_for_irr_ref_3:
60
refineIrregular3(cell, mesh);
62
case Cell::marked_for_irr_ref_4:
63
refineIrregular4(cell, mesh);
66
// We should not reach this case, cell cannot be
67
// marked_for_coarsening or marked_according_to_ref
68
dolfin_error("Inconsistent cell markers.");
71
//-----------------------------------------------------------------------------
72
bool TetMeshRefinement::checkRuleRegular(Cell& cell, int no_marked_edges)
74
// Check if cell should be regularly refined.
75
// A cell is refined regularly if all edges are marked.
77
if ( no_marked_edges != 6 )
80
cell.marker() = Cell::marked_for_reg_ref;
83
//-----------------------------------------------------------------------------
84
bool TetMeshRefinement::checkRuleIrregular1(Cell& cell, int no_marked_edges)
86
// Check if cell matches irregular refinement rule 1
88
if ( no_marked_edges != 3 )
91
if ( !markedEdgesOnSameFace(cell) )
94
cell.marker() = Cell::marked_for_irr_ref_1;
97
//-----------------------------------------------------------------------------
98
bool TetMeshRefinement::checkRuleIrregular2(Cell& cell, int no_marked_edges)
100
// Check if cell matches irregular refinement rule 2
102
if ( no_marked_edges != 1 )
105
cell.marker() = Cell::marked_for_irr_ref_2;
108
//-----------------------------------------------------------------------------
109
bool TetMeshRefinement::checkRuleIrregular3(Cell& cell, int no_marked_edges)
111
// Check if cell matches irregular refinement rule 3
113
if ( no_marked_edges != 2 )
116
if ( !markedEdgesOnSameFace(cell) )
119
cell.marker() = Cell::marked_for_irr_ref_3;
122
//-----------------------------------------------------------------------------
123
bool TetMeshRefinement::checkRuleIrregular4(Cell& cell, int no_marked_edges)
125
// Check if cell matches irregular refinement rule 4
127
if ( no_marked_edges != 3 )
130
// Note that this has already been checked by checkRule3(), but this
131
// way the algorithm is a little cleaner.
132
if ( markedEdgesOnSameFace(cell) )
135
cell.marker() = Cell::marked_for_irr_ref_4;
138
//-----------------------------------------------------------------------------
139
void TetMeshRefinement::refineNoRefine(Cell& cell, Mesh& mesh)
141
// Don't refine the tetrahedron and create a copy in the new mesh.
143
// Check that the cell is marked correctly
144
dolfin_assert(cell.marker() == Cell::marked_for_no_ref);
146
// Create new vertices with the same coordinates as existing vertices
147
Vertex& n0 = createVertex(cell.vertex(0), mesh, cell);
148
Vertex& n1 = createVertex(cell.vertex(1), mesh, cell);
149
Vertex& n2 = createVertex(cell.vertex(2), mesh, cell);
150
Vertex& n3 = createVertex(cell.vertex(3), mesh, cell);
153
cell.initChildren(1);
154
createCell(n0, n1, n2, n3, mesh, cell);
156
// Set marker of cell
157
cell.marker() = Cell::marked_according_to_ref;
159
// Set status of cell
160
cell.status() = Cell::unref;
162
//-----------------------------------------------------------------------------
163
void TetMeshRefinement::refineRegular(Cell& cell, Mesh& mesh)
165
// Refine 1 tetrahedron into 8 new ones, introducing new vertices
166
// at the midpoints of the edges.
168
// Check that cell's parent is not refined irregularly,
169
// since then it should be further refined
170
dolfin_assert(okToRefine(cell));
172
// Check that the cell is marked correctly
173
dolfin_assert(cell.marker() == Cell::marked_for_reg_ref);
175
// Create new vertices with the same coordinates as existing vertices
176
Vertex& n0 = createVertex(cell.vertex(0), mesh, cell);
177
Vertex& n1 = createVertex(cell.vertex(1), mesh, cell);
178
Vertex& n2 = createVertex(cell.vertex(2), mesh, cell);
179
Vertex& n3 = createVertex(cell.vertex(3), mesh, cell);
181
// Create new vertices with the new coordinates
182
Vertex& n01 = createVertex(cell.vertex(0).midpoint(cell.vertex(1)), mesh, cell);
183
Vertex& n02 = createVertex(cell.vertex(0).midpoint(cell.vertex(2)), mesh, cell);
184
Vertex& n03 = createVertex(cell.vertex(0).midpoint(cell.vertex(3)), mesh, cell);
185
Vertex& n12 = createVertex(cell.vertex(1).midpoint(cell.vertex(2)), mesh, cell);
186
Vertex& n13 = createVertex(cell.vertex(1).midpoint(cell.vertex(3)), mesh, cell);
187
Vertex& n23 = createVertex(cell.vertex(2).midpoint(cell.vertex(3)), mesh, cell);
190
cell.initChildren(8);
191
createCell(n0, n01, n02, n03, mesh, cell);
192
createCell(n01, n1, n12, n13, mesh, cell);
193
createCell(n02, n12, n2, n23, mesh, cell);
194
createCell(n03, n13, n23, n3, mesh, cell);
195
createCell(n01, n02, n03, n13, mesh, cell);
196
createCell(n01, n02, n12, n13, mesh, cell);
197
createCell(n02, n03, n13, n23, mesh, cell);
198
createCell(n02, n12, n13, n23, mesh, cell);
200
// Set marker of cell
201
cell.marker() = Cell::marked_according_to_ref;
203
// Set status of cell
204
cell.status() = Cell::ref_reg;
206
//-----------------------------------------------------------------------------
207
void TetMeshRefinement::refineIrregular1(Cell& cell, Mesh& mesh)
209
// Three edges are marked on the same face. Insert three new vertices
210
// at the midpoints on the marked edges, connect the new vertices to
211
// each other, as well as to the vertex that is not on the marked
212
// face. This gives 4 new tetrahedra.
214
// Check that cell's parent is not refined irregularly,
215
// since then it should be further refined
216
dolfin_assert(okToRefine(cell));
218
// Check that the cell is marked correctly
219
dolfin_assert(cell.marker() == Cell::marked_for_irr_ref_1);
221
// Sort vertices by the number of marked edges
222
PArray<Vertex*> vertices;
223
sortVertices(cell,vertices);
225
// Create new vertices with the same coordinates as the old vertices
226
Vertex& n0 = createVertex(*vertices(0), mesh, cell);
227
Vertex& n1 = createVertex(*vertices(1), mesh, cell);
228
Vertex& n2 = createVertex(*vertices(2), mesh, cell);
229
Vertex& nn = createVertex(*vertices(3), mesh, cell); // Not marked
232
Edge* e01 = cell.findEdge(*vertices(0), *vertices(1));
233
Edge* e02 = cell.findEdge(*vertices(0), *vertices(2));
234
Edge* e12 = cell.findEdge(*vertices(1), *vertices(2));
239
// Create new vertices on the edges of the marked face
240
Vertex& n01 = createVertex(e01->midpoint(), mesh, cell);
241
Vertex& n02 = createVertex(e02->midpoint(), mesh, cell);
242
Vertex& n12 = createVertex(e12->midpoint(), mesh, cell);
245
cell.initChildren(4);
246
createCell(nn, n01, n02, n12, mesh, cell);
247
createCell(nn, n01, n02, n0, mesh, cell);
248
createCell(nn, n01, n12, n1, mesh, cell);
249
createCell(nn, n02, n12, n2, mesh, cell);
251
// Set marker of cell
252
cell.marker() = Cell::marked_according_to_ref;
254
// Set status of cell
255
cell.status() = Cell::ref_irr;
257
//-----------------------------------------------------------------------------
258
void TetMeshRefinement::refineIrregular2(Cell& cell, Mesh& mesh)
260
// One edge is marked. Insert one new vertex at the midpoint of the
261
// marked edge, then connect this new vertex to the two vertices not on
262
// the marked edge. This gives 2 new tetrahedra.
264
// Check that cell's parent is not refined irregularly,
265
// since then it should be further refined
266
dolfin_assert(okToRefine(cell));
268
// Check that the cell is marked correctly
269
dolfin_assert(cell.marker() == Cell::marked_for_irr_ref_2);
271
// Sort vertices by the number of marked edges
272
PArray<Vertex*> vertices;
273
sortVertices(cell, vertices);
275
// Create new vertices with the same coordinates as the old vertices
276
Vertex& n0 = createVertex(*vertices(0), mesh, cell);
277
Vertex& n1 = createVertex(*vertices(1), mesh, cell);
278
Vertex& nn0 = createVertex(*vertices(2), mesh, cell); // Not marked
279
Vertex& nn1 = createVertex(*vertices(3), mesh, cell); // Not marked
281
// Find the marked edge
282
Edge* e = cell.findEdge(*vertices(0), *vertices(1));
285
// Create new vertex on marked edge
286
Vertex& ne = createVertex(e->midpoint(), mesh, cell);
289
cell.initChildren(2);
290
createCell(ne, nn0, nn1, n0, mesh, cell);
291
createCell(ne, nn0, nn1, n1, mesh, cell);
293
// Set marker of cell
294
cell.marker() = Cell::marked_according_to_ref;
296
// Set status of cell
297
cell.status() = Cell::ref_irr;
299
//-----------------------------------------------------------------------------
300
void TetMeshRefinement::refineIrregular3(Cell& cell, Mesh& mesh)
302
// Two edges are marked, both on the same face. There are two
303
// possibilities, and the chosen alternative must match the
304
// corresponding face of the neighbor tetrahedron. If this neighbor
305
// is marked for regular refinement, so is this tetrahedron.
307
// We insert two new vertices at the midpoints of the marked edges.
308
// Three new edges are created by connecting the two new vertices to
309
// each other and to the vertex opposite the face of the two marked
310
// edges. Finally, an edge is created by either
312
// (1) connecting new vertex 1 with the endvertex of marked edge 2,
313
// that is not common with marked edge 1, or
315
// (2) connecting new vertex 2 with the endvertex of marked edge 1,
316
// that is not common with marked edge 2.
318
// Check that cell's parent is not refined irregularly,
319
// since then it should be further refined
320
dolfin_assert(okToRefine(cell));
322
// Check that the cell is marked correctly
323
dolfin_assert(cell.marker() == Cell::marked_for_irr_ref_3);
325
// Sort vertices by the number of marked edges
326
PArray<Vertex*> vertices;
327
sortVertices(cell, vertices);
330
Edge* e0 = cell.findEdge(*vertices(0), *vertices(1));
331
Edge* e1 = cell.findEdge(*vertices(0), *vertices(2));
336
Face* f = cell.findFace(*e0, *e1);
339
// Find neighbor with common face
340
Cell* neighbor = findNeighbor(cell, *f);
341
dolfin_assert(neighbor);
343
if ( neighbor->marker() == Cell::marked_for_reg_ref ) {
344
// If neighbor is marked for regular refinement so is cell
345
refineIrregular31(cell,mesh);
347
else if ( neighbor->marker() == Cell::marked_for_irr_ref_3) {
348
// If neighbor is marked refinement by rule 3,
349
// just chose an orientation, and it will be up to
350
// the neighbor to make sure the common face match
351
refineIrregular32(cell,mesh,vertices);
353
else if ( neighbor->marker() == Cell::marked_according_to_ref ) {
354
// If neighbor has been refined irregular according to
355
// refinement rule 3, make sure the common face matches
356
refineIrregular33(cell, mesh, vertices, *neighbor);
359
// This case shouldn't happen
360
dolfin_error("Unable to handle refinement rule of neighbor.");
363
//-----------------------------------------------------------------------------
364
void TetMeshRefinement::refineIrregular4(Cell& cell, Mesh& mesh)
366
// Two edges are marked, opposite to each other. We insert two new
367
// vertices at the midpoints of the marked edges, insert a new edge
368
// between the two vertices, and insert four new edges by connecting
369
// the new vertices to the endpoints of the opposite edges.
371
// Check that cell's parent is not refined irregularly,
372
// since then it should be further refined
373
dolfin_assert(okToRefine(cell));
375
// Check that the cell is marked correctly
376
dolfin_assert(cell.marker() == Cell::marked_for_irr_ref_4);
378
// Find the two marked edges
379
PArray<Edge*> marked_edges(2);
382
for (EdgeIterator e(cell); !e.end(); ++e)
383
if (e->marked()) marked_edges(cnt++) = e;
385
// Create new vertices with the same coordinates as the old vertices
386
Vertex& n00 = createVertex(marked_edges(0)->vertex(0), mesh, cell);
387
Vertex& n01 = createVertex(marked_edges(0)->vertex(1), mesh, cell);
388
Vertex& n10 = createVertex(marked_edges(1)->vertex(0), mesh, cell);
389
Vertex& n11 = createVertex(marked_edges(1)->vertex(1), mesh, cell);
391
// Create new vertex on marked edge
392
Vertex& n_e0 = createVertex(marked_edges(0)->midpoint(), mesh, cell);
393
Vertex& n_e1 = createVertex(marked_edges(1)->midpoint(), mesh, cell);
396
cell.initChildren(4);
397
createCell(n_e0, n_e1, n00, n10, mesh, cell);
398
createCell(n_e0, n_e1, n00, n11, mesh, cell);
399
createCell(n_e0, n_e1, n01, n10, mesh, cell);
400
createCell(n_e0, n_e1, n01, n11, mesh, cell);
402
// Set marker of cell
403
cell.marker() = Cell::marked_according_to_ref;
405
// Set status of cell
406
cell.status() = Cell::ref_irr;
408
//-----------------------------------------------------------------------------
409
void TetMeshRefinement::refineIrregular31(Cell& cell, Mesh& mesh)
411
// If neighbor with common 2-marked-edges-face is marked for regular
412
// refinement do the same for cell
413
cell.marker() = Cell::marked_for_reg_ref;
414
refineRegular(cell,mesh);
416
//-----------------------------------------------------------------------------
417
void TetMeshRefinement::refineIrregular32(Cell& cell, Mesh& mesh,
418
PArray<Vertex*>& sorted_vertices)
420
// If neighbor is marked refinement by rule 3,
421
// just chose an orientation, and it will be up to
422
// the neighbor to make sure the common face match
424
// Create new vertices with the same coordinates as the old vertices
425
Vertex& n_dm = createVertex(*sorted_vertices(0), mesh, cell);
426
Vertex& n_m0 = createVertex(*sorted_vertices(1), mesh, cell);
427
Vertex& n_m1 = createVertex(*sorted_vertices(2), mesh, cell);
428
Vertex& n_nm = createVertex(*sorted_vertices(3), mesh, cell);
430
// Create new vertex on marked edge
431
Vertex& n_e0 = createVertex(sorted_vertices(0)->midpoint(*sorted_vertices(1)), mesh, cell);
432
Vertex& n_e1 = createVertex(sorted_vertices(0)->midpoint(*sorted_vertices(2)), mesh, cell);
435
cell.initChildren(3);
436
createCell(n_dm, n_e0, n_e1, n_nm, mesh, cell);
437
createCell(n_m0, n_m1, n_e1, n_nm, mesh, cell);
438
createCell(n_e0, n_e1, n_m0, n_nm, mesh, cell);
440
// Set marker of cell
441
cell.marker() = Cell::marked_according_to_ref;
443
// Set status of cell
444
cell.status() = Cell::ref_irr;
446
//-----------------------------------------------------------------------------
447
void TetMeshRefinement::refineIrregular33(Cell& cell, Mesh& mesh,
448
PArray<Vertex*>& sorted_vertices,
452
// Create new vertices with the same coordinates as the old vertices
453
Vertex& n_dm = createVertex(*sorted_vertices(0), mesh, cell);
454
Vertex& n_m0 = createVertex(*sorted_vertices(1), mesh, cell);
455
Vertex& n_m1 = createVertex(*sorted_vertices(2), mesh, cell);
456
Vertex& n_nm = createVertex(*sorted_vertices(3), mesh, cell);
458
// Create new vertex on marked edge
459
Vertex& n_e0 = createVertex(sorted_vertices(0)->midpoint(*sorted_vertices(1)), mesh, cell);
460
Vertex& n_e1 = createVertex(sorted_vertices(0)->midpoint(*sorted_vertices(2)), mesh, cell);
463
cell.initChildren(3);
464
createCell(n_dm, n_e0, n_e1, n_nm, mesh, cell);
466
// If neighbor has been refined irregular according to
467
// refinement rule 3, make sure the common face matches
469
for (int i = 0; i < face_neighbor.numChildren(); i++) {
470
c = face_neighbor.child(i);
471
if ( !(c->haveVertex(n_dm)) ){
472
if ( c->haveVertex(n_e0) && c->haveVertex(n_e1) ){
473
if ( c->haveVertex(n_m0) ){
474
createCell(n_e0, n_e1, n_m0, n_nm, mesh, cell);
475
createCell(n_m0, n_m1, n_e1, n_nm, mesh, cell);
479
createCell(n_e0, n_e1, n_m1, n_nm, mesh, cell);
480
createCell(n_m0, n_m1, n_e0, n_nm, mesh, cell);
487
// Set marker of cell
488
cell.marker() = Cell::marked_according_to_ref;
490
// Set status of cell
491
cell.status() = Cell::ref_irr;
493
//-----------------------------------------------------------------------------
494
bool TetMeshRefinement::markedEdgesOnSameFace(Cell& cell)
496
// Check if the marked edges of cell are on the same face:
498
// 0 marked edge -> false
499
// 1 marked edge -> true
500
// 2 marked edges -> true if edges have any common vertices
501
// 3 marked edges -> true if there is a face with the marked edges
502
// 4 marked edges -> false
503
// 5 marked edges -> false
504
// 6 marked edges -> false
506
// Count the number of marked edges
508
for (EdgeIterator e(cell); !e.end(); ++e)
509
if (e->marked()) cnt++;
511
// Case 0, 1, 4, 5, 6
512
dolfin_assert(cnt >= 0 && cnt <= 6);
513
if (cnt == 0) return false;
514
if (cnt == 1) return true;
515
if (cnt > 3) return false;
517
// Create a list of the marked edges
518
PArray<Edge*> marked_edges(cnt);
521
for (EdgeIterator e(cell); !e.end(); ++e)
522
if (e->marked()) marked_edges(cnt++) = e;
524
// Check that number of marked edges are consistent
525
dolfin_assert(cnt == 2 || cnt == 3);
529
if (marked_edges(0)->contains(marked_edges(1)->vertex(0)) ||
530
marked_edges(0)->contains(marked_edges(1)->vertex(1)))
537
for (FaceIterator f(cell); !f.end(); ++f){
538
if (f->equals(*marked_edges(0), *marked_edges(1), *marked_edges(2)))
544
// We shouldn't reach this case
545
dolfin_error("Inconsistent edge markers.");
548
//-----------------------------------------------------------------------------
549
Cell* TetMeshRefinement::findNeighbor(Cell& cell, Face& face)
551
// Find a cell neighbor sharing a common face
553
for (CellIterator c(cell); !c.end(); ++c) {
555
// Don't check the cell itself
556
if ( c->id() == cell.id() )
558
for (FaceIterator f(*c); !f.end(); ++f) {
559
if ( f->id() == face.id() ) {
565
// If no neighbor is found, return the cell itself
568
//-----------------------------------------------------------------------------
569
Cell& TetMeshRefinement::createCell(Vertex& n0, Vertex& n1, Vertex& n2, Vertex& n3,
570
Mesh& mesh, Cell& cell)
572
Cell& c = mesh.createCell(n0, n1, n2, n3);
578
//-----------------------------------------------------------------------------
579
Cell& TetMeshRefinement::createChildCopy(Cell& cell, Mesh& mesh)
581
return createCell(*cell.vertex(0).child(), *cell.vertex(1).child(),
582
*cell.vertex(2).child(), *cell.vertex(3).child(), mesh, cell);
584
//-----------------------------------------------------------------------------
added
removed
src/kernel/mesh/TetMeshRefinement.cpp
