applications/utilities/surface/surfaceCoarsen/surfaceCoarsen.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 Description
  25     Surface coarsening using 'bunnylod':
  26
  27         Polygon Reduction Demo
  28         By Stan Melax (c) 1998
  29         mailto:melax@cs.ualberta.ca
  30         http://www.cs.ualberta.ca/~melax
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "fileName.H"
  36 #include "triSurface.H"
  37 #include "OFstream.H"
  38 #include "triFace.H"
  39 #include "triFaceList.H"
  40
  41 // From bunnylod
  42 #include "progmesh.h"
  43
  44 using namespace Foam;
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 int mapVertex(::List<int>& collapse_map, int a, int mx)
  49 {
  50     if (mx <= 0)
  51     {
  52         return 0;
  53     }
  54     while (a >= mx)
  55     {
  56         a = collapse_map[a];
  57     }
  58     return a;
  59 }
  60
  61
  62 // Main program:
  63
  64 int main(int argc, char *argv[])
  65 {
  66     argList::noParallel();
  67     argList::validArgs.append("surfaceFile");
  68     argList::validArgs.append("reductionFactor");
  69     argList::validArgs.append("output surfaceFile");
  70     argList args(argc, argv);
  71
  72     const fileName inFileName = args[1];
  73     const scalar reduction = args.argRead<scalar>(2);
  74     const fileName outFileName = args[3];
  75
  76     if (reduction <= 0 || reduction > 1)
  77     {
  78         FatalErrorIn(args.executable())
  79             << "Reduction factor " << reduction
  80             << " should be within 0..1" << endl
  81             << "(it is the reduction in number of vertices)"
  82             << exit(FatalError);
  83     }
  84
  85     Info<< "Input surface   :" << inFileName << endl
  86         << "Reduction factor:" << reduction << endl
  87         << "Output surface  :" << outFileName << endl << endl;
  88
  89     const triSurface surf(inFileName);
  90
  91     Info<< "Surface:" << endl;
  92     surf.writeStats(Info);
  93     Info<< endl;
  94
  95
  96     ::List< ::Vector> vert;     // global list of vertices
  97     ::List< ::tridata> tri;     // global list of triangles
  98
  99
 100     // Convert triSurface to progmesh format. Note: can use global point
 101     // numbering since surface read in from file.
 102     const pointField& pts = surf.points();
 103
 104     forAll(pts, ptI)
 105     {
 106         const point& pt = pts[ptI];
 107
 108         vert.Add( ::Vector(pt.x(), pt.y(), pt.z()));
 109     }
 110
 111     forAll(surf, faceI)
 112     {
 113         const labelledTri& f = surf[faceI];
 114
 115         tridata td;
 116         td.v[0]=f[0];
 117         td.v[1]=f[1];
 118         td.v[2]=f[2];
 119         tri.Add(td);
 120     }
 121
 122     ::List<int> collapse_map;   // to which neighbor each vertex collapses
 123     ::List<int> permutation;
 124
 125     ::ProgressiveMesh(vert,tri,collapse_map,permutation);
 126
 127     // rearrange the vertex list
 128     ::List< ::Vector> temp_list;
 129     for (int i=0;i<vert.num;i++)
 130     {
 131         temp_list.Add(vert[i]);
 132     }
 133     for (int i=0;i<vert.num;i++)
 134     {
 135         vert[permutation[i]]=temp_list[i];
 136     }
 137
 138     // update the changes in the entries in the triangle list
 139     for (int i=0;i<tri.num;i++)
 140     {
 141         for (int j=0;j<3;j++)
 142         {
 143             tri[i].v[j] = permutation[tri[i].v[j]];
 144         }
 145     }
 146
 147     // Only get triangles with non-collapsed edges.
 148     int render_num = int(reduction * surf.nPoints());
 149
 150     Info<< "Reducing to " << render_num << " vertices" << endl;
 151
 152
 153     // Storage for new surface.
 154     Foam::List<labelledTri> newTris(surf.size());
 155
 156     label newI = 0;
 157
 158     for (int i=0; i<tri.num; i++)
 159     {
 160         int p0 = mapVertex(collapse_map, tri[i].v[0], render_num);
 161         int p1 = mapVertex(collapse_map, tri[i].v[1], render_num);
 162         int p2 = mapVertex(collapse_map, tri[i].v[2], render_num);
 163
 164         // note:  serious optimization opportunity here,
 165         //  by sorting the triangles the following "continue"
 166         //  could have been made into a "break" statement.
 167         if (p0 == p1 || p1 == p2 || p2 == p0)
 168         {
 169             continue;
 170         }
 171
 172         newTris[newI++] = labelledTri(p0, p1, p2, 0);
 173     }
 174     newTris.setSize(newI);
 175
 176     // Convert vert into pointField.
 177     pointField newPoints(vert.num);
 178
 179     for (int i=0; i<vert.num; i++)
 180     {
 181         const ::Vector & v = vert[i];
 182
 183         newPoints[i] = point(v.x, v.y, v.z);
 184     }
 185
 186     triSurface surf2(newTris, newPoints);
 187
 188     triSurface outSurf
 189     (
 190         surf2.localFaces(),
 191         surf2.patches(),
 192         surf2.localPoints()
 193     );
 194
 195     Info<< "Coarsened surface:" << endl;
 196     surf2.writeStats(Info);
 197     Info<< endl;
 198
 199     Info<< "Writing to file " << outFileName << endl << endl;
 200
 201     surf2.write(outFileName);
 202
 203     Info<< "End\n" << endl;
 204
 205     return 0;
 206 }
 207
 208
 209 // ************************************************************************* //