applications/utilities/surface/surfaceCoarsen/surfaceCoarsen.C

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