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