applications/utilities/surface/surfaceAdd/surfaceAdd.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 Description
  26     Add two surfaces. Does geometric merge on points. Does not check for
  27     overlapping/intersecting triangles.
  28
  29     Keeps patches separate by renumbering.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "argList.H"
  34 #include "fileName.H"
  35 #include "triSurface.H"
  36 #include "OFstream.H"
  37 #include "IFstream.H"
  38 #include "triFace.H"
  39 #include "triFaceList.H"
  40
  41 using namespace Foam;
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 // Main program:
  46
  47 int main(int argc, char *argv[])
  48 {
  49     argList::noParallel();
  50     argList::validArgs.clear();
  51     argList::validArgs.append("Foam surface file");
  52     argList::validArgs.append("Foam surface file");
  53     argList::validArgs.append("Foam output file");
  54     argList::validOptions.insert("points", "pointsFile");
  55     argList::validOptions.insert("mergeRegions", "");
  56     argList args(argc, argv);
  57
  58     fileName inFileName1(args.additionalArgs()[0]);
  59     fileName inFileName2(args.additionalArgs()[1]);
  60     fileName outFileName(args.additionalArgs()[2]);
  61
  62     bool addPoint     = args.optionFound("points");
  63     bool mergeRegions = args.optionFound("mergeRegions");
  64
  65     if (addPoint)
  66     {
  67         Info<< "Reading a surface and adding points from a file"
  68             << "; merging the points and writing the surface to another file"
  69             << nl << endl;
  70
  71         Info<< "Surface  : " << inFileName1<< nl
  72             << "Points   : " << args.option("points") << nl
  73             << "Writing  : " << outFileName << nl << endl;
  74     }
  75     else
  76     {
  77         Info<< "Reading two surfaces"
  78             << "; merging points and writing the surface to another file"
  79             << nl << endl;
  80
  81         if (mergeRegions)
  82         {
  83             Info<< "Regions from the two files will get merged" << nl
  84                 << "Do not use this option if you want to keep the regions"
  85                 << " separate" << nl << endl;
  86         }
  87         else
  88         {
  89             Info<< "Regions from the two files will not get merged" << nl
  90                 << "Regions from " << inFileName2 << " will get offset so"
  91                 << " as not to overlap with the regions in " << inFileName1
  92                 << nl << endl;
  93         }
  94
  95
  96         Info<< "Surface1 : " << inFileName1<< nl
  97             << "Surface2 : " << inFileName2<< nl
  98             << "Writing  : " << outFileName << nl << endl;
  99     }
 100
 101     const triSurface surface1(inFileName1);
 102
 103     Info<< "Surface1:" << endl;
 104     surface1.writeStats(Info);
 105     Info<< endl;
 106
 107     const pointField& points1 = surface1.points();
 108
 109     // Final surface
 110     triSurface combinedSurf;
 111
 112     if (addPoint)
 113     {
 114         IFstream pointsFile(args.option("points"));
 115         pointField extraPoints(pointsFile);
 116
 117         Info<< "Additional Points:" << extraPoints.size() << endl;
 118
 119         vectorField pointsAll(points1);
 120         label pointI = pointsAll.size();
 121         pointsAll.setSize(pointsAll.size() + extraPoints.size());
 122
 123         forAll(extraPoints, i)
 124         {
 125             pointsAll[pointI++] = extraPoints[i];
 126         }
 127
 128         combinedSurf = triSurface(surface1, surface1.patches(), pointsAll);
 129     }
 130     else
 131     {
 132         const triSurface surface2(inFileName2);
 133
 134         Info<< "Surface2:" << endl;
 135         surface2.writeStats(Info);
 136         Info<< endl;
 137
 138
 139         // Make new storage
 140         List<labelledTri> facesAll(surface1.size() + surface2.size());
 141
 142         const pointField& points2 = surface2.points();
 143
 144         vectorField pointsAll(points1.size() + points2.size());
 145
 146
 147         label pointi = 0;
 148         // Copy points1 into pointsAll
 149         forAll(points1, point1i)
 150         {
 151             pointsAll[pointi++] = points1[point1i];
 152         }
 153         // Add surface2 points
 154         forAll(points2, point2i)
 155         {
 156             pointsAll[pointi++] = points2[point2i];
 157         }
 158
 159
 160         label trianglei = 0;
 161         label maxRegion1 = labelMin;
 162
 163         // Copy triangles1 into trianglesAll
 164         // Determine max region.
 165
 166         forAll(surface1, faceI)
 167         {
 168             facesAll[trianglei] = surface1[faceI];
 169
 170             maxRegion1 = max(maxRegion1, facesAll[trianglei].region());
 171
 172             trianglei++;
 173         }
 174
 175         label nRegions1 = maxRegion1 + 1;
 176
 177         if (!mergeRegions)
 178         {
 179             Info<< "Surface " << inFileName1 << " has " << nRegions1 << " regions"
 180                 << nl
 181                 << "All region numbers in " << inFileName2 << " will be offset"
 182                 << " by this amount" << nl << endl;
 183         }
 184
 185         // Add (renumbered) surface2 triangles
 186         label maxRegion2 = labelMin;
 187
 188         forAll(surface2, faceI)
 189         {
 190             const labelledTri& tri = surface2[faceI];
 191
 192             labelledTri& destTri = facesAll[trianglei++];
 193
 194             destTri[0] = tri[0] + points1.size();
 195             destTri[1] = tri[1] + points1.size();
 196             destTri[2] = tri[2] + points1.size();
 197
 198             maxRegion2  = max(maxRegion2, tri.region());
 199
 200             if (mergeRegions)
 201             {
 202                 destTri.region() = tri.region();
 203             }
 204             else
 205             {
 206                 destTri.region() = tri.region() + nRegions1;
 207             }
 208         }
 209
 210         label nRegions2 = maxRegion2 + 1;
 211
 212         geometricSurfacePatchList newPatches;
 213
 214         if (mergeRegions)
 215         {
 216             // Overwrite
 217             newPatches.setSize(max(nRegions1, nRegions2));
 218
 219             forAll(surface1.patches(), patchI)
 220             {
 221                 newPatches[patchI] = surface1.patches()[ patchI];
 222             }
 223             forAll(surface2.patches(), patchI)
 224             {
 225                 newPatches[patchI] = surface2.patches()[ patchI];
 226             }
 227         }
 228         else
 229         {
 230             Info<< "Regions from " << inFileName2 << " have been renumbered:"
 231                 << nl
 232                 << "    old\tnew" << nl;
 233
 234             for (label regionI = 0; regionI < nRegions2; regionI++)
 235             {
 236                 Info<< "    " << regionI << '\t' << regionI+nRegions1
 237                     << nl;
 238             }
 239             Info<< nl;
 240
 241             newPatches.setSize(nRegions1 + nRegions2);
 242
 243             label newPatchI = 0;
 244
 245             forAll(surface1.patches(), patchI)
 246             {
 247                 newPatches[newPatchI++] = surface1.patches()[ patchI];
 248             }
 249
 250             forAll(surface2.patches(), patchI)
 251             {
 252                 newPatches[newPatchI++] = surface2.patches()[ patchI];
 253             }
 254         }
 255
 256
 257         Info<< "New patches:" << nl;
 258         forAll(newPatches, patchI)
 259         {
 260             Info<< "    " << patchI << '\t' << newPatches[patchI].name() << nl;
 261         }
 262         Info<< endl;
 263
 264
 265         // Construct new surface mesh
 266         combinedSurf = triSurface(facesAll, newPatches, pointsAll);
 267     }
 268
 269     // Merge all common points and do some checks
 270     combinedSurf.cleanup(true);
 271
 272     Info<< "Merged surface:" << endl;
 273
 274     combinedSurf.writeStats(Info);
 275
 276     Info<< endl;
 277
 278     Info << "Writing : " << outFileName << endl;
 279
 280     // No need to 'group' while writing since all in correct order anyway.
 281     combinedSurf.write(outFileName);
 282
 283     Info << "End\n" << endl;
 284
 285     return 0;
 286 }
 287
 288
 289 // ************************************************************************* //