applications/utilities/surface/surfaceCoarsen/surfaceCoarsen.C

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