src/dynamicMesh/perfectInterface/perfectInterface.C

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