fixed writing out entries in advective bc
[OpenFOAM-1.6-ext.git] / src / edgeMesh / edgeMesh.C
blob1eae3fba173205ba65938223f4b429f528842f77
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright held by original author
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
11     OpenFOAM is free software; you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by the
13     Free Software Foundation; either version 2 of the License, or (at your
14     option) any later version.
16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19     for more details.
21     You should have received a copy of the GNU General Public License
22     along with OpenFOAM; if not, write to the Free Software Foundation,
23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "edgeMesh.H"
28 #include "mergePoints.H"
29 #include "StaticHashTable.H"
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
35 void Foam::edgeMesh::calcPointEdges() const
37     if (pointEdgesPtr_.valid())
38     {
39         FatalErrorIn("edgeMesh::calcPointEdges() const")
40             << "pointEdges already calculated." << abort(FatalError);
41     }
43     pointEdgesPtr_.reset(new labelListList(points_.size()));
44     labelListList& pointEdges = pointEdgesPtr_();
46     // Count
47     labelList nEdgesPerPoint(points_.size(), 0);
49     forAll(edges_, edgeI)
50     {
51         const edge& e = edges_[edgeI];
53         nEdgesPerPoint[e[0]]++;
54         nEdgesPerPoint[e[1]]++;
55     }
57     // Size
58     forAll(pointEdges, pointI)
59     {
60         pointEdges[pointI].setSize(nEdgesPerPoint[pointI]);
61     }
63     // Fill
64     nEdgesPerPoint = 0;
66     forAll(edges_, edgeI)
67     {
68         const edge& e = edges_[edgeI];
70         label p0 = e[0];
71         pointEdges[p0][nEdgesPerPoint[p0]++] = edgeI;
72         label p1 = e[1];
73         pointEdges[p1][nEdgesPerPoint[p1]++] = edgeI;
74     }
78 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
80 // construct from components
81 Foam::edgeMesh::edgeMesh(const pointField& points, const edgeList& edges)
83     points_(points),
84     edges_(edges)
88 // construct as copy
89 Foam::edgeMesh::edgeMesh(const edgeMesh& em)
91     points_(em.points_),
92     edges_(em.edges_),
93     pointEdgesPtr_(NULL)
97 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
99 Foam::label Foam::edgeMesh::regions(labelList& edgeRegion) const
101     edgeRegion.setSize(edges_.size());
102     edgeRegion = -1;
104     label startEdgeI = 0;
106     label currentRegion = 0;
108     while (true)
109     {
110         while (startEdgeI < edges_.size() && edgeRegion[startEdgeI] != -1)
111         {
112             startEdgeI++;
113         }
115         if (startEdgeI == edges_.size())
116         {
117             break;
118         }
120         // Found edge that has not yet been assigned a region.
121         // Mark connected region with currentRegion starting at startEdgeI.
123         edgeRegion[startEdgeI] = currentRegion;
124         labelList edgesToVisit(1, startEdgeI);
126         while (edgesToVisit.size())
127         {
128             // neighbours of current edgesToVisit
129             DynamicList<label> newEdgesToVisit(edgesToVisit.size());
131             // Mark all point connected edges with current region.
132             forAll(edgesToVisit, i)
133             {
134                 label edgeI = edgesToVisit[i];
136                 // Mark connected edges
137                 const edge& e = edges_[edgeI];
139                 forAll(e, fp)
140                 {
141                     const labelList& pEdges = pointEdges()[e[fp]];
143                     forAll(pEdges, pEdgeI)
144                     {
145                         label nbrEdgeI = pEdges[pEdgeI];
147                         if (edgeRegion[nbrEdgeI] == -1)
148                         {
149                             edgeRegion[nbrEdgeI] = currentRegion;
150                             newEdgesToVisit.append(nbrEdgeI);
151                         }
152                     }
153                 }
154             }
156             edgesToVisit.transfer(newEdgesToVisit);
157         }
159         currentRegion++;
160     }
161     return currentRegion;
165 void Foam::edgeMesh::mergePoints(const scalar mergeDist)
167     pointField newPoints;
168     labelList pointMap;
170     bool hasMerged = Foam::mergePoints
171     (
172         points_,
173         mergeDist,
174         false,
175         pointMap,
176         newPoints,
177         vector::zero
178     );
180     if (hasMerged)
181     {
182         pointEdgesPtr_.clear();
184         points_.transfer(newPoints);
186         // Renumber and make sure e[0] < e[1] (not really nessecary)
187         forAll(edges_, edgeI)
188         {
189             edge& e = edges_[edgeI];
191             label p0 = pointMap[e[0]];
192             label p1 = pointMap[e[1]];
194             if (p0 < p1)
195             {
196                 e[0] = p0;
197                 e[1] = p1;
198             }
199             else
200             {
201                 e[0] = p1;
202                 e[1] = p0;
203             }
204         }
206         // Compact using a hashtable and commutative hash of edge.
207         StaticHashTable<label, edge, Hash<edge> > edgeToLabel
208         (
209             2*edges_.size()
210         );
212         label newEdgeI = 0;
214         forAll(edges_, edgeI)
215         {
216             const edge& e = edges_[edgeI];
218             if (e[0] != e[1])
219             {
220                 if (edgeToLabel.insert(e, newEdgeI))
221                 {
222                     newEdgeI++;
223                 }
224             }
225         }
227         edges_.setSize(newEdgeI);
229         for
230         (
231             StaticHashTable<label, edge, Hash<edge> >::const_iterator iter =
232                 edgeToLabel.begin();
233             iter != edgeToLabel.end();
234             ++iter
235         )
236         {
237             edges_[iter()] = iter.key();
238         }
239     }
243 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
245 void Foam::edgeMesh::operator=(const edgeMesh& rhs)
247     points_ = rhs.points_;
248     edges_ = rhs.edges_;
249     pointEdgesPtr_.reset(NULL);
253 // ************************************************************************* //