1
// Copyright (C) 2003-2006 Johan Hoffman and Anders Logg.
2
// Licensed under the GNU GPL Version 2.
5
// Last changed: 2006-02-20
7
#include <dolfin/dolfin_log.h>
8
#include <dolfin/Mesh.h>
9
#include <dolfin/BoundaryInit.h>
11
using namespace dolfin;
13
//-----------------------------------------------------------------------------
14
void BoundaryInit::init(Mesh& mesh)
16
// It is important that the computation of boundary data is done in
17
// the correct order: First compute faces, then use this information
18
// to compute edges, and finally compute the vertices from the edges.
21
dolfin_begin("Computing boundary:");
31
//-----------------------------------------------------------------------------
32
void BoundaryInit::clear(Mesh& mesh)
36
//-----------------------------------------------------------------------------
37
void BoundaryInit::initFaces(Mesh& mesh)
39
switch ( mesh.type() ) {
43
case Mesh::tetrahedra:
47
dolfin_error("Unknown cell type.");
50
//-----------------------------------------------------------------------------
51
void BoundaryInit::initEdges(Mesh& mesh)
53
switch ( mesh.type() ) {
57
case Mesh::tetrahedra:
61
dolfin_error("Unknown cell type.");
64
//-----------------------------------------------------------------------------
65
void BoundaryInit::initVertices(Mesh& mesh)
67
// Compute vertices from edges. A vertex is on the boundary if it belongs to
68
// an edge which is on the boundary. We need an extra list so that we don't
69
// add a vertex more than once.
71
PArray<bool> marker(mesh.numVertices());
74
// Mark all vertices which are on the boundary
75
for (PList<Edge*>::Iterator e(mesh.bd->edges); !e.end(); ++e) {
76
marker((*e)->vertex(0).id()) = true;
77
marker((*e)->vertex(1).id()) = true;
80
// Add all vertices on the boundary
81
for (VertexIterator n(mesh); !n.end(); ++n)
82
if ( marker(n->id()) )
86
cout << "Found " << mesh.bd->noVertices() << " vertices on the boundary." << endl;
88
//-----------------------------------------------------------------------------
89
void BoundaryInit::initFacesTri(Mesh& mesh)
93
//-----------------------------------------------------------------------------
94
void BoundaryInit::initFacesTet(Mesh& mesh)
96
// Go through all faces and for each face check if it is on the
97
// boundary. A face is on the boundary if it is contained in only
98
// one cell. A list is used to count the number of cell neighbors
99
// for all faces. Warning: may not work if some faces have been
102
PArray<int> cellcount(mesh.numFaces());
105
// Count the number of cell neighbors for each face
106
for (CellIterator c(mesh); !c.end(); ++c)
107
for (FaceIterator f(c); !f.end(); ++f)
108
cellcount(f->id()) += 1;
110
// Add faces with only one cell neighbor to the boundary
111
for (FaceIterator f(mesh); !f.end(); ++f) {
112
if ( cellcount(f->id()) == 1 )
114
else if ( cellcount(f->id()) != 2 )
115
dolfin_error1("Inconsistent mesh. Found face with %d cell neighbors.", cellcount(f->id()));
118
// Check that we found a boundary
119
if ( mesh.bd->noFaces() == 0 )
120
dolfin_error("Found no faces on the boundary.");
123
cout << "Found " << mesh.bd->noFaces() << " faces on the boundary." << endl;
125
//-----------------------------------------------------------------------------
126
void BoundaryInit::initEdgesTri(Mesh& mesh)
128
// Go through all edges and for each edge check if it is on the
129
// boundary. This is similar to what is done in initFacesTet() for
130
// tetrahedra. An edge is on the boundary if it is contained in
131
// only one cell. A list is used to count the number of cell
132
// neighbors for all edges. Warning: may not work if some edges
135
PArray<int> cellcount(mesh.numEdges());
138
// Count the number of cell neighbors for each edge
139
for (CellIterator c(mesh); !c.end(); ++c)
140
for (EdgeIterator e(c); !e.end(); ++e)
141
cellcount(e->id()) += 1;
143
// Add edges with only one cell neighbor to the boundary
144
for (EdgeIterator e(mesh); !e.end(); ++e) {
145
if ( cellcount(e->id()) == 1 )
147
else if ( cellcount(e->id()) != 2 )
148
dolfin_error1("Inconsistent mesh. Found edge with %d cell neighbors.", cellcount(e->id()));
151
// Check that we found a boundary
152
if ( mesh.bd->noEdges() == 0 )
153
dolfin_error("Found no edges on the boundary.");
156
cout << "Found " << mesh.bd->noEdges() << " edges on the boundary." << endl;
158
//-----------------------------------------------------------------------------
159
void BoundaryInit::initEdgesTet(Mesh& mesh)
161
// Compute edges from faces. An edge is on the boundary if it belongs to
162
// a face which is on the boundary. We need an extra list so that we don't
163
// add an edge more than once.
165
PArray<bool> marker(mesh.numEdges());
168
// Mark all edges which are on the boundary
169
for (PList<Face*>::Iterator f(mesh.bd->faces); !f.end(); ++f)
170
for (EdgeIterator e(**f); !e.end(); ++e)
171
marker(e->id()) = true;
173
// Add all edges on the boundary
174
for (EdgeIterator e(mesh); !e.end(); ++e)
175
if ( marker(e->id()) )
179
cout << "Found " << mesh.bd->noEdges() << " edges on the boundary." << endl;
181
//-----------------------------------------------------------------------------