src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011-2011 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
  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 \*---------------------------------------------------------------------------*/
  25
  26 #include "refinementFeatures.H"
  27 #include "Time.H"
  28
  29 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  30
  31 Foam::refinementFeatures::refinementFeatures
  32 (
  33     const objectRegistry& io,
  34     const PtrList<dictionary>& featDicts
  35 )
  36 :
  37     PtrList<featureEdgeMesh>(featDicts.size()),
  38     levels_(featDicts.size()),
  39     edgeTrees_(featDicts.size()),
  40     pointTrees_(featDicts.size())
  41 {
  42     // Read features
  43
  44     forAll(featDicts, i)
  45     {
  46         const dictionary& dict = featDicts[i];
  47
  48         fileName featFileName(dict.lookup("file"));
  49
  50         set
  51         (
  52             i,
  53             new featureEdgeMesh
  54             (
  55                 IOobject
  56                 (
  57                     featFileName,                       // name
  58                     io.time().constant(),               // instance
  59                     "triSurface",                       // local
  60                     io.time(),                          // registry
  61                     IOobject::MUST_READ,
  62                     IOobject::NO_WRITE,
  63                     false
  64                 )
  65             )
  66         );
  67
  68         const featureEdgeMesh& eMesh = operator[](i);
  69
  70         //eMesh.mergePoints(meshRefiner_.mergeDistance());
  71         levels_[i] = readLabel(dict.lookup("level"));
  72
  73         Info<< "Refinement level " << levels_[i]
  74             << " for all cells crossed by feature " << featFileName
  75             << " (" << eMesh.points().size() << " points, "
  76             << eMesh.edges().size() << " edges)." << endl;
  77     }
  78
  79
  80     // Search engines
  81
  82     forAll(*this, i)
  83     {
  84         const featureEdgeMesh& eMesh = operator[](i);
  85         const pointField& points = eMesh.points();
  86         const edgeList& edges = eMesh.edges();
  87
  88         // Calculate bb of all points
  89         treeBoundBox bb(points);
  90
  91         // Random number generator. Bit dodgy since not exactly random ;-)
  92         Random rndGen(65431);
  93
  94         // Slightly extended bb. Slightly off-centred just so on symmetric
  95         // geometry there are less face/edge aligned items.
  96         bb = bb.extend(rndGen, 1E-4);
  97         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
  98         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
  99
 100         edgeTrees_.set
 101         (
 102             i,
 103             new indexedOctree<treeDataEdge>
 104             (
 105                 treeDataEdge
 106                 (
 107                     false,                  // do not cache bb
 108                     edges,
 109                     points,
 110                     identity(edges.size())
 111                 ),
 112                 bb,     // overall search domain
 113                 8,      // maxLevel
 114                 10,     // leafsize
 115                 3.0     // duplicity
 116             )
 117         );
 118
 119
 120         // Detect feature points from edges.
 121         const labelListList& pointEdges = eMesh.pointEdges();
 122         DynamicList<label> featurePoints;
 123         forAll(pointEdges, pointI)
 124         {
 125             if (pointEdges[pointI].size() > 2)
 126             {
 127                 featurePoints.append(pointI);
 128             }
 129         }
 130
 131         Info<< "Detected " << featurePoints.size()
 132             << " featurePoints out of " << points.size()
 133             << " on feature " << eMesh.name() << endl;
 134
 135         pointTrees_.set
 136         (
 137             i,
 138             new indexedOctree<treeDataPoint>
 139             (
 140                 treeDataPoint(points, featurePoints),
 141                 bb,     // overall search domain
 142                 8,      // maxLevel
 143                 10,     // leafsize
 144                 3.0     // duplicity
 145             )
 146         );
 147     }
 148 }
 149
 150
 151 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 152
 153 void Foam::refinementFeatures::findNearestEdge
 154 (
 155     const pointField& samples,
 156     const scalarField& nearestDistSqr,
 157     labelList& nearFeature,
 158     List<pointIndexHit>& nearInfo
 159 ) const
 160 {
 161     nearFeature.setSize(samples.size());
 162     nearFeature = -1;
 163     nearInfo.setSize(samples.size());
 164
 165     forAll(edgeTrees_, featI)
 166     {
 167         const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
 168
 169         if (tree.shapes().size() > 0)
 170         {
 171             forAll(samples, sampleI)
 172             {
 173                 const point& sample = samples[sampleI];
 174
 175                 scalar distSqr;
 176                 if (nearInfo[sampleI].hit())
 177                 {
 178                     distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
 179                 }
 180                 else
 181                 {
 182                     distSqr = nearestDistSqr[sampleI];
 183                 }
 184
 185                 pointIndexHit info = tree.findNearest(sample, distSqr);
 186
 187                 if (info.hit())
 188                 {
 189                     nearInfo[sampleI] = info;
 190                     nearFeature[sampleI] = featI;
 191                 }
 192             }
 193         }
 194     }
 195 }
 196
 197
 198 void Foam::refinementFeatures::findNearestPoint
 199 (
 200     const pointField& samples,
 201     const scalarField& nearestDistSqr,
 202     labelList& nearFeature,
 203     labelList& nearIndex
 204 ) const
 205 {
 206     nearFeature.setSize(samples.size());
 207     nearFeature = -1;
 208     nearIndex.setSize(samples.size());
 209     nearIndex = -1;
 210
 211     forAll(pointTrees_, featI)
 212     {
 213         const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
 214
 215         if (tree.shapes().pointLabels().size() > 0)
 216         {
 217             forAll(samples, sampleI)
 218             {
 219                 const point& sample = samples[sampleI];
 220
 221                 scalar distSqr;
 222                 if (nearFeature[sampleI] != -1)
 223                 {
 224                     label nearFeatI = nearFeature[sampleI];
 225                     const indexedOctree<treeDataPoint>& nearTree =
 226                         pointTrees_[nearFeatI];
 227                     label featPointI =
 228                         nearTree.shapes().pointLabels()[nearIndex[sampleI]];
 229                     const point& featPt =
 230                         operator[](nearFeatI).points()[featPointI];
 231                     distSqr = magSqr(featPt-sample);
 232                 }
 233                 else
 234                 {
 235                     distSqr = nearestDistSqr[sampleI];
 236                 }
 237
 238                 pointIndexHit info = tree.findNearest(sample, distSqr);
 239
 240                 if (info.hit())
 241                 {
 242                     nearFeature[sampleI] = featI;
 243                     nearIndex[sampleI] = info.index();
 244                 }
 245             }
 246         }
 247     }
 248 }
 249
 250
 251 // ************************************************************************* //