src/dynamicMesh/polyMeshModifiers/perfectInterface/perfectInterface.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     Best thing is probably to look at attachDetach which does almost exactly
  27     the same but for the geometric matching of points and face centres.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "perfectInterface.H"
  32 #include "polyTopoChanger.H"
  33 #include "polyMesh.H"
  34 #include "polyTopoChange.H"
  35 #include "addToRunTimeSelectionTable.H"
  36 #include "mapPolyMesh.H"
  37 #include "matchPoints.H"
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43     defineTypeNameAndDebug(perfectInterface, 0);
  44     addToRunTimeSelectionTable
  45     (
  46         polyMeshModifier,
  47         perfectInterface,
  48         dictionary
  49     );
  50 }
  51
  52
  53 // Tolerance used as fraction of minimum edge length.
  54 const Foam::scalar Foam::perfectInterface::tol_ = 1e-3;
  55
  56
  57 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  58
  59 Foam::pointField Foam::perfectInterface::calcFaceCentres
  60 (
  61     const primitivePatch& pp
  62 )
  63 {
  64     const pointField& points = pp.points();
  65
  66     pointField ctrs(pp.size());
  67
  68     forAll(ctrs, patchFaceI)
  69     {
  70         ctrs[patchFaceI] = pp[patchFaceI].centre(points);
  71     }
  72
  73     return ctrs;
  74 }
  75
  76
  77 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  78
  79 // Construct from components
  80 Foam::perfectInterface::perfectInterface
  81 (
  82     const word& name,
  83     const label index,
  84     const polyTopoChanger& mme,
  85     const word& faceZoneName,
  86     const word& masterPatchName,
  87     const word& slavePatchName
  88 )
  89 :
  90     polyMeshModifier(name, index, mme, true),
  91     faceZoneID_(faceZoneName, mme.mesh().faceZones()),
  92     masterPatchID_(masterPatchName, mme.mesh().boundaryMesh()),
  93     slavePatchID_(slavePatchName, mme.mesh().boundaryMesh())
  94 {}
  95
  96
  97 // Construct from dictionary
  98 Foam::perfectInterface::perfectInterface
  99 (
 100     const word& name,
 101     const dictionary& dict,
 102     const label index,
 103     const polyTopoChanger& mme
 104 )
 105 :
 106     polyMeshModifier(name, index, mme, readBool(dict.lookup("active"))),
 107     faceZoneID_
 108     (
 109         dict.lookup("faceZoneName"),
 110         mme.mesh().faceZones()
 111     ),
 112     masterPatchID_
 113     (
 114         dict.lookup("masterPatchName"),
 115         mme.mesh().boundaryMesh()
 116     ),
 117     slavePatchID_
 118     (
 119         dict.lookup("slavePatchName"),
 120         mme.mesh().boundaryMesh()
 121     )
 122 {}
 123
 124
 125 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 126
 127 Foam::perfectInterface::~perfectInterface()
 128 {}
 129
 130
 131 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 132
 133 bool Foam::perfectInterface::changeTopology() const
 134 {
 135     // If modifier is inactive, skip change
 136     if (!active())
 137     {
 138         if (debug)
 139         {
 140             Pout<< "bool perfectInterface::changeTopology() const "
 141                 << "for object " << name() << " : "
 142                 << "Inactive" << endl;
 143         }
 144
 145         return false;
 146     }
 147     else
 148     {
 149         // I want topo change every time step.
 150         return true;
 151     }
 152 }
 153
 154
 155 void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
 156 {
 157     if (debug)
 158     {
 159         Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
 160             << "for object " << name() << " : "
 161             << "masterPatchID_:" << masterPatchID_
 162             << " slavePatchID_:" << slavePatchID_
 163             << " faceZoneID_:" << faceZoneID_ << endl;
 164     }
 165
 166     if
 167     (
 168         masterPatchID_.active()
 169      && slavePatchID_.active()
 170      && faceZoneID_.active()
 171     )
 172     {
 173         const polyMesh& mesh = topoChanger().mesh();
 174
 175         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 176
 177         const polyPatch& pp0 = patches[masterPatchID_.index()];
 178         const polyPatch& pp1 = patches[slavePatchID_.index()];
 179
 180         // Some aliases
 181         const edgeList& edges0 = pp0.edges();
 182         const pointField& pts0 = pp0.localPoints();
 183         const pointField& pts1 = pp1.localPoints();
 184         const labelList& meshPts0 = pp0.meshPoints();
 185         const labelList& meshPts1 = pp1.meshPoints();
 186
 187
 188         // Get local dimension as fraction of minimum edge length
 189
 190         scalar minLen = GREAT;
 191
 192         forAll(edges0, edgeI)
 193         {
 194             minLen = min(minLen, edges0[edgeI].mag(pts0));
 195         }
 196         scalar typDim = tol_*minLen;
 197
 198         if (debug)
 199         {
 200             Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
 201                 << " pts0:" << pts0.size() << " pts1:" << pts1.size()
 202                 << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
 203         }
 204
 205
 206         // Determine pointMapping in mesh point labels. Uses geometric
 207         // comparison to find correspondence between patch points.
 208
 209         labelList renumberPoints(mesh.allPoints().size());
 210         forAll(renumberPoints, i)
 211         {
 212             renumberPoints[i] = i;
 213         }
 214         {
 215             labelList from1To0Points(pts1.size());
 216
 217             bool matchOk = matchPoints
 218             (
 219                 pts1,
 220                 pts0,
 221                 scalarField(pts1.size(), typDim),   // tolerance
 222                 true,                               // verbose
 223                 from1To0Points
 224             );
 225
 226             if (!matchOk)
 227             {
 228                 FatalErrorIn
 229                 (
 230                     "perfectInterface::setRefinement(polyTopoChange& ref) const"
 231                 )   << "Points on patches " << pp0.name() << " and "
 232                     << pp1.name() << " do not match to within tolerance "
 233                     << typDim << exit(FatalError);
 234             }
 235
 236             forAll(pts1, i)
 237             {
 238                 renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
 239             }
 240         }
 241
 242
 243
 244         // Calculate correspondence between patch faces
 245
 246         labelList from0To1Faces(pp1.size());
 247
 248         bool matchOk = matchPoints
 249         (
 250             calcFaceCentres(pp0),
 251             calcFaceCentres(pp1),
 252             scalarField(pp0.size(), typDim),    // tolerance
 253             true,                               // verbose
 254             from0To1Faces
 255         );
 256
 257         if (!matchOk)
 258         {
 259             FatalErrorIn
 260             (
 261                 "perfectInterface::setRefinement(polyTopoChange& ref) const"
 262             )   << "Face centres of patches " << pp0.name() << " and "
 263                 << pp1.name() << " do not match to within tolerance " << typDim
 264                 << exit(FatalError);
 265         }
 266
 267
 268
 269         // Now
 270         // - renumber faces using pts1 (except patch1 faces)
 271         // - remove patch1 faces. Remember cell label on owner side.
 272         // - modify patch0 faces to be internal.
 273
 274         // 1. Get faces to be renumbered
 275         labelHashSet affectedFaces(2*pp1.size());
 276         forAll(meshPts1, i)
 277         {
 278             label meshPointI = meshPts1[i];
 279
 280             if (meshPointI != renumberPoints[meshPointI])
 281             {
 282                 const labelList& pFaces = mesh.pointFaces()[meshPointI];
 283
 284                 forAll(pFaces, pFaceI)
 285                 {
 286                     affectedFaces.insert(pFaces[pFaceI]);
 287                 }
 288             }
 289         }
 290         forAll(pp1, i)
 291         {
 292             affectedFaces.erase(pp1.start() + i);
 293         }
 294         // Remove patch0 from renumbered faces. Should not be nessecary since
 295         // patch0 and 1 should not share any point (if created by mergeMeshing)
 296         // so affectedFaces should not contain any patch0 faces but you can
 297         // never be sure what the user is doing.
 298         forAll(pp0, i)
 299         {
 300             if (affectedFaces.erase(pp0.start() + i))
 301             {
 302                 WarningIn
 303                 (
 304                     "perfectInterface::setRefinement(polyTopoChange&) const"
 305                 )   << "Found face " << pp0.start() + i << " vertices "
 306                     << mesh.faces()[pp0.start() + i] << " whose points are"
 307                     << " used both by master patch " << pp0.name()
 308                     << " and slave patch " << pp1.name()
 309                     << endl;
 310             }
 311         }
 312
 313
 314         // 2. Renumber (non patch0/1) faces.
 315         for
 316         (
 317             labelHashSet::const_iterator iter = affectedFaces.begin();
 318             iter != affectedFaces.end();
 319             ++iter
 320         )
 321         {
 322             label faceI = iter.key();
 323
 324             const face& f = mesh.faces()[faceI];
 325
 326             face newFace(f.size());
 327
 328             forAll(newFace, fp)
 329             {
 330                 newFace[fp] = renumberPoints[f[fp]];
 331             }
 332
 333             label nbr = -1;
 334
 335             label patchI = -1;
 336
 337             if (mesh.isInternalFace(faceI))
 338             {
 339                 nbr = mesh.faceNeighbour()[faceI];
 340             }
 341             else
 342             {
 343                 patchI = patches.whichPatch(faceI);
 344             }
 345
 346             label zoneID = mesh.faceZones().whichZone(faceI);
 347
 348             bool zoneFlip = false;
 349
 350             if (zoneID >= 0)
 351             {
 352                 const faceZone& fZone = mesh.faceZones()[zoneID];
 353
 354                 zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 355             }
 356
 357             ref.setAction
 358             (
 359                 polyModifyFace
 360                 (
 361                     newFace,                    // modified face
 362                     faceI,                      // label of face being modified
 363                     mesh.faceOwner()[faceI],    // owner
 364                     nbr,                        // neighbour
 365                     false,                      // face flip
 366                     patchI,                     // patch for face
 367                     false,                      // remove from zone
 368                     zoneID,                     // zone for face
 369                     zoneFlip                    // face flip in zone
 370                 )
 371             );
 372         }
 373
 374
 375         // 3. Remove patch1 points
 376         forAll(meshPts1, i)
 377         {
 378             label meshPointI = meshPts1[i];
 379
 380             if (meshPointI != renumberPoints[meshPointI])
 381             {
 382                 ref.setAction(polyRemovePoint(meshPointI));
 383             }
 384         }
 385
 386
 387         // 4. Remove patch1 faces
 388         forAll(pp1, i)
 389         {
 390             ref.setAction(polyRemoveFace(pp1.start() + i));
 391         }
 392
 393
 394         // 5. Modify patch0 faces for new points (not really nessecary; see
 395         // comment above about patch1 and patch0 never sharing points) and
 396         // becoming internal.
 397         const boolList& mfFlip =
 398             mesh.faceZones()[faceZoneID_.index()].flipMap();
 399
 400         forAll(pp0, i)
 401         {
 402             label faceI = pp0.start() + i;
 403
 404             const face& f = mesh.faces()[faceI];
 405
 406             face newFace(f.size());
 407
 408             forAll(newFace, fp)
 409             {
 410                 newFace[fp] = renumberPoints[f[fp]];
 411             }
 412
 413             label own = mesh.faceOwner()[faceI];
 414
 415             label pp1FaceI = pp1.start() + from0To1Faces[i];
 416
 417             label nbr = mesh.faceOwner()[pp1FaceI];
 418
 419             if (own < nbr)
 420             {
 421                 ref.setAction
 422                 (
 423                     polyModifyFace
 424                     (
 425                         newFace,                // modified face
 426                         faceI,                  // label of face being modified
 427                         own,                    // owner
 428                         nbr,                    // neighbour
 429                         false,                  // face flip
 430                         -1,                     // patch for face
 431                         false,                  // remove from zone
 432                         faceZoneID_.index(),    // zone for face
 433                         mfFlip[i]               // face flip in zone
 434                     )
 435                 );
 436             }
 437             else
 438             {
 439                 ref.setAction
 440                 (
 441                     polyModifyFace
 442                     (
 443                         newFace.reverseFace(),  // modified face
 444                         faceI,                  // label of face being modified
 445                         nbr,                    // owner
 446                         own,                    // neighbour
 447                         false,                  // face flip
 448                         -1,                     // patch for face
 449                         false,                  // remove from zone
 450                         faceZoneID_.index(),    // zone for face
 451                         !mfFlip[i]              // face flip in zone
 452                     )
 453                 );
 454             }
 455         }
 456     }
 457 }
 458
 459
 460 void Foam::perfectInterface::modifyMotionPoints(pointField& motionPoints) const
 461 {
 462     // Update only my points. Nothing to be done here as points already
 463     // shared by now.
 464 }
 465
 466
 467 void Foam::perfectInterface::updateMesh(const mapPolyMesh& morphMap)
 468 {
 469     // Mesh has changed topologically.  Update local topological data
 470     const polyMesh& mesh = topoChanger().mesh();
 471
 472     faceZoneID_.update(mesh.faceZones());
 473     masterPatchID_.update(mesh.boundaryMesh());
 474     slavePatchID_.update(mesh.boundaryMesh());
 475 }
 476
 477
 478 void Foam::perfectInterface::write(Ostream& os) const
 479 {
 480     os  << nl << type() << nl
 481         << name()<< nl
 482         << faceZoneID_.name() << nl
 483         << masterPatchID_.name() << nl
 484         << slavePatchID_.name() << endl;
 485 }
 486
 487
 488 void Foam::perfectInterface::writeDict(Ostream& os) const
 489 {
 490     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
 491
 492         << "    type " << type()
 493         << token::END_STATEMENT << nl
 494
 495         << "    active " << active()
 496         << token::END_STATEMENT << nl
 497
 498         << "    faceZoneName " << faceZoneID_.name()
 499         << token::END_STATEMENT << nl
 500
 501         << "    masterPatchName " << masterPatchID_.name()
 502         << token::END_STATEMENT << nl
 503
 504         << "    slavePatchName " << slavePatchID_.name()
 505         << token::END_STATEMENT << nl
 506
 507         << token::END_BLOCK << endl;
 508 }
 509
 510
 511 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 512
 513
 514 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 515
 516
 517 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 518
 519
 520 // ************************************************************************* //