src/dynamicMesh/meshMotion/RBFMotionSolver/RBFMotionSolver.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "RBFMotionSolver.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(RBFMotionSolver, 0);
  34
  35     addToRunTimeSelectionTable
  36     (
  37         motionSolver,
  38         RBFMotionSolver,
  39         dictionary
  40     );
  41 }
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 void Foam::RBFMotionSolver::makeControlIDs()
  46 {
  47     // Points that are neither on moving nor on static patches
  48     // will be marked with 0
  49     labelList markedPoints(mesh().nPoints(), 0);
  50
  51     // Mark all points on moving patches with 1
  52     label nMovingPoints = 0;
  53
  54     forAll (movingPatches_, patchI)
  55     {
  56         // Find the patch in boundary
  57         label patchIndex =
  58             mesh().boundaryMesh().findPatchID(movingPatches_[patchI]);
  59
  60         if (patchIndex < 0)
  61         {
  62             FatalErrorIn("void RBFMotionSolver::makeControlIDs()")
  63                 << "Patch " << movingPatches_[patchI] << " not found.  "
  64                 << "valid patch names: " << mesh().boundaryMesh().names()
  65                 << abort(FatalError);
  66         }
  67
  68         const labelList& mp = mesh().boundaryMesh()[patchIndex].meshPoints();
  69
  70         forAll (mp, i)
  71         {
  72             markedPoints[mp[i]] = 1;
  73             nMovingPoints++;
  74         }
  75     }
  76
  77     // Mark moving points and select control points from moving patches
  78     movingIDs_.setSize(nMovingPoints);
  79
  80     Info<< "Total points on moving boundaries: " << nMovingPoints << endl;
  81
  82     const pointField& points = mesh().points();
  83
  84     // Re-use counter to count moving points
  85     // Note: the control points also hold static points in the second part
  86     // of the list if static patches are included in the RBF
  87     // HJ, 24/Mar/2011
  88     nMovingPoints = 0;
  89
  90     // Count moving points first
  91     forAll (markedPoints, i)
  92     {
  93         if (markedPoints[i] == 1)
  94         {
  95             // Grab internal point
  96             movingIDs_[nMovingPoints] = i;
  97             nMovingPoints++;
  98         }
  99     }
 100
 101     movingIDs_.setSize(nMovingPoints);
 102
 103     // Actual location of moving points will be set later on request
 104     // HJ, 19/Dec/2008
 105     movingPoints_.setSize(nMovingPoints, vector::zero);
 106
 107     // Mark all points on static patches with -1
 108     label nStaticPoints = 0;
 109
 110     forAll (staticPatches_, patchI)
 111     {
 112         // Find the patch in boundary
 113         label patchIndex =
 114             mesh().boundaryMesh().findPatchID(staticPatches_[patchI]);
 115
 116         if (patchIndex < 0)
 117         {
 118             FatalErrorIn("void RBFMotionSolver::makeControlPoints()")
 119                 << "Patch " << staticPatches_[patchI] << " not found.  "
 120                 << "valid patch names: " << mesh().boundaryMesh().names()
 121                 << abort(FatalError);
 122         }
 123
 124         const labelList& mp = mesh().boundaryMesh()[patchIndex].meshPoints();
 125
 126         forAll (mp, i)
 127         {
 128             markedPoints[mp[i]] = -1;
 129             nStaticPoints++;
 130         }
 131     }
 132
 133     Info<< "Total points on static boundaries: " << nStaticPoints << endl;
 134     staticIDs_.setSize(nStaticPoints);
 135
 136     // Re-use counter
 137     nStaticPoints = 0;
 138
 139     // Count total number of control points
 140     forAll (markedPoints, i)
 141     {
 142         if (markedPoints[i] == -1)
 143         {
 144             staticIDs_[nStaticPoints] = i;
 145             nStaticPoints++;
 146         }
 147     }
 148
 149     staticIDs_.setSize(nStaticPoints);
 150
 151     // Control IDs also potentially include points on static patches
 152     // HJ, 24/Mar/2011
 153     controlIDs_.setSize(movingIDs_.size() + staticIDs_.size());
 154     motion_.setSize(controlIDs_.size(), vector::zero);
 155
 156     label nControlPoints = 0;
 157
 158     forAll (movingPatches_, patchI)
 159     {
 160         // Find the patch in boundary
 161         label patchIndex =
 162             mesh().boundaryMesh().findPatchID(movingPatches_[patchI]);
 163
 164         const labelList& mp = mesh().boundaryMesh()[patchIndex].meshPoints();
 165
 166         for
 167         (
 168             label pickedPoint = 0;
 169             pickedPoint < mp.size();
 170             pickedPoint += coarseningRatio_
 171         )
 172         {
 173             // Pick point as control point
 174             controlIDs_[nControlPoints] = mp[pickedPoint];
 175
 176             // Mark the point as picked
 177             markedPoints[mp[pickedPoint]] = 2;
 178             nControlPoints++;
 179         }
 180     }
 181
 182     Info<< "Selected " << nControlPoints
 183         << " control points on moving boundaries" << endl;
 184
 185     if (includeStaticPatches_)
 186     {
 187         forAll (staticPatches_, patchI)
 188         {
 189             // Find the patch in boundary
 190             label patchIndex =
 191                 mesh().boundaryMesh().findPatchID(staticPatches_[patchI]);
 192
 193             const labelList& mp =
 194                 mesh().boundaryMesh()[patchIndex].meshPoints();
 195
 196             for
 197             (
 198                 label pickedPoint = 0;
 199                 pickedPoint < mp.size();
 200                 pickedPoint += coarseningRatio_
 201             )
 202             {
 203                 // Pick point as control point
 204                 controlIDs_[nControlPoints] = mp[pickedPoint];
 205
 206                 // Mark the point as picked
 207                 markedPoints[mp[pickedPoint]] = 2;
 208                 nControlPoints++;
 209             }
 210         }
 211
 212         Info<< "Selected " << nControlPoints
 213             << " total control points" << endl;
 214     }
 215
 216     // Resize control IDs
 217     controlIDs_.setSize(nControlPoints);
 218
 219     // Pick up point locations
 220     controlPoints_.setSize(nControlPoints);
 221
 222     // Set control points
 223     forAll (controlIDs_, i)
 224     {
 225         controlPoints_[i] = points[controlIDs_[i]];
 226     }
 227
 228     // Pick up all internal points
 229     internalIDs_.setSize(points.size());
 230     internalPoints_.setSize(points.size());
 231
 232     // Count internal points
 233     label nInternalPoints = 0;
 234
 235     forAll (markedPoints, i)
 236     {
 237         if (markedPoints[i] == 0)
 238         {
 239             // Grab internal point
 240             internalIDs_[nInternalPoints] = i;
 241             internalPoints_[nInternalPoints] = points[i];
 242             nInternalPoints++;
 243         }
 244     }
 245
 246     Info << "Number of internal points: " << nInternalPoints << endl;
 247
 248     // Resize the lists
 249     internalIDs_.setSize(nInternalPoints);
 250     internalPoints_.setSize(nInternalPoints);
 251 }
 252
 253
 254 void Foam::RBFMotionSolver::setMovingPoints() const
 255 {
 256     const pointField& points = mesh().points();
 257
 258     // Set moving points
 259     forAll (movingIDs_, i)
 260     {
 261         movingPoints_[i] = points[movingIDs_[i]];
 262     }
 263 }
 264
 265
 266 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 267
 268 Foam::RBFMotionSolver::RBFMotionSolver
 269 (
 270     const polyMesh& mesh,
 271     Istream&
 272 )
 273 :
 274     motionSolver(mesh),
 275     movingPatches_(lookup("movingPatches")),
 276     staticPatches_(lookup("staticPatches")),
 277     coarseningRatio_(readLabel(lookup("coarseningRatio"))),
 278     includeStaticPatches_(lookup("includeStaticPatches")),
 279     frozenInterpolation_(lookup("frozenInterpolation")),
 280     movingIDs_(0),
 281     movingPoints_(0),
 282     staticIDs_(0),
 283     controlIDs_(0),
 284     controlPoints_(0),
 285     internalIDs_(0),
 286     internalPoints_(0),
 287     motion_(0),
 288     interpolation_
 289     (
 290         subDict("interpolation"),
 291         controlPoints_,
 292         internalPoints_
 293     )
 294 {
 295     makeControlIDs();
 296 }
 297
 298
 299 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 300
 301 Foam::RBFMotionSolver::~RBFMotionSolver()
 302 {}
 303
 304
 305 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 306
 307 void Foam::RBFMotionSolver::setMotion(const vectorField& m)
 308 {
 309     if (m.size() != movingIDs_.size())
 310     {
 311         FatalErrorIn
 312         (
 313             "void RBFMotionSolver::setMotion(const vectorField& m)"
 314         )   << "Incorrect size of motion points: m = " << m.size()
 315             << " movingIDs = " << movingIDs_.size()
 316             << abort(FatalError);
 317     }
 318
 319     // Motion of static points is zero and moving points are first
 320     // in the list.  HJ, 24/Mar/2011
 321     motion_ = vector::zero;
 322
 323     forAll (m, i)
 324     {
 325         motion_[i] = m[i];
 326     }
 327
 328     if (!frozenInterpolation_)
 329     {
 330         // Set control points
 331         const pointField& points = mesh().points();
 332
 333         forAll (controlIDs_, i)
 334         {
 335             controlPoints_[i] = points[controlIDs_[i]];
 336         }
 337
 338         // Re-calculate interpolation
 339         interpolation_.movePoints();
 340     }
 341 }
 342
 343
 344 const Foam::vectorField& Foam::RBFMotionSolver::movingPoints() const
 345 {
 346     // Update moving points based on current mesh
 347     setMovingPoints();
 348
 349     return movingPoints_;
 350 }
 351
 352
 353 Foam::tmp<Foam::pointField> Foam::RBFMotionSolver::curPoints() const
 354 {
 355     // Prepare new points: same as old point
 356     tmp<pointField> tcurPoints
 357     (
 358         new vectorField(mesh().nPoints(), vector::zero)
 359     );
 360     pointField& curPoints = tcurPoints();
 361
 362     // Add motion to existing points
 363
 364     // 1. Insert prescribed motion of moving points
 365     forAll (movingIDs_, i)
 366     {
 367         curPoints[movingIDs_[i]] = motion_[i];
 368     }
 369
 370     // 2. Insert zero motion of static points
 371     forAll (staticIDs_, i)
 372     {
 373         curPoints[staticIDs_[i]] = vector::zero;
 374     }
 375
 376     // Set motion of control
 377     vectorField motionOfControl(controlIDs_.size());
 378
 379     // 2. Capture positions of control points
 380     forAll (controlIDs_, i)
 381     {
 382         motionOfControl[i] = curPoints[controlIDs_[i]];
 383     }
 384
 385     // Call interpolation
 386     vectorField interpolatedMotion =
 387         interpolation_.interpolate(motionOfControl);
 388
 389     // 3. Insert RBF interpolated motion
 390     forAll (internalIDs_, i)
 391     {
 392         curPoints[internalIDs_[i]] = interpolatedMotion[i];
 393     }
 394
 395     // 4. Add old point positions
 396     curPoints += mesh().points();
 397
 398     twoDCorrectPoints(tcurPoints());
 399
 400     return tcurPoints;
 401 }
 402
 403
 404 void Foam::RBFMotionSolver::solve()
 405 {}
 406
 407
 408 void Foam::RBFMotionSolver::updateMesh(const mapPolyMesh&)
 409 {
 410     // Recalculate control point IDs
 411     makeControlIDs();
 412 }
 413
 414
 415 // ************************************************************************* //