src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C

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