src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "surfaceSlipDisplacementPointPatchVectorField.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "Time.H"
  30 #include "transformField.H"
  31 #include "fvMesh.H"
  32 #include "displacementFvMotionSolver.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 template<>
  42 const char*
  43 NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>::
  44 names[] =
  45 {
  46     "nearest",
  47     "pointNormal",
  48     "fixedNormal"
  49 };
  50
  51 const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>
  52     surfaceSlipDisplacementPointPatchVectorField::projectModeNames_;
  53
  54
  55 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  56
  57 void surfaceSlipDisplacementPointPatchVectorField::calcProjection
  58 (
  59     vectorField& displacement
  60 ) const
  61 {
  62     const polyMesh& mesh = patch().boundaryMesh().mesh()();
  63     const pointField& localPoints = patch().localPoints();
  64     const labelList& meshPoints = patch().meshPoints();
  65
  66     //const scalar deltaT = mesh.time().deltaT().value();
  67
  68     // Construct large enough vector in direction of projectDir so
  69     // we're guaranteed to hit something.
  70
  71     //- Per point projection vector:
  72     const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
  73
  74     // For case of fixed projection vector:
  75     vector projectVec;
  76     if (projectMode_ == FIXEDNORMAL)
  77     {
  78         vector n = projectDir_/mag(projectDir_);
  79         projectVec = projectLen*n;
  80     }
  81
  82
  83     // Get fixed points (bit of a hack)
  84     const pointZone* zonePtr = NULL;
  85
  86     if (frozenPointsZone_.size() > 0)
  87     {
  88         const pointZoneMesh& pZones = mesh.pointZones();
  89
  90         zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)];
  91
  92         Pout<< "surfaceSlipDisplacementPointPatchVectorField : Fixing all "
  93             << zonePtr->size() << " points in pointZone " << zonePtr->name()
  94             << endl;
  95     }
  96
  97     // Get the starting locations from the motionSolver
  98     const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
  99     (
 100         "dynamicMeshDict"
 101     ).points0();
 102
 103
 104     pointField start(meshPoints.size());
 105     forAll(start, i)
 106     {
 107         start[i] = points0[meshPoints[i]] + displacement[i];
 108     }
 109
 110     label nNotProjected = 0;
 111
 112     if (projectMode_ == NEAREST)
 113     {
 114         List<pointIndexHit> nearest;
 115         labelList hitSurfaces;
 116         surfaces().findNearest
 117         (
 118             start,
 119             scalarField(start.size(), sqr(projectLen)),
 120             hitSurfaces,
 121             nearest
 122         );
 123
 124         forAll(nearest, i)
 125         {
 126             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 127             {
 128                 // Fixed point. Reset to point0 location.
 129                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 130             }
 131             else if (nearest[i].hit())
 132             {
 133                 displacement[i] =
 134                     nearest[i].hitPoint()
 135                   - points0[meshPoints[i]];
 136             }
 137             else
 138             {
 139                 nNotProjected++;
 140
 141                 if (debug)
 142                 {
 143                     Pout<< "    point:" << meshPoints[i]
 144                         << " coord:" << localPoints[i]
 145                         << "  did not find any surface within " << projectLen
 146                         << endl;
 147                 }
 148             }
 149         }
 150     }
 151     else
 152     {
 153         // Do tests on all points. Combine later on.
 154
 155         // 1. Check if already on surface
 156         List<pointIndexHit> nearest;
 157         {
 158             labelList nearestSurface;
 159             surfaces().findNearest
 160             (
 161                 start,
 162                 scalarField(start.size(), sqr(SMALL)),
 163                 nearestSurface,
 164                 nearest
 165             );
 166         }
 167
 168         // 2. intersection. (combined later on with information from nearest
 169         // above)
 170         vectorField projectVecs(start.size(), projectVec);
 171
 172         if (projectMode_ == POINTNORMAL)
 173         {
 174             projectVecs = projectLen*patch().pointNormals();
 175         }
 176
 177         // Knock out any wedge component
 178         scalarField offset(start.size(), 0.0);
 179         if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 180         {
 181             forAll(offset, i)
 182             {
 183                 offset[i] = start[i][wedgePlane_];
 184                 start[i][wedgePlane_] = 0;
 185                 projectVecs[i][wedgePlane_] = 0;
 186             }
 187         }
 188
 189         List<pointIndexHit> rightHit;
 190         {
 191             labelList rightSurf;
 192             surfaces().findAnyIntersection
 193             (
 194                 start,
 195                 start+projectVecs,
 196                 rightSurf,
 197                 rightHit
 198             );
 199         }
 200
 201         List<pointIndexHit> leftHit;
 202         {
 203             labelList leftSurf;
 204             surfaces().findAnyIntersection
 205             (
 206                 start,
 207                 start-projectVecs,
 208                 leftSurf,
 209                 leftHit
 210             );
 211         }
 212
 213         // 3. Choose either -fixed, nearest, right, left.
 214         forAll(displacement, i)
 215         {
 216             if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
 217             {
 218                 // Fixed point. Reset to point0 location.
 219                 displacement[i] = points0[meshPoints[i]] - localPoints[i];
 220             }
 221             else if (nearest[i].hit())
 222             {
 223                 // Found nearest.
 224                 displacement[i] =
 225                     nearest[i].hitPoint()
 226                   - points0[meshPoints[i]];
 227             }
 228             else
 229             {
 230                 pointIndexHit interPt;
 231
 232                 if (rightHit[i].hit())
 233                 {
 234                     if (leftHit[i].hit())
 235                     {
 236                         if
 237                         (
 238                             magSqr(rightHit[i].hitPoint()-start[i])
 239                           < magSqr(leftHit[i].hitPoint()-start[i])
 240                         )
 241                         {
 242                             interPt = rightHit[i];
 243                         }
 244                         else
 245                         {
 246                             interPt = leftHit[i];
 247                         }
 248                     }
 249                     else
 250                     {
 251                         interPt = rightHit[i];
 252                     }
 253                 }
 254                 else
 255                 {
 256                     if (leftHit[i].hit())
 257                     {
 258                         interPt = leftHit[i];
 259                     }
 260                 }
 261
 262
 263                 if (interPt.hit())
 264                 {
 265                     if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
 266                     {
 267                         interPt.rawPoint()[wedgePlane_] += offset[i];
 268                     }
 269                     displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
 270                 }
 271                 else
 272                 {
 273                     nNotProjected++;
 274
 275                     if (debug)
 276                     {
 277                         Pout<< "    point:" << meshPoints[i]
 278                             << " coord:" << localPoints[i]
 279                             << "  did not find any intersection between"
 280                             << " ray from " << start[i]-projectVecs[i]
 281                             << " to " << start[i]+projectVecs[i] << endl;
 282                     }
 283                 }
 284             }
 285         }
 286     }
 287
 288     reduce(nNotProjected, sumOp<label>());
 289
 290     if (nNotProjected > 0)
 291     {
 292         Info<< "surfaceSlipDisplacement :"
 293             << " on patch " << patch().name()
 294             << " did not project " << nNotProjected
 295             << " out of " << returnReduce(localPoints.size(), sumOp<label>())
 296             << " points." << endl;
 297     }
 298 }
 299
 300
 301 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 302
 303 surfaceSlipDisplacementPointPatchVectorField::
 304 surfaceSlipDisplacementPointPatchVectorField
 305 (
 306     const pointPatch& p,
 307     const DimensionedField<vector, pointMesh>& iF
 308 )
 309 :
 310     pointPatchVectorField(p, iF),
 311     projectMode_(NEAREST),
 312     projectDir_(vector::zero),
 313     wedgePlane_(-1)
 314 {}
 315
 316
 317 surfaceSlipDisplacementPointPatchVectorField::
 318 surfaceSlipDisplacementPointPatchVectorField
 319 (
 320     const pointPatch& p,
 321     const DimensionedField<vector, pointMesh>& iF,
 322     const dictionary& dict
 323 )
 324 :
 325     pointPatchVectorField(p, iF, dict),
 326     surfacesDict_(dict.subDict("geometry")),
 327     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
 328     projectDir_(dict.lookup("projectDirection")),
 329     wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
 330     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 331 {}
 332
 333
 334 surfaceSlipDisplacementPointPatchVectorField::
 335 surfaceSlipDisplacementPointPatchVectorField
 336 (
 337     const surfaceSlipDisplacementPointPatchVectorField& ppf,
 338     const pointPatch& p,
 339     const DimensionedField<vector, pointMesh>& iF,
 340     const pointPatchFieldMapper&
 341 )
 342 :
 343     pointPatchVectorField(p, iF),
 344     surfacesDict_(ppf.surfacesDict_),
 345     projectMode_(ppf.projectMode_),
 346     projectDir_(ppf.projectDir_),
 347     wedgePlane_(ppf.wedgePlane_),
 348     frozenPointsZone_(ppf.frozenPointsZone_)
 349 {}
 350
 351
 352 surfaceSlipDisplacementPointPatchVectorField::
 353 surfaceSlipDisplacementPointPatchVectorField
 354 (
 355     const surfaceSlipDisplacementPointPatchVectorField& ppf
 356 )
 357 :
 358     pointPatchVectorField(ppf),
 359     surfacesDict_(ppf.surfacesDict_),
 360     projectMode_(ppf.projectMode_),
 361     projectDir_(ppf.projectDir_),
 362     wedgePlane_(ppf.wedgePlane_),
 363     frozenPointsZone_(ppf.frozenPointsZone_)
 364 {}
 365
 366
 367 surfaceSlipDisplacementPointPatchVectorField::
 368 surfaceSlipDisplacementPointPatchVectorField
 369 (
 370     const surfaceSlipDisplacementPointPatchVectorField& ppf,
 371     const DimensionedField<vector, pointMesh>& iF
 372 )
 373 :
 374     pointPatchVectorField(ppf, iF),
 375     surfacesDict_(ppf.surfacesDict_),
 376     projectMode_(ppf.projectMode_),
 377     projectDir_(ppf.projectDir_),
 378     wedgePlane_(ppf.wedgePlane_),
 379     frozenPointsZone_(ppf.frozenPointsZone_)
 380 {}
 381
 382
 383 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 384
 385 const searchableSurfaces&
 386 surfaceSlipDisplacementPointPatchVectorField::surfaces() const
 387 {
 388     if (surfacesPtr_.empty())
 389     {
 390         surfacesPtr_.reset
 391         (
 392             new searchableSurfaces
 393             (
 394                 IOobject
 395                 (
 396                     "abc",                              // dummy name
 397                     db().time().constant(),             // directory
 398                     "triSurface",                       // instance
 399                     db().time(),                        // registry
 400                     IOobject::MUST_READ,
 401                     IOobject::NO_WRITE
 402                 ),
 403                 surfacesDict_
 404             )
 405         );
 406     }
 407     return surfacesPtr_();
 408 }
 409
 410
 411 void surfaceSlipDisplacementPointPatchVectorField::evaluate
 412 (
 413     const Pstream::commsTypes commsType
 414 )
 415 {
 416     vectorField displacement(this->patchInternalField());
 417
 418     // Calculate displacement to project points onto surface
 419     calcProjection(displacement);
 420
 421     // Get internal field to insert values into
 422     Field<vector>& iF = const_cast<Field<vector>&>(this->internalField());
 423
 424     //setInInternalField(iF, motionU);
 425     setInInternalField(iF, displacement);
 426
 427     pointPatchVectorField::evaluate(commsType);
 428 }
 429
 430
 431 void surfaceSlipDisplacementPointPatchVectorField::write(Ostream& os) const
 432 {
 433     pointPatchVectorField::write(os);
 434     os.writeKeyword("geometry") << surfacesDict_
 435         << token::END_STATEMENT << nl;
 436     os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
 437         << token::END_STATEMENT << nl;
 438     os.writeKeyword("projectDirection") << projectDir_
 439         << token::END_STATEMENT << nl;
 440     os.writeKeyword("wedgePlane") << wedgePlane_
 441         << token::END_STATEMENT << nl;
 442     if (frozenPointsZone_ != word::null)
 443     {
 444         os.writeKeyword("frozenPointsZone") << frozenPointsZone_
 445             << token::END_STATEMENT << nl;
 446     }
 447 }
 448
 449
 450 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 451
 452 makePointPatchTypeField
 453 (
 454     pointPatchVectorField,
 455     surfaceSlipDisplacementPointPatchVectorField
 456 );
 457
 458
 459 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 460
 461 } // End namespace Foam
 462
 463 // ************************************************************************* //