src/OpenFOAM/meshes/meshTools/mergePoints.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "mergePoints.H"
  28 #include "SortableList.H"
  29 #include "ListOps.H"
  30
  31 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  32
  33 bool Foam::mergePoints
  34 (
  35     const UList<point>& points,
  36     const scalar mergeTol,
  37     const bool verbose,
  38     labelList& pointMap,
  39     List<point>& newPoints,
  40     const point& origin
  41 )
  42 {
  43     point compareOrigin = origin;
  44
  45     if (origin == point(VGREAT, VGREAT, VGREAT))
  46     {
  47         if (points.size())
  48         {
  49             compareOrigin = sum(points)/points.size();
  50         }
  51     }
  52
  53     // Create a old to new point mapping array
  54     pointMap.setSize(points.size());
  55     pointMap = -1;
  56
  57     // Storage for merged points
  58     newPoints.setSize(points.size());
  59
  60     if (points.empty())
  61     {
  62         return false;
  63     }
  64
  65
  66     const scalar mergeTolSqr = sqr(mergeTol);
  67
  68     // Sort points by magSqr
  69     SortableList<scalar> sortedMagSqr(magSqr(points - compareOrigin));
  70
  71     bool hasMerged = false;
  72
  73     label newPointI = 0;
  74
  75
  76     // Handle 0th point separately (is always unique)
  77     label pointI = sortedMagSqr.indices()[0];
  78     pointMap[pointI] = newPointI;
  79     newPoints[newPointI++] = points[pointI];
  80
  81
  82     for (label sortI = 1; sortI < sortedMagSqr.size(); sortI++)
  83     {
  84         // Get original point index
  85         label pointI = sortedMagSqr.indices()[sortI];
  86
  87         // Compare to previous points to find equal one.
  88         label equalPointI = -1;
  89
  90         for
  91         (
  92             label prevSortI = sortI - 1;
  93             prevSortI >= 0
  94          && mag
  95             (
  96                 sortedMagSqr[prevSortI]
  97              -  sortedMagSqr[sortI]
  98             ) <= mergeTolSqr;
  99             prevSortI--
 100         )
 101         {
 102             label prevPointI = sortedMagSqr.indices()[prevSortI];
 103
 104             if (magSqr(points[pointI] - points[prevPointI]) <= mergeTolSqr)
 105             {
 106                 // Found match.
 107                 equalPointI = prevPointI;
 108
 109                 break;
 110             }
 111         }
 112
 113
 114         if (equalPointI != -1)
 115         {
 116             // Same coordinate as equalPointI. Map to same new point.
 117             pointMap[pointI] = pointMap[equalPointI];
 118
 119             hasMerged = true;
 120
 121             if (verbose)
 122             {
 123                 Pout<< "Foam::mergePoints : Merging points "
 124                     << pointI << " and " << equalPointI
 125                     << " with coordinates:" << points[pointI]
 126                     << " and " << points[equalPointI]
 127                     << endl;
 128             }
 129         }
 130         else
 131         {
 132             // Differs. Store new point.
 133             pointMap[pointI] = newPointI;
 134             newPoints[newPointI++] = points[pointI];
 135         }
 136     }
 137
 138     newPoints.setSize(newPointI);
 139
 140     return hasMerged;
 141 }
 142
 143
 144 // ************************************************************************* //