src/dynamicMesh/perfectInterface/perfectInterface.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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 #include "indirectPrimitivePatch.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 indirectPrimitivePatch& 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
 159 (
 160     const indirectPrimitivePatch& pp0,
 161     const indirectPrimitivePatch& pp1,
 162     polyTopoChange& ref
 163 ) const
 164 {
 165     const polyMesh& mesh = topoChanger().mesh();
 166
 167     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 168
 169     // Some aliases
 170     const edgeList& edges0 = pp0.edges();
 171     const pointField& pts0 = pp0.localPoints();
 172     const pointField& pts1 = pp1.localPoints();
 173     const labelList& meshPts0 = pp0.meshPoints();
 174     const labelList& meshPts1 = pp1.meshPoints();
 175
 176
 177     // Get local dimension as fraction of minimum edge length
 178
 179     scalar minLen = GREAT;
 180
 181     forAll(edges0, edgeI)
 182     {
 183         minLen = min(minLen, edges0[edgeI].mag(pts0));
 184     }
 185     scalar typDim = tol_*minLen;
 186
 187     if (debug)
 188     {
 189         Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
 190             << " pts0:" << pts0.size() << " pts1:" << pts1.size()
 191             << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
 192     }
 193
 194
 195     // Determine pointMapping in mesh point labels. Uses geometric
 196     // comparison to find correspondence between patch points.
 197
 198     labelList renumberPoints(mesh.points().size());
 199     forAll(renumberPoints, i)
 200     {
 201         renumberPoints[i] = i;
 202     }
 203     {
 204         labelList from1To0Points(pts1.size());
 205
 206         bool matchOk = matchPoints
 207         (
 208             pts1,
 209             pts0,
 210             scalarField(pts1.size(), typDim),   // tolerance
 211             true,                               // verbose
 212             from1To0Points
 213         );
 214
 215         if (!matchOk)
 216         {
 217             FatalErrorIn
 218             (
 219                 "perfectInterface::setRefinement(polyTopoChange& ref) const"
 220             )   << "Points on patch sides do not match to within tolerance "
 221                 << typDim << exit(FatalError);
 222         }
 223
 224         forAll(pts1, i)
 225         {
 226             renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
 227         }
 228     }
 229
 230
 231
 232     // Calculate correspondence between patch faces
 233
 234     labelList from0To1Faces(pp1.size());
 235
 236     bool matchOk = matchPoints
 237     (
 238         calcFaceCentres(pp0),
 239         calcFaceCentres(pp1),
 240         scalarField(pp0.size(), typDim),    // tolerance
 241         true,                               // verbose
 242         from0To1Faces
 243     );
 244
 245     if (!matchOk)
 246     {
 247         FatalErrorIn
 248         (
 249             "perfectInterface::setRefinement(polyTopoChange& ref) const"
 250         )   << "Face centres of patch sides do not match to within tolerance "
 251             << typDim << exit(FatalError);
 252     }
 253
 254
 255
 256     // Now
 257     // - renumber faces using pts1 (except patch1 faces)
 258     // - remove patch1 faces. Remember cell label on owner side.
 259     // - modify patch0 faces to be internal.
 260
 261     // 1. Get faces to be renumbered
 262     labelHashSet affectedFaces(2*pp1.size());
 263     forAll(meshPts1, i)
 264     {
 265         label meshPointI = meshPts1[i];
 266
 267         if (meshPointI != renumberPoints[meshPointI])
 268         {
 269             const labelList& pFaces = mesh.pointFaces()[meshPointI];
 270
 271             forAll(pFaces, pFaceI)
 272             {
 273                 affectedFaces.insert(pFaces[pFaceI]);
 274             }
 275         }
 276     }
 277     forAll(pp1, i)
 278     {
 279         affectedFaces.erase(pp1.addressing()[i]);
 280     }
 281     // Remove patch0 from renumbered faces. Should not be nessecary since
 282     // patch0 and 1 should not share any point (if created by mergeMeshing)
 283     // so affectedFaces should not contain any patch0 faces but you can
 284     // never be sure what the user is doing.
 285     forAll(pp0, i)
 286     {
 287         label faceI = pp0.addressing()[i];
 288
 289         if (affectedFaces.erase(faceI))
 290         {
 291             WarningIn
 292             (
 293                 "perfectInterface::setRefinement(polyTopoChange&) const"
 294             )   << "Found face " << faceI << " vertices "
 295                 << mesh.faces()[faceI] << " whose points are"
 296                 << " used both by master patch and slave patch" << endl;
 297         }
 298     }
 299
 300
 301     // 2. Renumber (non patch0/1) faces.
 302     forAllConstIter(labelHashSet, affectedFaces, iter)
 303     {
 304         const label faceI = iter.key();
 305         const face& f = mesh.faces()[faceI];
 306
 307         face newFace(f.size());
 308
 309         forAll(newFace, fp)
 310         {
 311             newFace[fp] = renumberPoints[f[fp]];
 312         }
 313
 314         label nbr = -1;
 315
 316         label patchI = -1;
 317
 318         if (mesh.isInternalFace(faceI))
 319         {
 320             nbr = mesh.faceNeighbour()[faceI];
 321         }
 322         else
 323         {
 324             patchI = patches.whichPatch(faceI);
 325         }
 326
 327         label zoneID = mesh.faceZones().whichZone(faceI);
 328
 329         bool zoneFlip = false;
 330
 331         if (zoneID >= 0)
 332         {
 333             const faceZone& fZone = mesh.faceZones()[zoneID];
 334
 335             zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 336         }
 337
 338         ref.setAction
 339         (
 340             polyModifyFace
 341             (
 342                 newFace,                    // modified face
 343                 faceI,                      // label of face being modified
 344                 mesh.faceOwner()[faceI],    // owner
 345                 nbr,                        // neighbour
 346                 false,                      // face flip
 347                 patchI,                     // patch for face
 348                 false,                      // remove from zone
 349                 zoneID,                     // zone for face
 350                 zoneFlip                    // face flip in zone
 351             )
 352         );
 353     }
 354
 355
 356     // 3. Remove patch1 points
 357     forAll(meshPts1, i)
 358     {
 359         label meshPointI = meshPts1[i];
 360
 361         if (meshPointI != renumberPoints[meshPointI])
 362         {
 363             ref.setAction(polyRemovePoint(meshPointI));
 364         }
 365     }
 366
 367
 368     // 4. Remove patch1 faces
 369     forAll(pp1, i)
 370     {
 371         label faceI = pp1.addressing()[i];
 372         ref.setAction(polyRemoveFace(faceI));
 373     }
 374
 375
 376     // 5. Modify patch0 faces for new points (not really nessecary; see
 377     // comment above about patch1 and patch0 never sharing points) and
 378     // becoming internal.
 379     const boolList& mfFlip =
 380         mesh.faceZones()[faceZoneID_.index()].flipMap();
 381
 382     forAll(pp0, i)
 383     {
 384         label faceI = pp0.addressing()[i];
 385
 386         const face& f = mesh.faces()[faceI];
 387
 388         face newFace(f.size());
 389
 390         forAll(newFace, fp)
 391         {
 392             newFace[fp] = renumberPoints[f[fp]];
 393         }
 394
 395         label own = mesh.faceOwner()[faceI];
 396
 397         label pp1FaceI = pp1.addressing()[from0To1Faces[i]];
 398
 399         label nbr = mesh.faceOwner()[pp1FaceI];
 400
 401         if (own < nbr)
 402         {
 403             ref.setAction
 404             (
 405                 polyModifyFace
 406                 (
 407                     newFace,                // modified face
 408                     faceI,                  // label of face being modified
 409                     own,                    // owner
 410                     nbr,                    // neighbour
 411                     false,                  // face flip
 412                     -1,                     // patch for face
 413                     false,                  // remove from zone
 414                     faceZoneID_.index(),    // zone for face
 415                     mfFlip[i]               // face flip in zone
 416                 )
 417             );
 418         }
 419         else
 420         {
 421             ref.setAction
 422             (
 423                 polyModifyFace
 424                 (
 425                     newFace.reverseFace(),  // modified face
 426                     faceI,                  // label of face being modified
 427                     nbr,                    // owner
 428                     own,                    // neighbour
 429                     true,                   // 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     }
 438 }
 439
 440
 441 void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
 442 {
 443     if (debug)
 444     {
 445         Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
 446             << "for object " << name() << " : "
 447             << "masterPatchID_:" << masterPatchID_
 448             << " slavePatchID_:" << slavePatchID_
 449             << " faceZoneID_:" << faceZoneID_ << endl;
 450     }
 451
 452     if
 453     (
 454         masterPatchID_.active()
 455      && slavePatchID_.active()
 456      && faceZoneID_.active()
 457     )
 458     {
 459         const polyMesh& mesh = topoChanger().mesh();
 460
 461         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 462         const polyPatch& patch0 = patches[masterPatchID_.index()];
 463         const polyPatch& patch1 = patches[slavePatchID_.index()];
 464
 465
 466         labelList pp0Labels(identity(patch0.size())+patch0.start());
 467         indirectPrimitivePatch pp0
 468         (
 469             IndirectList<face>(mesh.faces(), pp0Labels),
 470             mesh.points()
 471         );
 472
 473         labelList pp1Labels(identity(patch1.size())+patch1.start());
 474         indirectPrimitivePatch pp1
 475         (
 476             IndirectList<face>(mesh.faces(), pp1Labels),
 477             mesh.points()
 478         );
 479
 480         setRefinement(pp0, pp1, ref);
 481     }
 482 }
 483
 484
 485 void Foam::perfectInterface::modifyMotionPoints(pointField& motionPoints) const
 486 {
 487     // Update only my points. Nothing to be done here as points already
 488     // shared by now.
 489 }
 490
 491
 492 void Foam::perfectInterface::updateMesh(const mapPolyMesh& morphMap)
 493 {
 494     // Mesh has changed topologically.  Update local topological data
 495     const polyMesh& mesh = topoChanger().mesh();
 496
 497     faceZoneID_.update(mesh.faceZones());
 498     masterPatchID_.update(mesh.boundaryMesh());
 499     slavePatchID_.update(mesh.boundaryMesh());
 500 }
 501
 502
 503 void Foam::perfectInterface::write(Ostream& os) const
 504 {
 505     os  << nl << type() << nl
 506         << name()<< nl
 507         << faceZoneID_.name() << nl
 508         << masterPatchID_.name() << nl
 509         << slavePatchID_.name() << endl;
 510 }
 511
 512
 513 void Foam::perfectInterface::writeDict(Ostream& os) const
 514 {
 515     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
 516
 517         << "    type " << type()
 518         << token::END_STATEMENT << nl
 519
 520         << "    active " << active()
 521         << token::END_STATEMENT << nl
 522
 523         << "    faceZoneName " << faceZoneID_.name()
 524         << token::END_STATEMENT << nl
 525
 526         << "    masterPatchName " << masterPatchID_.name()
 527         << token::END_STATEMENT << nl
 528
 529         << "    slavePatchName " << slavePatchID_.name()
 530         << token::END_STATEMENT << nl
 531
 532         << token::END_BLOCK << endl;
 533 }
 534
 535
 536 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 537
 538
 539 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 540
 541
 542 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 543
 544
 545 // ************************************************************************* //