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