src/edgeMesh/edgeMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "edgeMesh.H"
  27 #include "mergePoints.H"
  28 #include "StaticHashTable.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 void Foam::edgeMesh::calcPointEdges() const
  35 {
  36     if (pointEdgesPtr_.valid())
  37     {
  38         FatalErrorIn("edgeMesh::calcPointEdges() const")
  39             << "pointEdges already calculated." << abort(FatalError);
  40     }
  41
  42     pointEdgesPtr_.reset(new labelListList(points_.size()));
  43     labelListList& pointEdges = pointEdgesPtr_();
  44
  45     // Count
  46     labelList nEdgesPerPoint(points_.size(), 0);
  47
  48     forAll(edges_, edgeI)
  49     {
  50         const edge& e = edges_[edgeI];
  51
  52         nEdgesPerPoint[e[0]]++;
  53         nEdgesPerPoint[e[1]]++;
  54     }
  55
  56     // Size
  57     forAll(pointEdges, pointI)
  58     {
  59         pointEdges[pointI].setSize(nEdgesPerPoint[pointI]);
  60     }
  61
  62     // Fill
  63     nEdgesPerPoint = 0;
  64
  65     forAll(edges_, edgeI)
  66     {
  67         const edge& e = edges_[edgeI];
  68
  69         label p0 = e[0];
  70         pointEdges[p0][nEdgesPerPoint[p0]++] = edgeI;
  71         label p1 = e[1];
  72         pointEdges[p1][nEdgesPerPoint[p1]++] = edgeI;
  73     }
  74 }
  75
  76
  77 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  78
  79 // construct from components
  80 Foam::edgeMesh::edgeMesh(const pointField& points, const edgeList& edges)
  81 :
  82     points_(points),
  83     edges_(edges)
  84 {}
  85
  86
  87 // construct as copy
  88 Foam::edgeMesh::edgeMesh(const edgeMesh& em)
  89 :
  90     points_(em.points_),
  91     edges_(em.edges_),
  92     pointEdgesPtr_(NULL)
  93 {}
  94
  95
  96 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  97
  98 Foam::label Foam::edgeMesh::regions(labelList& edgeRegion) const
  99 {
 100     edgeRegion.setSize(edges_.size());
 101     edgeRegion = -1;
 102
 103     label startEdgeI = 0;
 104
 105     label currentRegion = 0;
 106
 107     while (true)
 108     {
 109         while (startEdgeI < edges_.size() && edgeRegion[startEdgeI] != -1)
 110         {
 111             startEdgeI++;
 112         }
 113
 114         if (startEdgeI == edges_.size())
 115         {
 116             break;
 117         }
 118
 119         // Found edge that has not yet been assigned a region.
 120         // Mark connected region with currentRegion starting at startEdgeI.
 121
 122         edgeRegion[startEdgeI] = currentRegion;
 123         labelList edgesToVisit(1, startEdgeI);
 124
 125         while (edgesToVisit.size())
 126         {
 127             // neighbours of current edgesToVisit
 128             DynamicList<label> newEdgesToVisit(edgesToVisit.size());
 129
 130             // Mark all point connected edges with current region.
 131             forAll(edgesToVisit, i)
 132             {
 133                 label edgeI = edgesToVisit[i];
 134
 135                 // Mark connected edges
 136                 const edge& e = edges_[edgeI];
 137
 138                 forAll(e, fp)
 139                 {
 140                     const labelList& pEdges = pointEdges()[e[fp]];
 141
 142                     forAll(pEdges, pEdgeI)
 143                     {
 144                         label nbrEdgeI = pEdges[pEdgeI];
 145
 146                         if (edgeRegion[nbrEdgeI] == -1)
 147                         {
 148                             edgeRegion[nbrEdgeI] = currentRegion;
 149                             newEdgesToVisit.append(nbrEdgeI);
 150                         }
 151                     }
 152                 }
 153             }
 154
 155             edgesToVisit.transfer(newEdgesToVisit);
 156         }
 157
 158         currentRegion++;
 159     }
 160     return currentRegion;
 161 }
 162
 163
 164 void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 165 {
 166     pointField newPoints;
 167     labelList pointMap;
 168
 169     bool hasMerged = Foam::mergePoints
 170     (
 171         points_,
 172         mergeDist,
 173         false,
 174         pointMap,
 175         newPoints,
 176         vector::zero
 177     );
 178
 179     if (hasMerged)
 180     {
 181         pointEdgesPtr_.clear();
 182
 183         points_.transfer(newPoints);
 184
 185         // Renumber and make sure e[0] < e[1] (not really nessecary)
 186         forAll(edges_, edgeI)
 187         {
 188             edge& e = edges_[edgeI];
 189
 190             label p0 = pointMap[e[0]];
 191             label p1 = pointMap[e[1]];
 192
 193             if (p0 < p1)
 194             {
 195                 e[0] = p0;
 196                 e[1] = p1;
 197             }
 198             else
 199             {
 200                 e[0] = p1;
 201                 e[1] = p0;
 202             }
 203         }
 204
 205         // Compact using a hashtable and commutative hash of edge.
 206         StaticHashTable<label, edge, Hash<edge> > edgeToLabel
 207         (
 208             2*edges_.size()
 209         );
 210
 211         label newEdgeI = 0;
 212
 213         forAll(edges_, edgeI)
 214         {
 215             const edge& e = edges_[edgeI];
 216
 217             if (e[0] != e[1])
 218             {
 219                 if (edgeToLabel.insert(e, newEdgeI))
 220                 {
 221                     newEdgeI++;
 222                 }
 223             }
 224         }
 225
 226         edges_.setSize(newEdgeI);
 227
 228         for
 229         (
 230             StaticHashTable<label, edge, Hash<edge> >::const_iterator iter =
 231                 edgeToLabel.begin();
 232             iter != edgeToLabel.end();
 233             ++iter
 234         )
 235         {
 236             edges_[iter()] = iter.key();
 237         }
 238     }
 239 }
 240
 241
 242 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 243
 244 void Foam::edgeMesh::operator=(const edgeMesh& rhs)
 245 {
 246     points_ = rhs.points_;
 247     edges_ = rhs.edges_;
 248     pointEdgesPtr_.reset(NULL);
 249 }
 250
 251
 252 // ************************************************************************* //