applications/test/extendedStencil/Test-ExtendedStencil.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 Application
  25     testExtendedStencil
  26
  27 Description
  28     Test app for determining extended stencil.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "argList.H"
  33 #include "fvMesh.H"
  34 #include "volFields.H"
  35 #include "Time.H"
  36 //#include "mapDistribute.H"
  37 #include "OFstream.H"
  38 #include "meshTools.H"
  39 //#include "FECCellToFaceStencil.H"
  40 //#include "CFCCellToFaceStencil.H"
  41 //#include "CPCCellToFaceStencil.H"
  42 //#include "CECCellToFaceStencil.H"
  43 //#include "extendedCentredCellToFaceStencil.H"
  44 //#include "extendedUpwindCellToFaceStencil.H"
  45
  46 //#include "centredCFCCellToFaceStencilObject.H"
  47 //#include "centredFECCellToFaceStencilObject.H"
  48 //#include "centredCPCCellToFaceStencilObject.H"
  49 //#include "centredCECCellToFaceStencilObject.H"
  50
  51 //#include "upwindFECCellToFaceStencilObject.H"
  52 //#include "upwindCPCCellToFaceStencilObject.H"
  53 //#include "upwindCECCellToFaceStencilObject.H"
  54
  55 //#include "upwindCFCCellToFaceStencilObject.H"
  56 #include "centredCFCFaceToCellStencilObject.H"
  57
  58 using namespace Foam;
  59
  60 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  61
  62 void writeStencilOBJ
  63 (
  64     const fileName& fName,
  65     const point& fc,
  66     const List<point>& stencilCc
  67 )
  68 {
  69     OFstream str(fName);
  70     label vertI = 0;
  71
  72     meshTools::writeOBJ(str, fc);
  73     vertI++;
  74
  75     forAll(stencilCc, i)
  76     {
  77         meshTools::writeOBJ(str, stencilCc[i]);
  78         vertI++;
  79         str << "l 1 " << vertI << nl;
  80     }
  81 }
  82
  83
  84 // Stats
  85 void writeStencilStats(const labelListList& stencil)
  86 {
  87     label sumSize = 0;
  88     label nSum = 0;
  89     label minSize = labelMax;
  90     label maxSize = labelMin;
  91
  92     forAll(stencil, i)
  93     {
  94         const labelList& sCells = stencil[i];
  95
  96         if (sCells.size() > 0)
  97         {
  98             sumSize += sCells.size();
  99             nSum++;
 100             minSize = min(minSize, sCells.size());
 101             maxSize = max(maxSize, sCells.size());
 102         }
 103     }
 104     reduce(sumSize, sumOp<label>());
 105     reduce(nSum, sumOp<label>());
 106     sumSize /= nSum;
 107
 108     reduce(minSize, minOp<label>());
 109     reduce(maxSize, maxOp<label>());
 110
 111     Info<< "Stencil size :" << nl
 112         << "    average : " << sumSize << nl
 113         << "    min     : " << minSize << nl
 114         << "    max     : " << maxSize << nl
 115         << endl;
 116 }
 117
 118
 119 // Main program:
 120
 121 int main(int argc, char *argv[])
 122 {
 123 #   include "addTimeOptions.H"
 124 #   include "setRootCase.H"
 125 #   include "createTime.H"
 126
 127     // Get times list
 128     instantList Times = runTime.times();
 129 #   include "checkTimeOptions.H"
 130     runTime.setTime(Times[startTime], startTime);
 131 #   include "createMesh.H"
 132
 133     // Force calculation of extended edge addressing
 134     const labelListList& edgeFaces = mesh.edgeFaces();
 135     const labelListList& edgeCells = mesh.edgeCells();
 136     const labelListList& pointCells = mesh.pointCells();
 137     Info<< "dummy:" << edgeFaces.size() + edgeCells.size() + pointCells.size()
 138         << endl;
 139
 140
 141     // Centred, semi-extended stencil (edge cells only)
 142     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 143
 144 //    {
 145 //        //const FECCellToFaceStencil cfcStencil(mesh);
 146 //        //const extendedCentredCellToFaceStencil addressing
 147 //        //(
 148 //        //    cfcStencil
 149 //        //);
 150 //        const extendedCentredStencil& addressing =
 151 //        centredFECCellToFaceStencilObject::New
 152 //        (
 153 //            mesh
 154 //        );
 155 //
 156 //        Info<< "faceEdgeCell:" << endl;
 157 //        writeStencilStats(addressing.stencil());
 158 //
 159 //        // Collect stencil cell centres
 160 //        List<List<point> > stencilPoints(mesh.nFaces());
 161 //        addressing.collectData
 162 //        (
 163 //            mesh.C(),
 164 //            stencilPoints
 165 //        );
 166 //
 167 //        forAll(stencilPoints, faceI)
 168 //        {
 169 //            writeStencilOBJ
 170 //            (
 171 //                runTime.path()/"faceEdgeCell" + Foam::name(faceI) + ".obj",
 172 //                mesh.faceCentres()[faceI],
 173 //                stencilPoints[faceI]
 174 //            );
 175 //        }
 176 //    }
 177
 178
 179
 180
 181 //    // Centred, face stencil
 182 //    // ~~~~~~~~~~~~~~~~~~~~~
 183 //
 184 //    {
 185 //        const extendedCentredCellToFaceStencil& addressing =
 186 //        centredCFCCellToFaceStencilObject::New
 187 //        (
 188 //            mesh
 189 //        );
 190 //
 191 //        Info<< "cellFaceCell:" << endl;
 192 //        writeStencilStats(addressing.stencil());
 193 //
 194 //
 195 //        //// Do some interpolation.
 196 //        //{
 197 //        //    const labelListList& stencil = addressing.stencil();
 198 //        //    List<List<scalar> > stencilWeights(stencil.size());
 199 //        //    forAll(stencil, faceI)
 200 //        //    {
 201 //        //        const labelList& fStencil = stencil[faceI];
 202 //        //
 203 //        //        if (fStencil.size() > 0)
 204 //        //        {
 205 //        //            // Uniform weights
 206 //        //            stencilWeights[faceI] = scalarList
 207 //        //            (
 208 //        //                fStencil.size(),
 209 //        //                1.0/fStencil.size()
 210 //        //            );
 211 //        //        }
 212 //        //    }
 213 //        //
 214 //        //    tmp<surfaceVectorField> tfc
 215 //        //    (
 216 //        //        addressing.weightedSum(mesh.C(), stencilWeights)
 217 //        //    );
 218 //        //}
 219 //
 220 //
 221 //        // Collect stencil cell centres
 222 //        List<List<point> > stencilPoints(mesh.nFaces());
 223 //        addressing.collectData
 224 //        (
 225 //            mesh.C(),
 226 //            stencilPoints
 227 //        );
 228 //
 229 //        forAll(stencilPoints, faceI)
 230 //        {
 231 //            if (stencilPoints[faceI].size() >= 15)
 232 //            {
 233 //                writeStencilOBJ
 234 //                (
 235 //                    runTime.path()/"centredFace" + Foam::name(faceI) + ".obj",
 236 //                    mesh.faceCentres()[faceI],
 237 //                    stencilPoints[faceI]
 238 //                );
 239 //            }
 240 //        }
 241 //    }
 242
 243
 244 //    // Centred, point stencil
 245 //    // ~~~~~~~~~~~~~~~~~~~~~~
 246 //
 247 //    {
 248 //        //const extendedCentredCellToFaceStencil& addressing =
 249 //        //centredCPCStencilObject::New
 250 //        //(
 251 //        //    mesh
 252 //        //);
 253 //        //
 254 //        //Info<< "cellPointCell:" << endl;
 255 //        //writeStencilStats(addressing.stencil());
 256 //        //
 257 //        //
 258 //        //// Collect stencil cell centres
 259 //        //List<List<point> > stencilPoints(mesh.nFaces());
 260 //        //addressing.collectData
 261 //        //(
 262 //        //    mesh.C(),
 263 //        //    stencilPoints
 264 //        //);
 265 //        //
 266 //        //forAll(stencilPoints, faceI)
 267 //        //{
 268 //        //    writeStencilOBJ
 269 //        //    (
 270 //        //        runTime.path()/"centredPoint" + Foam::name(faceI) + ".obj",
 271 //        //        mesh.faceCentres()[faceI],
 272 //        //        stencilPoints[faceI]
 273 //        //    );
 274 //        //}
 275 //    }
 276
 277
 278
 279 //    // Centred, edge stencil
 280 //    // ~~~~~~~~~~~~~~~~~~~~~~
 281 //
 282 //    {
 283 //        //const extendedCentredCellToFaceStencil& addressing =
 284 //        //centredCECStencilObject::New
 285 //        //(
 286 //        //    mesh
 287 //        //);
 288 //        //
 289 //        //Info<< "cellEdgeCell:" << endl;
 290 //        //writeStencilStats(addressing.stencil());
 291 //        //
 292 //        //
 293 //        //// Collect stencil cell centres
 294 //        //List<List<point> > stencilPoints(mesh.nFaces());
 295 //        //addressing.collectData
 296 //        //(
 297 //        //    mesh.C(),
 298 //        //    stencilPoints
 299 //        //);
 300 //        //
 301 //        //forAll(stencilPoints, faceI)
 302 //        //{
 303 //        //    writeStencilOBJ
 304 //        //    (
 305 //        //        runTime.path()/"centredEdge" + Foam::name(faceI) + ".obj",
 306 //        //        mesh.faceCentres()[faceI],
 307 //        //        stencilPoints[faceI]
 308 //        //    );
 309 //        //}
 310 //    }
 311
 312
 313
 314     // Upwind, semi-extended stencil
 315     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 316
 317     //{
 318     //    const extendedUpwindCellToFaceStencil& addressing =
 319     //    upwindFECCellToFaceStencilObject::New
 320     //    (
 321     //        mesh,
 322     //        0.5
 323     //    );
 324     //
 325     //    Info<< "upwind-faceEdgeCell:" << endl;
 326     //    writeStencilStats(addressing.ownStencil());
 327     //
 328     //    {
 329     //        // Collect stencil cell centres
 330     //        List<List<point> > ownPoints(mesh.nFaces());
 331     //        addressing.collectData
 332     //        (
 333     //            addressing.ownMap(),
 334     //            addressing.ownStencil(),
 335     //            mesh.C(),
 336     //            ownPoints
 337     //        );
 338     //
 339     //        forAll(ownPoints, faceI)
 340     //        {
 341     //            writeStencilOBJ
 342     //            (
 343     //                runTime.path()/"ownFEC" + Foam::name(faceI) + ".obj",
 344     //                mesh.faceCentres()[faceI],
 345     //                ownPoints[faceI]
 346     //            );
 347     //        }
 348     //    }
 349     //    {
 350     //        // Collect stencil cell centres
 351     //        List<List<point> > neiPoints(mesh.nFaces());
 352     //        addressing.collectData
 353     //        (
 354     //            addressing.neiMap(),
 355     //            addressing.neiStencil(),
 356     //            mesh.C(),
 357     //            neiPoints
 358     //        );
 359     //
 360     //        forAll(neiPoints, faceI)
 361     //        {
 362     //            writeStencilOBJ
 363     //            (
 364     //                runTime.path()/"neiFEC" + Foam::name(faceI) + ".obj",
 365     //                mesh.faceCentres()[faceI],
 366     //                neiPoints[faceI]
 367     //            );
 368     //        }
 369     //    }
 370     //}
 371
 372
 373
 374     // Upwind, extended stencil
 375     // ~~~~~~~~~~~~~~~~~~~~~~~~
 376
 377     //{
 378     //    const extendedUpwindCellToFaceStencil& addressing =
 379     //    upwindCFCCellToFaceStencilObject::New
 380     //    (
 381     //        mesh,
 382     //        0.5
 383     //    );
 384     //
 385     //    Info<< "upwind-cellFaceCell:" << endl;
 386     //    writeStencilStats(addressing.ownStencil());
 387     //
 388     //    {
 389     //        // Collect stencil cell centres
 390     //        List<List<point> > ownPoints(mesh.nFaces());
 391     //        addressing.collectData
 392     //        (
 393     //            addressing.ownMap(),
 394     //            addressing.ownStencil(),
 395     //            mesh.C(),
 396     //            ownPoints
 397     //        );
 398     //
 399     //        forAll(ownPoints, faceI)
 400     //        {
 401     //            writeStencilOBJ
 402     //            (
 403     //                runTime.path()/"ownCFC" + Foam::name(faceI) + ".obj",
 404     //                mesh.faceCentres()[faceI],
 405     //                ownPoints[faceI]
 406     //            );
 407     //        }
 408     //    }
 409     //    {
 410     //        // Collect stencil cell centres
 411     //        List<List<point> > neiPoints(mesh.nFaces());
 412     //        addressing.collectData
 413     //        (
 414     //            addressing.neiMap(),
 415     //            addressing.neiStencil(),
 416     //            mesh.C(),
 417     //            neiPoints
 418     //        );
 419     //
 420     //        forAll(neiPoints, faceI)
 421     //        {
 422     //            writeStencilOBJ
 423     //            (
 424     //                runTime.path()/"neiCFC" + Foam::name(faceI) + ".obj",
 425     //                mesh.faceCentres()[faceI],
 426     //                neiPoints[faceI]
 427     //            );
 428     //        }
 429     //    }
 430     //}
 431
 432
 433
 434     //---- CELL CENTRED STENCIL -----
 435
 436     // Centred, cell stencil
 437     // ~~~~~~~~~~~~~~~~~~~~~
 438
 439     {
 440         const extendedCentredFaceToCellStencil& addressing =
 441         centredCFCFaceToCellStencilObject::New
 442         (
 443             mesh
 444         );
 445
 446         Info<< "cellFaceCell:" << endl;
 447         writeStencilStats(addressing.stencil());
 448
 449         // Collect stencil face centres
 450         List<List<point> > stencilPoints(mesh.nCells());
 451         addressing.collectData
 452         (
 453             mesh.Cf(),
 454             stencilPoints
 455         );
 456
 457         forAll(stencilPoints, cellI)
 458         {
 459             writeStencilOBJ
 460             (
 461                 runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
 462                 mesh.cellCentres()[cellI],
 463                 stencilPoints[cellI]
 464             );
 465         }
 466     }
 467
 468
 469 //XXXXXX
 470 //    // Evaluate
 471 //    List<List<scalar> > stencilData(faceStencils.size());
 472 //    collectStencilData
 473 //    (
 474 //        distMap,
 475 //        faceStencils,
 476 //        vf,
 477 //        stencilData
 478 //    );
 479 //    for (label faci = 0; faci < mesh.nInternalFaces(); faci++)
 480 //    {
 481 //        const scalarList& stData = stencilData[faceI];
 482 //        const scalarList& stWeight = fit[faceI];
 483 //
 484 //        forAll(stData, i)
 485 //        {
 486 //            sf[faceI] += stWeight[i]*stData[i];
 487 //        }
 488 //    }
 489 //    See finiteVolume/lnInclude/leastSquaresGrad.C
 490 //XXXXXX
 491
 492     Pout<< "End\n" << endl;
 493
 494     return 0;
 495 }
 496
 497
 498 // ************************************************************************* //