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) 2011 OpenFOAM Foundation
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "mergePoints.H"
  27 #include "SortableList.H"
  28 #include "ListOps.H"
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 bool Foam::mergePoints
  33 (
  34     const UList<point>& points,
  35     const scalar mergeTol,
  36     const bool verbose,
  37     labelList& pointMap,
  38     List<point>& newPoints,
  39     const point& origin
  40 )
  41 {
  42     point compareOrigin = origin;
  43
  44     if (origin == point(VGREAT, VGREAT, VGREAT))
  45     {
  46         if (points.size())
  47         {
  48             compareOrigin = sum(points)/points.size();
  49         }
  50     }
  51
  52     // Create a old to new point mapping array
  53     pointMap.setSize(points.size());
  54     pointMap = -1;
  55
  56     // Storage for merged points
  57     newPoints.setSize(points.size());
  58
  59     if (points.empty())
  60     {
  61         return false;
  62     }
  63
  64     // We're comparing distance squared to origin first.
  65     // Say if starting from two close points:
  66     //     x, y, z
  67     //     x+mergeTol, y+mergeTol, z+mergeTol
  68     // Then the magSqr of both will be
  69     //     x^2+y^2+z^2
  70     //     x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
  71     // so the difference will be 2*mergeTol*(x+y+z)
  72
  73     const scalar mergeTolSqr = sqr(mergeTol);
  74
  75     // Sort points by magSqr
  76     const pointField d(points - compareOrigin);
  77     SortableList<scalar> sortedMagSqr(magSqr(d));
  78     scalarField sortedTol(points.size());
  79     forAll(sortedMagSqr.indices(), sortI)
  80     {
  81         const point& pt = d[sortedMagSqr.indices()[sortI]];
  82         sortedTol[sortI] = 2*mergeTol*(mag(pt.x())+mag(pt.y())+mag(pt.z()));
  83     }
  84
  85     bool hasMerged = false;
  86
  87     label newPointI = 0;
  88
  89
  90     // Handle 0th point separately (is always unique)
  91     label pointI = sortedMagSqr.indices()[0];
  92     pointMap[pointI] = newPointI;
  93     newPoints[newPointI++] = points[pointI];
  94
  95
  96     for (label sortI = 1; sortI < sortedMagSqr.size(); sortI++)
  97     {
  98         // Get original point index
  99         label pointI = sortedMagSqr.indices()[sortI];
 100
 101         // Compare to previous points to find equal one.
 102         label equalPointI = -1;
 103
 104         for
 105         (
 106             label prevSortI = sortI - 1;
 107             prevSortI >= 0
 108          && mag
 109             (
 110                 sortedMagSqr[prevSortI]
 111              -  sortedMagSqr[sortI]
 112             ) <= sortedTol[sortI];
 113             prevSortI--
 114         )
 115         {
 116             label prevPointI = sortedMagSqr.indices()[prevSortI];
 117
 118             if (magSqr(points[pointI] - points[prevPointI]) <= mergeTolSqr)
 119             {
 120                 // Found match.
 121                 equalPointI = prevPointI;
 122
 123                 break;
 124             }
 125         }
 126
 127
 128         if (equalPointI != -1)
 129         {
 130             // Same coordinate as equalPointI. Map to same new point.
 131             pointMap[pointI] = pointMap[equalPointI];
 132
 133             hasMerged = true;
 134
 135             if (verbose)
 136             {
 137                 Pout<< "Foam::mergePoints : Merging points "
 138                     << pointI << " and " << equalPointI
 139                     << " with coordinates:" << points[pointI]
 140                     << " and " << points[equalPointI]
 141                     << endl;
 142             }
 143         }
 144         else
 145         {
 146             // Differs. Store new point.
 147             pointMap[pointI] = newPointI;
 148             newPoints[newPointI++] = points[pointI];
 149         }
 150     }
 151
 152     newPoints.setSize(newPointI);
 153
 154     return hasMerged;
 155 }
 156
 157
 158 // ************************************************************************* //