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