src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchFaces.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "enrichedPatch.H"
  27 #include "DynamicList.H"
  28
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 const Foam::label Foam::enrichedPatch::enrichedFaceRatio_ = 3;
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 void Foam::enrichedPatch::calcEnrichedFaces
  37 (
  38     const labelListList& pointsIntoMasterEdges,
  39     const labelListList& pointsIntoSlaveEdges,
  40     const pointField& projectedSlavePoints
  41 )
  42 {
  43     if (enrichedFacesPtr_)
  44     {
  45         FatalErrorIn
  46         (
  47             "void enrichedPatch::calcEnrichedFaces\n"
  48             "(\n"
  49             "    const labelListList& pointsIntoMasterEdges,\n"
  50             "    const labelListList& pointsIntoSlaveEdges,\n"
  51             "    const pointField& projectedSlavePoints\n"
  52             ")"
  53         )   << "Enriched faces already calculated."
  54             << abort(FatalError);
  55     }
  56
  57     // Create a list of enriched faces
  58     // Algorithm:
  59     // 1) Grab the original face and start from point zero.
  60     // 2) If the point has been merged away, grab the merge label;
  61     //    otherwise, keep the original label.
  62     // 3) Go to the next edge. Collect all the points to be added along
  63     //    the edge; order them in the necessary direction and insert onto the
  64     //    face.
  65     // 4) Grab the next point and return on step 2.
  66     enrichedFacesPtr_ = new faceList(masterPatch_.size() + slavePatch_.size());
  67     faceList& enrichedFaces = *enrichedFacesPtr_;
  68
  69     label nEnrichedFaces = 0;
  70
  71     const pointField& masterLocalPoints = masterPatch_.localPoints();
  72     const faceList& masterLocalFaces = masterPatch_.localFaces();
  73     const labelListList& masterFaceEdges = masterPatch_.faceEdges();
  74
  75     const faceList& slaveLocalFaces = slavePatch_.localFaces();
  76     const labelListList& slaveFaceEdges = slavePatch_.faceEdges();
  77
  78     // For correct functioning of the enrichedPatch class, the slave
  79     // faces need to be inserted first.  See comments in
  80     // enrichedPatch.H
  81
  82     // Get reference to the point merge map
  83     const Map<label>& pmm = pointMergeMap();
  84
  85     // Add slave faces into the enriched faces list
  86
  87     forAll(slavePatch_, faceI)
  88     {
  89         const face oldFace = slavePatch_[faceI];
  90         const face oldLocalFace = slaveLocalFaces[faceI];
  91 //         Info<< "old slave face " << faceI << ": " << oldFace << endl;
  92         const labelList& curEdges = slaveFaceEdges[faceI];
  93
  94         DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
  95
  96         // Note: The number of points and edges in a face is always identical
  97         // so both can be done is the same loop
  98         forAll(oldFace, i)
  99         {
 100             // Add the point
 101             Map<label>::const_iterator mpIter =
 102                 pmm.find(oldFace[i]);
 103
 104             if (mpIter == pmm.end())
 105             {
 106                 // Point not mapped
 107                 newFace.append(oldFace[i]);
 108
 109                 // Add the projected point into the patch support
 110                 pointMap().insert
 111                 (
 112                     oldFace[i],    // Global label of point
 113                     projectedSlavePoints[oldLocalFace[i]] // Projected position
 114                 );
 115             }
 116             else
 117             {
 118                 // Point mapped
 119                 newFace.append(mpIter());
 120
 121                 // Add the projected point into the patch support
 122                 pointMap().insert
 123                 (
 124                     mpIter(),    // Merged global label of point
 125                     projectedSlavePoints[oldLocalFace[i]] // Projected position
 126                 );
 127             }
 128
 129             // Grab the edge points
 130
 131             const labelList& slavePointsOnEdge =
 132                 pointsIntoSlaveEdges[curEdges[i]];
 133
 134             // Info<< "slavePointsOnEdge for "
 135             //     << curEdges[i] << ": " << slavePointsOnEdge
 136             //     << endl;
 137
 138             // If there are no points on the edge, skip everything
 139             // If there is only one point, no need for sorting
 140             if (slavePointsOnEdge.size())
 141             {
 142                 // Sort edge points in order
 143                 scalarField edgePointWeights(slavePointsOnEdge.size());
 144                 const point& startPoint = projectedSlavePoints[oldLocalFace[i]];
 145
 146                 vector e =
 147                     projectedSlavePoints[oldLocalFace.nextLabel(i)]
 148                   - startPoint;
 149
 150                 scalar magSqrE = magSqr(e);
 151
 152                 if (magSqrE > SMALL)
 153                 {
 154                     e /= magSqrE;
 155                 }
 156                 else
 157                 {
 158                     FatalErrorIn
 159                     (
 160                         "void enrichedPatch::calcEnrichedFaces\n"
 161                         "(\n"
 162                         "    const labelListList& pointsIntoMasterEdges,\n"
 163                         "    const labelListList& pointsIntoSlaveEdges,\n"
 164                         "    const pointField& projectedSlavePoints\n"
 165                         ")"
 166                     )   << "Zero length edge in slave patch for face " << i
 167                         << ".  This is not allowed."
 168                         << abort(FatalError);
 169                 }
 170
 171                 pointField slavePosOnEdge(slavePointsOnEdge.size());
 172
 173                 forAll(slavePointsOnEdge, edgePointI)
 174                 {
 175                     slavePosOnEdge[edgePointI] =
 176                         pointMap().find(slavePointsOnEdge[edgePointI])();
 177
 178                     edgePointWeights[edgePointI] =
 179                         (e & (slavePosOnEdge[edgePointI] - startPoint));
 180                 }
 181
 182                 if (debug)
 183                 {
 184                     // Check weights: all new points should be on the edge
 185                     if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
 186                     {
 187                         FatalErrorIn
 188                         (
 189                             "void enrichedPatch::calcEnrichedFaces\n"
 190                             "(\n"
 191                             "    const labelListList& pointsIntoMasterEdges,\n"
 192                             "    const labelListList& pointsIntoSlaveEdges,\n"
 193                             "    const pointField& projectedSlavePoints\n"
 194                             ")"
 195                         )   << "Invalid point edge weights.  Some of points are"
 196                             << " not on the edge for edge " << curEdges[i]
 197                             << " of face " << faceI << " in slave patch." << nl
 198                             << "Min weight: " << min(edgePointWeights)
 199                             << " Max weight: " << max(edgePointWeights)
 200                             << abort(FatalError);
 201                     }
 202                 }
 203
 204                 // Go through the points and collect them based on
 205                 // weights from lower to higher.  This gives the
 206                 // correct order of points along the edge.
 207                 forAll(edgePointWeights, passI)
 208                 {
 209                     // Max weight can only be one, so the sorting is
 210                     // done by elimination.
 211                     label nextPoint = -1;
 212                     scalar dist = 2;
 213
 214                     forAll(edgePointWeights, wI)
 215                     {
 216                         if (edgePointWeights[wI] < dist)
 217                         {
 218                             dist = edgePointWeights[wI];
 219                             nextPoint = wI;
 220                         }
 221                     }
 222
 223                     // Insert the next point and reset its weight to exclude it
 224                     // from future picks
 225                     newFace.append(slavePointsOnEdge[nextPoint]);
 226                     edgePointWeights[nextPoint] = GREAT;
 227
 228                     // Add the point into patch support
 229                     pointMap().insert
 230                     (
 231                         slavePointsOnEdge[nextPoint],
 232                         slavePosOnEdge[nextPoint]
 233                     );
 234                 }
 235             }
 236         }
 237         // Info<< "New slave face " << faceI << ": " << newFace << endl;
 238
 239         // Add the new face to the list
 240         enrichedFaces[nEnrichedFaces].transfer(newFace);
 241         nEnrichedFaces++;
 242     }
 243
 244     // Add master faces into the enriched faces list
 245
 246     forAll(masterPatch_, faceI)
 247     {
 248         const face& oldFace = masterPatch_[faceI];
 249         const face& oldLocalFace = masterLocalFaces[faceI];
 250 //         Info<< "old master face: " << oldFace << endl;
 251         const labelList& curEdges = masterFaceEdges[faceI];
 252
 253         DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
 254
 255         // Note: The number of points and edges in a face is always identical
 256         // so both can be done is the same loop
 257         forAll(oldFace, i)
 258         {
 259             // Add the point
 260             Map<label>::const_iterator mpIter =
 261                 pmm.find(oldFace[i]);
 262
 263             if (mpIter == pmm.end())
 264             {
 265                 // Point not mapped
 266                 newFace.append(oldFace[i]);
 267
 268                 // Add the point into patch support
 269                 pointMap().insert
 270                 (
 271                     oldFace[i],
 272                     masterLocalPoints[oldLocalFace[i]]
 273                 );
 274             }
 275             else
 276             {
 277                 // Point mapped
 278                 newFace.append(mpIter());
 279
 280                 // Add the point into support
 281                 pointMap().insert(mpIter(), masterLocalPoints[oldLocalFace[i]]);
 282             }
 283
 284             // Grab the edge points
 285
 286             const labelList& masterPointsOnEdge =
 287                 pointsIntoMasterEdges[curEdges[i]];
 288
 289             // If there are no points on the edge, skip everything
 290             // If there is only one point, no need for sorting
 291             if (masterPointsOnEdge.size())
 292             {
 293                 // Sort edge points in order
 294                 scalarField edgePointWeights(masterPointsOnEdge.size());
 295                 const point& startPoint = masterLocalPoints[oldLocalFace[i]];
 296
 297                 vector e =
 298                     masterLocalPoints[oldLocalFace.nextLabel(i)]
 299                   - startPoint;
 300
 301                 scalar magSqrE = magSqr(e);
 302
 303                 if (magSqrE > SMALL)
 304                 {
 305                     e /= magSqrE;
 306                 }
 307                 else
 308                 {
 309                     FatalErrorIn
 310                     (
 311                         "void enrichedPatch::calcEnrichedFaces\n"
 312                         "(\n"
 313                         "    const labelListList& pointsIntoMasterEdges,\n"
 314                         "    const labelListList& pointsIntoSlaveEdges,\n"
 315                         "    const pointField& projectedSlavePoints\n"
 316                         ")"
 317                     )   << "Zero length edge in master patch for face " << i
 318                         << ".  This is not allowed."
 319                         << abort(FatalError);
 320                 }
 321
 322                 pointField masterPosOnEdge(masterPointsOnEdge.size());
 323
 324                 forAll(masterPointsOnEdge, edgePointI)
 325                 {
 326                     masterPosOnEdge[edgePointI] =
 327                         pointMap().find(masterPointsOnEdge[edgePointI])();
 328
 329                     edgePointWeights[edgePointI] =
 330                         (e & (masterPosOnEdge[edgePointI] - startPoint));
 331                 }
 332
 333                 if (debug)
 334                 {
 335                     // Check weights: all new points should be on the edge
 336                     if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
 337                     {
 338                         FatalErrorIn
 339                         (
 340                             "void enrichedPatch::calcEnrichedFaces\n"
 341                             "(\n"
 342                             "    const labelListList& pointsIntoMasterEdges,\n"
 343                             "    const labelListList& pointsIntoSlaveEdges,\n"
 344                             "    const pointField& projectedSlavePoints\n"
 345                             ")"
 346                         )   << "Invalid point edge weights.  Some of points are"
 347                             << " not on the edge for edge " << curEdges[i]
 348                             << " of face " << faceI << " in master patch." << nl
 349                             << "Min weight: " << min(edgePointWeights)
 350                             << " Max weight: " << max(edgePointWeights)
 351                             << abort(FatalError);
 352                     }
 353                 }
 354
 355                 // Go through the points and collect them based on
 356                 // weights from lower to higher.  This gives the
 357                 // correct order of points along the edge.
 358                 forAll(edgePointWeights, passI)
 359                 {
 360                     // Max weight can only be one, so the sorting is
 361                     // done by elimination.
 362                     label nextPoint = -1;
 363                     scalar dist = 2;
 364
 365                     forAll(edgePointWeights, wI)
 366                     {
 367                         if (edgePointWeights[wI] < dist)
 368                         {
 369                             dist = edgePointWeights[wI];
 370                             nextPoint = wI;
 371                         }
 372                     }
 373
 374                     // Insert the next point and reset its weight to exclude it
 375                     // from future picks
 376                     newFace.append(masterPointsOnEdge[nextPoint]);
 377                     edgePointWeights[nextPoint] = GREAT;
 378
 379                     // Add the point into patch support
 380                     pointMap().insert
 381                     (
 382                         masterPointsOnEdge[nextPoint],
 383                         masterPosOnEdge[nextPoint]
 384                     );
 385                 }
 386             }
 387         }
 388         // Info<< "New master face: " << newFace << endl;
 389
 390         // Add the new face to the list
 391         enrichedFaces[nEnrichedFaces].transfer(newFace);
 392         nEnrichedFaces++;
 393     }
 394
 395     // Check the support for the enriched patch
 396     if (debug)
 397     {
 398         if (!checkSupport())
 399         {
 400             Info<< "Enriched patch support OK. Slave faces: "
 401                 << slavePatch_.size() << " Master faces: "
 402                 << masterPatch_.size() << endl;
 403         }
 404         else
 405         {
 406             FatalErrorIn
 407             (
 408                 "void enrichedPatch::calcEnrichedFaces\n"
 409                 "(\n"
 410                 "    const labelListList& pointsIntoMasterEdges,\n"
 411                 "    const labelListList& pointsIntoSlaveEdges,\n"
 412                 "    const pointField& projectedSlavePoints\n"
 413                 ")"
 414             )   << "Error in enriched patch support"
 415                 << abort(FatalError);
 416         }
 417     }
 418 }
 419
 420
 421 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 422
 423 const Foam::faceList& Foam::enrichedPatch::enrichedFaces() const
 424 {
 425     if (!enrichedFacesPtr_)
 426     {
 427         FatalErrorIn("const faceList& enrichedPatch::enrichedFaces() const")
 428             << "Enriched faces not available yet.  Please use "
 429             << "void enrichedPatch::calcEnrichedFaces\n"
 430             << "(\n"
 431             << "    const labelListList& pointsIntoMasterEdges,\n"
 432             << "    const labelListList& pointsIntoSlaveEdges,\n"
 433             << "    const pointField& projectedSlavePoints\n"
 434             << ")"
 435             << " before trying to access faces."
 436             << abort(FatalError);
 437     }
 438
 439     return *enrichedFacesPtr_;
 440 }
 441
 442
 443 // ************************************************************************* //