applications/test/findCell-octree/findSubSet.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
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "fvCFD.H"
  30 #include "IStringStream.H"
  31
  32 #include "myBoundBox.H"
  33 #include "myBoundBoxList.H"
  34 #include "octree.H"
  35 #include "octreeData.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38 // Main program:
  39
  40
  41
  42 int main(int argc, char *argv[])
  43 {
  44     argList::validOptions.insert("x1", "X1");
  45     argList::validOptions.insert("y1", "Y1");
  46     argList::validOptions.insert("z1", "Z1");
  47
  48     argList::validOptions.insert("x2", "X2");
  49     argList::validOptions.insert("y2", "Y2");
  50     argList::validOptions.insert("z2", "Z2");
  51
  52 #   include "setRootCase.H"
  53 #   include "createTime.H"
  54 #   include "createMesh.H"
  55
  56     // Calculate BB of all cells
  57
  58
  59     myBoundBoxList allBb(mesh.nCells());
  60
  61     const pointField& allPoints = mesh.points();
  62
  63      vectorField bbMin(mesh.nCells());
  64      bbMin = vector(GREAT, GREAT, GREAT);
  65      vectorField bbMax(mesh.nCells());
  66      bbMax = vector(-GREAT, -GREAT, -GREAT);
  67
  68      const labelListList& pCells = mesh.pointCells();
  69
  70      forAll(pCells, pointi)
  71      {
  72          const point& vertCoord = allPoints[pointi];
  73
  74          const labelList& cells = pCells[pointi];
  75
  76          forAll(cells, celli)
  77          {
  78              label cellNum = cells[celli];
  79
  80              bbMin[cellNum].x() = min(bbMin[cellNum].x(), vertCoord.x());
  81              bbMin[cellNum].y() = min(bbMin[cellNum].y(), vertCoord.y());
  82              bbMin[cellNum].z() = min(bbMin[cellNum].z(), vertCoord.z());
  83
  84              bbMax[cellNum].x() = max(bbMax[cellNum].x(), vertCoord.x());
  85              bbMax[cellNum].y() = max(bbMax[cellNum].y(), vertCoord.y());
  86              bbMax[cellNum].z() = max(bbMax[cellNum].z(), vertCoord.z());
  87          }
  88      }
  89
  90
  91     forAll(allBb, celli)
  92     {
  93         allBb[celli] = myBoundBox(bbMin[celli], bbMax[celli]);
  94     }
  95
  96
  97     myBoundBox meshBb(allPoints);
  98
  99     scalar typDim = meshBb.minDim()/111;
 100
 101     myBoundBox shiftedBb
 102     (
 103         meshBb.min(),
 104         point
 105         (
 106             meshBb.max().x() + typDim,
 107             meshBb.max().y() + typDim,
 108             meshBb.max().z() + typDim
 109         )
 110     );
 111
 112
 113     Info<< "Mesh" << endl;
 114     Info<< "   bounding box     :" << shiftedBb << endl;
 115     Info<< "   typical dimension:" << shiftedBb.typDim() << endl;
 116
 117
 118     /*
 119      * Now we have allBb and shiftedBb
 120      */
 121
 122
 123
 124     // Construct table of subset of cells
 125
 126     labelList cellIndices(10);
 127
 128     cellIndices[0] = 1433;
 129     cellIndices[1] = 1434;
 130     cellIndices[2] = 1435;
 131     cellIndices[3] = 1436;
 132     cellIndices[4] = 1437;
 133     cellIndices[5] = 1438;
 134     cellIndices[6] = 1439;
 135     cellIndices[7] = 1440;
 136     cellIndices[8] = 1441;
 137     cellIndices[9] = 1442;
 138
 139     // Get the corresponding bounding boxes
 140
 141     forAll(cellIndices, i)
 142     {
 143         allBb[i] = allBb[cellIndices[i]];
 144     }
 145     allBb.setSize(cellIndices.size());
 146
 147
 148
 149     // Wrap indices and mesh information into helper object
 150     octreeData shapes(mesh, cellIndices);
 151
 152     octree oc
 153     (
 154         shiftedBb,  // overall bounding box
 155         shapes,     // all information needed to do checks on cells
 156         allBb,      // bounding boxes of cells
 157         10.0        // maximum ratio of cubes v.s. cells
 158     );
 159
 160 //    scalar x1(readScalar(IStringStream(args.options()["x1"])()));
 161 //    scalar y1(readScalar(IStringStream(args.options()["y1"])()));
 162 //    scalar z1(readScalar(IStringStream(args.options()["z1"])()));
 163
 164 //    scalar x2(readScalar(IStringStream(args.options()["x2"])()));
 165 //    scalar y2(readScalar(IStringStream(args.options()["y2"])()));
 166 //    scalar z2(readScalar(IStringStream(args.options()["z2"])()));
 167
 168
 169
 170     label nFound = 0;
 171
 172     scalar x = -5.0;
 173     for(int i = 0; i < 100; i++)
 174     {
 175         scalar y = -7.0;
 176         for(int j = 0; j < 10; j++)
 177         {
 178             scalar z = -12.0;
 179             for (int k = 0; k < 10; k++)
 180             {
 181                 point sample(x, y, z);
 182
 183                 label index = oc.find(sample);
 184
 185                 // Convert index into shapes back into cellindex.
 186                 label cell;
 187                 if (index != -1)
 188                 {
 189                     cell = cellIndices[index];
 190                 }
 191                 else
 192                 {
 193                     cell = -1;
 194                 }
 195                 Info<< "Point:" << sample
 196                     << " is in cell " << cell << "(octree)  "
 197                     << mesh.findCell(sample) << "(linear)"
 198                     << endl;
 199
 200                 z += 1.2;
 201             }
 202             y += 0.9;
 203         }
 204         x += 0.1;
 205     }
 206
 207
 208     Info<< "nFound=" << nFound << endl;
 209
 210     Info<< "End\n" << endl;
 211
 212
 213     Info<< "Statistics:" << endl
 214         << "  nCells   :" << allBb.size() << endl
 215         << "  nNodes   :" << oc.nNodes() << endl
 216         << "  nLeaves  :" << oc.nLeaves() << endl
 217         << "  nEntries :" << oc.nEntries() << endl
 218         << "  Cells per leaf :"
 219         << oc.nEntries()/oc.nLeaves()
 220         << endl
 221         << "  Every cell in  :"
 222         << oc.nEntries()/allBb.size() << " cubes"
 223         << endl;
 224
 225     return 0;
 226 }
 227
 228
 229 // ************************************************************************* //