~njansson/dolfin/hpc

« back to all changes in this revision

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

  • Committer: Anders Logg
  • Date: 2007-01-10 09:04:44 UTC
  • mfrom: (1689.1.221 trunk)
  • Revision ID: logg@simula.no-20070110090444-ecyux3n1qnei4i8h
RemoveĀ oldĀ head

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
// Copyright (C) 2003-2006 Johan Hoffman and Anders Logg.
2
 
// Licensed under the GNU GPL Version 2.
3
 
//
4
 
// First added:  2003
5
 
// Last changed: 2006-02-20
6
 
 
7
 
#include <dolfin/dolfin_log.h>
8
 
#include <dolfin/Mesh.h>
9
 
#include <dolfin/BoundaryInit.h>
10
 
 
11
 
using namespace dolfin;
12
 
 
13
 
//-----------------------------------------------------------------------------
14
 
void BoundaryInit::init(Mesh& mesh)
15
 
{
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.
19
 
 
20
 
  // Write a message
21
 
  dolfin_begin("Computing boundary:");
22
 
 
23
 
  clear(mesh);
24
 
 
25
 
  initFaces(mesh);
26
 
  initEdges(mesh);
27
 
  initVertices(mesh);
28
 
  
29
 
  dolfin_end();
30
 
}
31
 
//-----------------------------------------------------------------------------
32
 
void BoundaryInit::clear(Mesh& mesh)
33
 
{
34
 
  mesh.bd->clear();
35
 
}
36
 
//----------------------------------------------------------------------------- 
37
 
void BoundaryInit::initFaces(Mesh& mesh)
38
 
{
39
 
  switch ( mesh.type() ) {
40
 
  case Mesh::triangles:
41
 
    initFacesTri(mesh);
42
 
    break;
43
 
  case Mesh::tetrahedra:
44
 
    initFacesTet(mesh);
45
 
    break;
46
 
  default:
47
 
    dolfin_error("Unknown cell type.");
48
 
  }
49
 
}
50
 
//-----------------------------------------------------------------------------
51
 
void BoundaryInit::initEdges(Mesh& mesh)
52
 
{
53
 
  switch ( mesh.type() ) {
54
 
  case Mesh::triangles:
55
 
    initEdgesTri(mesh);
56
 
    break;
57
 
  case Mesh::tetrahedra:
58
 
    initEdgesTet(mesh);
59
 
    break;
60
 
  default:
61
 
    dolfin_error("Unknown cell type.");
62
 
  }
63
 
}
64
 
//-----------------------------------------------------------------------------
65
 
void BoundaryInit::initVertices(Mesh& mesh)
66
 
{
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.
70
 
 
71
 
  PArray<bool> marker(mesh.numVertices());
72
 
  marker = false;
73
 
 
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;
78
 
  }
79
 
 
80
 
  // Add all vertices on the boundary
81
 
  for (VertexIterator n(mesh); !n.end(); ++n)
82
 
    if ( marker(n->id()) )
83
 
      mesh.bd->add(*n);
84
 
  
85
 
  // Write a message
86
 
  cout << "Found " << mesh.bd->noVertices() << " vertices on the boundary." << endl;
87
 
}
88
 
//-----------------------------------------------------------------------------
89
 
void BoundaryInit::initFacesTri(Mesh& mesh)
90
 
{
91
 
  // Do nothing
92
 
}
93
 
//-----------------------------------------------------------------------------
94
 
void BoundaryInit::initFacesTet(Mesh& mesh)
95
 
{
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
100
 
  // removed
101
 
 
102
 
  PArray<int> cellcount(mesh.numFaces());
103
 
  cellcount = 0;
104
 
 
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;
109
 
  
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 )
113
 
      mesh.bd->add(*f);
114
 
    else if ( cellcount(f->id()) != 2 )
115
 
      dolfin_error1("Inconsistent mesh. Found face with %d cell neighbors.", cellcount(f->id()));
116
 
  }
117
 
  
118
 
  // Check that we found a boundary
119
 
  if ( mesh.bd->noFaces() == 0 )
120
 
    dolfin_error("Found no faces on the boundary.");
121
 
  
122
 
  // Write a message
123
 
  cout << "Found " << mesh.bd->noFaces() << " faces on the boundary." << endl;
124
 
}
125
 
//-----------------------------------------------------------------------------
126
 
void BoundaryInit::initEdgesTri(Mesh& mesh)
127
 
{
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
133
 
  // have been removed
134
 
 
135
 
  PArray<int> cellcount(mesh.numEdges());
136
 
  cellcount = 0;
137
 
 
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;
142
 
  
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 )
146
 
      mesh.bd->add(*e);
147
 
    else if ( cellcount(e->id()) != 2 )
148
 
      dolfin_error1("Inconsistent mesh. Found edge with %d cell neighbors.", cellcount(e->id()));
149
 
  }
150
 
  
151
 
  // Check that we found a boundary
152
 
  if ( mesh.bd->noEdges() == 0 )
153
 
    dolfin_error("Found no edges on the boundary.");
154
 
  
155
 
  // Write a message
156
 
  cout << "Found " << mesh.bd->noEdges() << " edges on the boundary." << endl;
157
 
}
158
 
//-----------------------------------------------------------------------------
159
 
void BoundaryInit::initEdgesTet(Mesh& mesh)
160
 
{
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.
164
 
 
165
 
  PArray<bool> marker(mesh.numEdges());
166
 
  marker = false;
167
 
 
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;
172
 
 
173
 
  // Add all edges on the boundary
174
 
  for (EdgeIterator e(mesh); !e.end(); ++e)
175
 
    if ( marker(e->id()) )
176
 
      mesh.bd->add(*e);
177
 
 
178
 
  // Write a message
179
 
  cout << "Found " << mesh.bd->noEdges() << " edges on the boundary." << endl;
180
 
}
181
 
//-----------------------------------------------------------------------------