applications/utilities/mesh/conversion/starToFoam/mergeCoupleFacePoints.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  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.
  15
  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.
  20
  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
  24
  25 Description
  26     Merge duplicate points created by arbitrary face coupling and remove unused
  27     points.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "starMesh.H"
  32 #include "demandDrivenData.H"
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 void starMesh::mergeCoupleFacePoints()
  37 {
  38     // mark all used points by looping through all faces in two goes.
  39     // First, go into every cell and find min edge length.  Use a
  40     // fraction of that as a merge tolerance.  Loop through all the
  41     // points of the cell and query a merge against every other point
  42     // of the same cell based on the current tolerance.  If two points
  43     // merge, find out if any of them already belongs to a "merge set"
  44     // (group of points that merge together.  If so, add the current
  45     // points into the merge set and make sure that the merge label
  46     // for the whole set is the lowest of all encountered vertices.
  47     // This is stored in a mergeSetID list, under the index of the
  48     // merge set.  Once all cells (and thus points) are visited, go
  49     // through the renumbering list and for each merging point use the
  50     // label of the merge set as the new point label.
  51     // This is VERY fancy. Use care if/when changing.
  52
  53     Info << endl << "Creating merge sets" << endl;
  54
  55     // Create a renumbering list for points
  56     // In the first instance the renumbering list is used as a
  57     // mergeSetID storage
  58     labelList renumberPoints(points_.size(), -1);
  59
  60     // create storage of mergeSetIDs
  61     label mergeIncrement = points_.size()/10;
  62     labelList mergeSetID(mergeIncrement, -1);
  63
  64     label nMergeSets = 0;
  65
  66     forAll (cellFaces_, cellI)
  67     {
  68         const faceList& curFaces = cellFaces_[cellI];
  69
  70         // create a list of all points for the cell with duplicates
  71         // and find the shortest edge length
  72
  73         label nPointsInCell = 0;
  74
  75         scalar pointMergeTol = GREAT;
  76
  77         forAll (curFaces, faceI)
  78         {
  79             nPointsInCell += curFaces[faceI].size();
  80
  81             edgeList curEdges = curFaces[faceI].edges();
  82
  83             forAll (curEdges, edgeI)
  84             {
  85                 scalar length = curEdges[edgeI].mag(points_);
  86
  87                 if (length < pointMergeTol)
  88                 {
  89                     pointMergeTol = length;
  90                 }
  91             }
  92         }
  93
  94         // set merge tolerance as a fraction of the shortest edge
  95         pointMergeTol /= 100.0;
  96
  97         // create a list of points
  98         labelList cellPoints(nPointsInCell);
  99         label nAddedPoints = 0;
 100
 101         forAll (curFaces, faceI)
 102         {
 103             const face& f = curFaces[faceI];
 104
 105             forAll (f, fI)
 106             {
 107                 cellPoints[nAddedPoints] = f[fI];
 108                 nAddedPoints++;
 109             }
 110         }
 111
 112         // loop n-squared through the local points and merge them up
 113         for (label firstPI = 0; firstPI < cellPoints.size() - 1; firstPI++)
 114         {
 115             for
 116             (
 117                 label otherPI = firstPI;
 118                 otherPI < cellPoints.size();
 119                 otherPI++
 120             )
 121             {
 122                 if (cellPoints[firstPI] != cellPoints[otherPI])
 123                 {
 124                     label a = cellPoints[firstPI];
 125                     label b = cellPoints[otherPI];
 126
 127                     if (edge (a, b).mag(points_) < pointMergeTol)
 128                     {
 129                         // found a pair of points to merge
 130 #                       ifdef DEBUG_MERGE
 131                         Info<< "Merging points " << a << " and " << b << endl;
 132 #                       endif
 133
 134                         // are the two points in a merge group?
 135                         label mergeSetA = -1;
 136                         label mergeSetB = -1;
 137
 138                         if (renumberPoints[a] > -1)
 139                         {
 140                             mergeSetA = renumberPoints[a];
 141                         }
 142
 143                         if (renumberPoints[b] > -1)
 144                         {
 145                             mergeSetB = renumberPoints[b];
 146                         }
 147
 148                         if (mergeSetA == -1 && mergeSetB == -1)
 149                         {
 150                             // add new merge group
 151 #                           ifdef DEBUG_MERGE
 152                             Info<< "adding new merge group " << nMergeSets
 153                                 << endl;
 154 #                           endif
 155
 156                             // mark points as belonging to a new merge set
 157                             renumberPoints[a] = nMergeSets;
 158                             renumberPoints[b] = nMergeSets;
 159
 160                             mergeSetID[nMergeSets] = min(a, b);
 161                             nMergeSets++;
 162
 163                             if (nMergeSets >= mergeSetID.size())
 164                             {
 165                                 Info << "Resizing mergeSetID" << endl;
 166
 167                                 mergeSetID.setSize
 168                                     (mergeSetID.size() + mergeIncrement);
 169                             }
 170                         }
 171                         else if (mergeSetA == -1 && mergeSetB != -1)
 172                         {
 173 #                           ifdef DEBUG_MERGE
 174                             Info<< "adding point a into the merge set of b. "
 175                                 << "a: " << a << endl;
 176 #                           endif
 177
 178                             // add point a into the merge set of b
 179                             renumberPoints[a] = mergeSetB;
 180
 181                             // reset the min label of the mergeSet
 182                             mergeSetID[mergeSetB] =
 183                                 min(a, mergeSetID[mergeSetB]);
 184                         }
 185                         else if  (mergeSetA != -1 && mergeSetB == -1)
 186                         {
 187 #                           ifdef DEBUG_MERGE
 188                             Info<< "adding point b into the merge set of a. "
 189                                 << "b: " << b << endl;
 190 #                           endif
 191
 192                             // add point b into the merge set of a
 193                             renumberPoints[b] = mergeSetA;
 194
 195                             // reset the min label of the mergeSet
 196                             mergeSetID[mergeSetA] =
 197                                 min(b, mergeSetID[mergeSetA]);
 198                         }
 199                         else if (mergeSetA != mergeSetB)
 200                         {
 201                             // Points already belong to two different merge
 202                             // sets. Eliminate the higher merge set
 203                             label minMerge = min(mergeSetA, mergeSetB);
 204                             label maxMerge = max(mergeSetA, mergeSetB);
 205
 206 #                           ifdef DEBUG_MERGE
 207                             Info<< "Points already belong to two "
 208                                 << "different merge sets. "
 209                                 << "Eliminate the higher merge set. Sets: "
 210                                 << minMerge << " and " << maxMerge << endl;
 211 #                           endif
 212
 213                             forAll (renumberPoints, elimI)
 214                             {
 215                                 if (renumberPoints[elimI] == maxMerge)
 216                                 {
 217                                     renumberPoints[elimI] = minMerge;
 218                                 }
 219                             }
 220
 221                             // set the lower label
 222                             mergeSetID[minMerge] =
 223                                 min(mergeSetID[minMerge], mergeSetID[maxMerge]);
 224                         }
 225                     }
 226                 }
 227             }
 228         }
 229     }
 230
 231     mergeSetID.setSize(nMergeSets);
 232
 233     Info<< "Finished creating merge sets.  Number of merge sets: "
 234         << nMergeSets << "." << endl;
 235
 236     // Insert the primary point renumbering into the list
 237     // Take care of possibly unused points in the list
 238     forAll (renumberPoints, pointI)
 239     {
 240         if (renumberPoints[pointI] < 0)
 241         {
 242             // point not merged
 243             renumberPoints[pointI] = pointI;
 244         }
 245         else
 246         {
 247             renumberPoints[pointI] = mergeSetID[renumberPoints[pointI]];
 248         }
 249     }
 250
 251     // Now every point carries either its own label or the lowest of
 252     // the labels of the points it is merged with. Now do PRELIMINARY
 253     // renumbering of all faces.  This will only be used to see which
 254     // points are still used!
 255
 256     forAll (cellFaces_, cellI)
 257     {
 258         faceList& prelimFaces = cellFaces_[cellI];
 259
 260         forAll (prelimFaces, faceI)
 261         {
 262             face oldFacePoints = prelimFaces[faceI];
 263
 264             face& prelimFacePoints = prelimFaces[faceI];
 265
 266             forAll (prelimFacePoints, pointI)
 267             {
 268                 if (renumberPoints[oldFacePoints[pointI]] < 0)
 269                 {
 270                     FatalErrorIn("starMesh::mergeCoupleFacePoints()")
 271                         << "Error in point renumbering. Old face: "
 272                         << oldFacePoints << endl
 273                         << "prelim face: " << prelimFacePoints
 274                         << abort(FatalError);
 275                 }
 276
 277                 prelimFacePoints[pointI] =
 278                     renumberPoints[oldFacePoints[pointI]];
 279             }
 280         }
 281     }
 282
 283     // First step complete. Reset the renumbering list, mark all used points,
 284     // re-create the point list and renumber the whole lot
 285     renumberPoints = 0;
 286
 287     forAll (cellFaces_, cellI)
 288     {
 289         const faceList& curFaces = cellFaces_[cellI];
 290
 291         forAll (curFaces, faceI)
 292         {
 293             const face& curFacePoints = curFaces[faceI];
 294
 295             forAll (curFacePoints, pointI)
 296             {
 297                 renumberPoints[curFacePoints[pointI]]++;
 298             }
 299         }
 300     }
 301
 302     forAll (cellShapes_, cellI)
 303     {
 304         const labelList& curLabels = cellShapes_[cellI];
 305
 306         forAll (curLabels, pointI)
 307         {
 308             if (renumberPoints[curLabels[pointI]] == 0)
 309             {
 310                 FatalErrorIn("starMesh::mergeCoupleFacePoints()")
 311                     << "Error in point merging for cell "
 312                     << cellI << ". STAR index: " << starCellID_[cellI]
 313                     << ". " << endl
 314                     << "Point index: " << curLabels[pointI] << " STAR index "
 315                     << starPointID_[curLabels[pointI]] << endl
 316                     << "Please check the geometry for the cell." << endl;
 317             }
 318         }
 319     }
 320
 321     label nUsedPoints = 0;
 322
 323     forAll (renumberPoints, pointI)
 324     {
 325         if (renumberPoints[pointI] > 0)
 326         {
 327             // This should be OK as the compressed points list will always
 328             // have less points that the original lists.  Even if there is
 329             // no points removed, this will copy the list back onto itself
 330             //
 331             renumberPoints[pointI] = nUsedPoints;
 332             points_[nUsedPoints] = points_[pointI];
 333
 334             nUsedPoints++;
 335         }
 336         else
 337         {
 338             renumberPoints[pointI] = -1;
 339         }
 340     }
 341
 342     Info<< "Total number of points: " << points_.size() << endl
 343         << "Number of used points: " << nUsedPoints << endl;
 344
 345     // reset number of points which need to be sorted
 346     points_.setSize(nUsedPoints);
 347
 348     Info << "Renumbering all faces" << endl;
 349
 350     forAll (cellFaces_, cellI)
 351     {
 352         faceList& newFaces = cellFaces_[cellI];
 353
 354         forAll (newFaces, faceI)
 355         {
 356             face oldFacePoints = newFaces[faceI];
 357
 358             face& newFacePoints = newFaces[faceI];
 359
 360             forAll (newFacePoints, pointI)
 361             {
 362                 if (renumberPoints[oldFacePoints[pointI]] < 0)
 363                 {
 364                     FatalErrorIn("starMesh::mergeCoupleFacePoints()")
 365                         << "Error in point renumbering for point "
 366                         << oldFacePoints[pointI]
 367                         << ". Renumbering index is -1." << endl
 368                         << "Old face: " << oldFacePoints << endl
 369                         << "New face: " << newFacePoints << abort(FatalError);
 370                 }
 371
 372                 newFacePoints[pointI] = renumberPoints[oldFacePoints[pointI]];
 373             }
 374         }
 375     }
 376
 377     Info << "Renumbering all cell shapes" << endl;
 378
 379     forAll (cellShapes_, cellI)
 380     {
 381         labelList oldLabels = cellShapes_[cellI];
 382
 383         labelList& curLabels = cellShapes_[cellI];
 384
 385         forAll (curLabels, pointI)
 386         {
 387             if (renumberPoints[curLabels[pointI]] < 0)
 388             {
 389                 FatalErrorIn("starMesh::mergeCoupleFacePoints()")
 390                     << "Error in point renumbering for cell "
 391                     << cellI << ". STAR index: " << starCellID_[cellI]
 392                     << ". " << endl
 393                     << "Point index: " << curLabels[pointI] << " STAR index "
 394                     << starPointID_[curLabels[pointI]] << " returns invalid "
 395                     << "renumbering index: "
 396                     << renumberPoints[curLabels[pointI]] << "." << endl
 397                     << "Old cellShape:  " << oldLabels << endl
 398                     << "New cell shape: " << curLabels << abort(FatalError);
 399                 }
 400
 401             curLabels[pointI] = renumberPoints[oldLabels[pointI]];
 402         }
 403     }
 404
 405     Info << "Renumbering STAR point lookup" << endl;
 406
 407     labelList oldStarPointID = starPointID_;
 408
 409     starPointID_ = -1;
 410
 411     forAll (starPointID_, pointI)
 412     {
 413         if (renumberPoints[pointI] > -1)
 414         {
 415             starPointID_[renumberPoints[pointI]] = oldStarPointID[pointI];
 416         }
 417     }
 418
 419     // point-cell addressing has changed. Force it to be re-created
 420     deleteDemandDrivenData(pointCellsPtr_);
 421 }
 422
 423
 424 // ************************************************************************* //