src/foam/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolationPatches.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     Mixing plane class dealing with transfer of data between two
  26     primitivePatches
  27
  28 Author
  29     Martin Beaudoin, Hydro-Quebec, 2009.  All rights reserved
  30
  31 Contributor
  32     Hrvoje Jasak, Wikki Ltd.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "MixingPlaneInterpolationTemplate.H"
  37 #include "demandDrivenData.H"
  38 #include "PrimitivePatchTemplate.H"
  39 #include "IOmanip.H"
  40 #include "transform.H"
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 //- Find and fix the patch faces that have cylindrical points coordinates in
  50 //  contact with or split on both sides of the -180/+180 degree
  51 //  axis eg: [-179,+179]
  52 //
  53 //  The cure is simple: we first find the problematic faces, (there might be
  54 //  none), then we modify only their angle coordinate in
  55 //  order to bring the whole face across the angle = +0 axis instead.
  56 //  This is a simple shift in cylindrical coordinates,
  57 //  but also a pure rotation in cartesian space for the faces.
  58 //
  59 //  Since the only objective for these localCoordFaces is to create a patch for
  60 //  constructing a temporary GGI in cylindrical coordinate
  61 //  space, a single shift will have no effect whatsoever on the GGI weights
  62 //   because we are using 360 degrees face ribbons for the other
  63 //  GGI patch.
  64 //
  65 //  By doing so, we will obvioulsy mess up the clean topology of the original
  66 //  cylindrical coordinates patch by shifting around a few faces;
  67 //  but keep in mind the only topology necessary for the GGI is the faces
  68 //  topology. A bunch of arbitrary located faces is as good as a bunch
  69 //  of cleanly laid-out faces into a regular patch. By only shifting the face
  70 //  along the angle axis, the face topology is preserved.
  71 //
  72 // Another interesting side-effect of this procedure is that it will cut
  73 //  a 360 degrees patch along the -180/+180 axis, leaving two un-connected
  74 // "right-side" and "left-side" edges when switching the patch points to
  75 //  cylindrical coordinates.
  76 template<class MasterPatch, class SlavePatch>
  77 void
  78 MixingPlaneInterpolation<MasterPatch, SlavePatch>::correctStraddlingFaces
  79 (
  80     faceList& localCoordFaces,
  81     pointField& localCoordFacesPoint
  82 ) const
  83 {
  84     // Memorize the list of displaced points so we can add them once at the
  85     // end. This minimizes the
  86     // memory reallocation necessary if we would do it one at a time instead.
  87
  88     typedef std::map<label, point> labelPointMap;
  89
  90     labelPointMap newFacePoints;
  91
  92     label newPointLabel = localCoordFacesPoint.size();
  93
  94     // Collect data on span limiting from coordinate system
  95     coordinateSystem::spanInfo spanLimited = cs_.spanLimited();
  96     boundBox spanBounds = cs_.spanBounds();
  97
  98     const direction sweepDir = sweepAxisSwitch();
  99
 100     // If there are no limits on both sides of the span,
 101     // adjustment can be skipped
 102     if
 103     (
 104         !spanLimited[sweepDir].first()
 105      || !spanLimited[sweepDir].second()
 106     )
 107     {
 108         // Nothing to do
 109         return;
 110     }
 111
 112     // Check if coordinate system is limited in spanwise range
 113     scalar maxSpan = spanBounds.max().component(sweepDir);
 114
 115     if (debug)
 116     {
 117         Info<< "Fixing span straddle in direction " << sweepDir
 118             << " for span " << maxSpan << endl;
 119     }
 120
 121     // Do all faces
 122     forAll (localCoordFaces, sFi)
 123     {
 124         bool faceIsOk(true);
 125
 126         scalarField pointsAngleValues =
 127             localCoordFaces[sFi].points
 128             (
 129                 localCoordFacesPoint
 130             ).component(sweepDir);
 131
 132         for (label i = 0; i < pointsAngleValues.size(); i++)
 133         {
 134             for (label j = i + 1; j < pointsAngleValues.size(); j++)
 135             {
 136                 // We shift away any faces in contact with the
 137                 // straddling axis
 138                 if
 139                 (
 140                     (
 141                         mag(pointsAngleValues[i] - pointsAngleValues[j])
 142                       > maxSpan
 143                     )
 144                  || (mag(mag(pointsAngleValues[i]) - maxSpan)) < VSMALL
 145                  || (mag(mag(pointsAngleValues[j]) - maxSpan)) < VSMALL
 146                 )
 147                 {
 148                     // We need to correct this
 149                     // Grab the original points labels
 150                     labelList pointLbls = localCoordFaces[sFi];
 151
 152                     forAll (pointLbls, ptI)
 153                     {
 154                         point ptCoord = localCoordFacesPoint[pointLbls[ptI]];
 155
 156                         // Switch point across the +0 angle axis
 157                         if (ptCoord[sweepDir] < 0.0)
 158                         {
 159                             ptCoord[sweepDir] = ptCoord[sweepDir] + maxSpan;
 160                         }
 161                         else
 162                         {
 163                             ptCoord[sweepDir] = ptCoord[sweepDir] - maxSpan;
 164                         }
 165
 166                         // Memorize data in order to reallocate
 167                         // localCoordFacesPoint only once
 168                         newFacePoints.insert
 169                         (
 170                             std::pair<label, point>(newPointLabel, ptCoord)
 171                         );
 172
 173                         // Memorize the new point label right away for this
 174                         // face and increment index for the next
 175                         localCoordFaces[sFi][ptI] = newPointLabel++;
 176                     }
 177                     // Done with this face
 178                     faceIsOk = false;
 179
 180                     break;
 181                 }
 182             }
 183
 184             if (!faceIsOk) break;
 185         }
 186     }
 187
 188     // Transfer the new points information into localCoordFacesPoint using
 189     // only one reallocation
 190     // We probably have point duplicates here; no big deal...
 191     localCoordFacesPoint.setSize
 192     (
 193         localCoordFacesPoint.size() + newFacePoints.size()
 194     );
 195
 196     forAllIter(labelPointMap, newFacePoints, nPi)
 197     {
 198           // First == label, Second == point
 199         localCoordFacesPoint[nPi->first] = nPi->second;
 200     }
 201 }
 202
 203
 204 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 205
 206 template<class MasterPatch, class SlavePatch>
 207 void
 208 MixingPlaneInterpolation<MasterPatch, SlavePatch>::
 209 calcTransformedPatches() const
 210 {
 211     if
 212     (
 213         transformedMasterPatchPtr_
 214      || transformedShadowPatchPtr_
 215     )
 216     {
 217         FatalErrorIn
 218         (
 219             "void MixingPlaneInterpolation::calcTransformedPatches() const"
 220         )   << "Patches already calculated"
 221             << abort(FatalError);
 222     }
 223
 224     // Duplicate the master/slave patch faces
 225     faceList masterFaces = masterPatch_.localFaces();
 226     faceList slaveFaces = slavePatch_.localFaces();
 227
 228     // Let's compute the patches coordinates values in mixing coordinates
 229     pointField masterPointsLocalCoord =
 230         cs_.localPosition(masterPatch_.localPoints());
 231
 232     pointField slavePointsLocalCoord =
 233         cs_.localPosition(slavePatch_.localPoints());
 234
 235     // Next, we need to find and fix the patch faces that have straddled
 236     // the span
 237     // Note: The face indices and the local coordinates and the
 238     // number of points may be modified within correctStraddlingFaces
 239     // MB, 27/Jan/2011
 240
 241     correctStraddlingFaces(masterFaces, masterPointsLocalCoord);
 242     correctStraddlingFaces(slaveFaces, slavePointsLocalCoord);
 243
 244     if(debug)
 245     {
 246         InfoIn
 247         (
 248             "MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 249             "calcTransformedPatches()"
 250         )   << "masterPointsLocalCoord: "
 251             << masterPointsLocalCoord << nl
 252             << "slavePointsLocalCoord: "
 253             << slavePointsLocalCoord << endl;
 254     }
 255
 256     // Create the local coords patches
 257
 258     transformedMasterPatchPtr_ =
 259         new standAlonePatch
 260         (
 261             masterFaces,
 262             masterPointsLocalCoord
 263         );
 264
 265     transformedShadowPatchPtr_ =
 266         new standAlonePatch
 267         (
 268             slaveFaces,
 269             slavePointsLocalCoord
 270         );
 271 }
 272
 273
 274 template<class MasterPatch, class SlavePatch>
 275 void
 276 MixingPlaneInterpolation<MasterPatch, SlavePatch>::calcMixingPlanePatch() const
 277 {
 278     if (mixingPlanePatchPtr_)
 279     {
 280         FatalErrorIn
 281         (
 282             "void MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 283             "calcMixingPlanePatch() const"
 284         )   << "Circumferential average patch alreacy calculated"
 285             << abort(FatalError);
 286     }
 287
 288     const pointField& profile = this->interpolationProfile();
 289
 290     // Create patch from profile
 291
 292     // Grab the bounding boxes of master and slave
 293
 294     boundBox masterBB
 295     (
 296         transformedMasterPatch().localPoints(),
 297         false
 298     );
 299
 300     boundBox slaveBB
 301     (
 302         transformedShadowPatch().localPoints(),
 303         false
 304     );
 305
 306     // Collect data on span limiting from coordinate system
 307     coordinateSystem::spanInfo spanLimited = cs_.spanLimited();
 308     boundBox spanBounds = cs_.spanBounds();
 309
 310     const direction sweepDir = sweepAxisSwitch();
 311
 312     // Get span from bounding boxes
 313     scalar minSpan =
 314         Foam::min
 315         (
 316             masterBB.min()[sweepDir],
 317             slaveBB.min()[sweepDir]
 318         ) - SMALL;
 319
 320     scalar maxSpan =
 321         Foam::max
 322         (
 323             masterBB.max()[sweepDir],
 324             slaveBB.max()[sweepDir]
 325         ) + SMALL;
 326
 327     if (debug)
 328     {
 329         InfoIn
 330         (
 331             "MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 332             "calcMixingPlanePatch() const"
 333         )   << "minSpan from patch BB : " << minSpan << nl
 334             << "maxSpan from patch BB : " << maxSpan << endl;
 335     }
 336
 337     // Correct for limited span
 338     if (spanLimited[sweepDir].first())
 339     {
 340         minSpan = spanBounds.min()[sweepDir];
 341     }
 342
 343     if (spanLimited[sweepDir].second())
 344     {
 345         maxSpan = spanBounds.max()[sweepDir];
 346     }
 347
 348     if (debug)
 349     {
 350         InfoIn
 351         (
 352             "MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 353             "calcMixingPlanePatch() const"
 354         )   << "minSpan after checking spanLimited: " << minSpan << nl
 355             << "maxSpan after checking spanLimited: " << maxSpan << endl;
 356     }
 357
 358     label nRibbons = profile.size() - 1;
 359
 360     // Points for averaging patch in the local coordinate system
 361     pointField mixingPatchPoints(2*nRibbons + 2);
 362     label nextPointID = 0;
 363
 364     faceList mixingPatchFaces(nRibbons);
 365     scalarField mixingPatchFacesArea(nRibbons);
 366
 367     // Insert lower bound points and expand to bounds
 368     mixingPatchPoints[0] = profile[0];
 369     mixingPatchPoints[0].replace(sweepDir, minSpan);
 370
 371     mixingPatchPoints[1] = profile[0];
 372     mixingPatchPoints[1].replace(sweepDir, maxSpan);
 373     nextPointID = 2;
 374
 375
 376     // Create patch faces
 377     forAll (mixingPatchFaces, fI)
 378     {
 379         label fi2 = 2*fI;
 380         // Add top bound points and expand to bounds
 381         mixingPatchPoints[nextPointID] = profile[fI + 1];
 382         mixingPatchPoints[nextPointID].replace(sweepDir, minSpan);
 383         nextPointID++;
 384
 385         mixingPatchPoints[nextPointID] = profile[fI + 1];
 386         mixingPatchPoints[nextPointID].replace(sweepDir, maxSpan);
 387         nextPointID++;
 388
 389         face curFace(4);
 390         curFace[0] = fi2;
 391         curFace[1] = fi2 + 1;
 392         curFace[2] = fi2 + 3;
 393         curFace[3] = fi2 + 2;
 394
 395         mixingPatchFaces[fI] = curFace;
 396     }
 397
 398     mixingPlanePatchPtr_ =
 399         new standAlonePatch
 400         (
 401             mixingPatchFaces,
 402             mixingPatchPoints
 403         );
 404
 405
 406     if (debug > 0)
 407     {
 408         InfoIn
 409         (
 410             "MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
 411             "calcMixingPlanePatch() const"
 412         )   << "mixingPatch: "
 413             << *mixingPlanePatchPtr_ << nl
 414             << "mixingPatch.points : "
 415             << mixingPlanePatchPtr_->points() << endl;
 416     }
 417 }
 418
 419
 420 template<class MasterPatch, class SlavePatch>
 421 void
 422 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearTransfomedPatches()
 423 {
 424     deleteDemandDrivenData(transformedMasterPatchPtr_);
 425     deleteDemandDrivenData(transformedShadowPatchPtr_);
 426 }
 427
 428
 429 template<class MasterPatch, class SlavePatch>
 430 void
 431 MixingPlaneInterpolation<MasterPatch, SlavePatch>::clearMixingPlanePatch()
 432 {
 433     deleteDemandDrivenData(mixingPlanePatchPtr_);
 434 }
 435
 436
 437 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 438
 439 } // End namespace Foam
 440
 441 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 442
 443
 444 // ************************************************************************* //