FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

cyclicPolyPatch.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "cyclicPolyPatch.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/polyBoundaryMesh.H>
00029 #include <OpenFOAM/polyMesh.H>
00030 #include <OpenFOAM/demandDrivenData.H>
00031 #include <OpenFOAM/OFstream.H>
00032 #include <OpenFOAM/patchZones.H>
00033 #include <OpenFOAM/matchPoints.H>
00034 #include <OpenFOAM/EdgeMap.H>
00035 #include <OpenFOAM/Time.H>
00036 #include <OpenFOAM/transformList.H>
00037 
00038 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00039 
00040 namespace Foam
00041 {
00042     defineTypeNameAndDebug(cyclicPolyPatch, 0);
00043 
00044     addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word);
00045     addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary);
00046 }
00047 
00048 
00049 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00050 
00051 Foam::label Foam::cyclicPolyPatch::findMaxArea
00052 (
00053     const pointField& points,
00054     const faceList& faces
00055 )
00056 {
00057     label maxI = -1;
00058     scalar maxAreaSqr = -GREAT;
00059 
00060     forAll(faces, faceI)
00061     {
00062         scalar areaSqr = magSqr(faces[faceI].normal(points));
00063 
00064         if (areaSqr > maxAreaSqr)
00065         {
00066             maxAreaSqr = areaSqr;
00067             maxI = faceI;
00068         }
00069     }
00070     return maxI;
00071 }
00072 
00073 
00074 void Foam::cyclicPolyPatch::calcTransforms()
00075 {
00076     if (size())
00077     {
00078         const pointField& points = this->points();
00079 
00080         // Determine geometric quantities on the two halves
00081         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00082 
00083         primitivePatch half0
00084         (
00085             SubList<face>
00086             (
00087                 *this,
00088                 size()/2
00089             ),
00090             points
00091         );
00092 
00093         pointField half0Ctrs(calcFaceCentres(half0, half0.points()));
00094 
00095         scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs));
00096 
00097         primitivePatch half1
00098         (
00099             SubList<face>
00100             (
00101                 *this,
00102                 size()/2,
00103                 size()/2
00104             ),
00105             points
00106         );
00107         pointField half1Ctrs(calcFaceCentres(half1, half1.points()));
00108 
00109         // Dump halves
00110         if (debug)
00111         {
00112             fileName casePath(boundaryMesh().mesh().time().path());
00113 
00114             fileName nm0(casePath/name()+"_half0_faces.obj");
00115             Pout<< "cyclicPolyPatch::calcTransforms : Writing half0"
00116                 << " faces to OBJ file " << nm0 << endl;
00117             writeOBJ(nm0, half0, half0.points());
00118 
00119             fileName nm1(casePath/name()+"_half1_faces.obj");
00120             Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
00121                 << " faces to OBJ file " << nm1 << endl;
00122             writeOBJ(nm1, half1, half1.points());
00123 
00124             OFstream str(casePath/name()+"_half0_to_half1.obj");
00125             label vertI = 0;
00126             Pout<< "cyclicPolyPatch::calcTransforms :"
00127                 << " Writing coupled face centres as lines to " << str.name()
00128                 << endl;
00129             forAll(half0Ctrs, i)
00130             {
00131                 const point& p0 = half0Ctrs[i];
00132                 str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
00133                 vertI++;
00134                 const point& p1 = half1Ctrs[i];
00135                 str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
00136                 vertI++;
00137                 str << "l " << vertI-1 << ' ' << vertI << nl;
00138             }
00139         }
00140 
00141         vectorField half0Normals(half0.size());
00142         vectorField half1Normals(half1.size());
00143 
00144         scalar maxAreaDiff = -GREAT;
00145         label maxAreaFacei = -1;
00146 
00147         for (label facei = 0; facei < size()/2; facei++)
00148         {
00149             half0Normals[facei] = operator[](facei).normal(points);
00150             label nbrFacei = facei+size()/2;
00151             half1Normals[facei] = operator[](nbrFacei).normal(points);
00152 
00153             scalar magSf = mag(half0Normals[facei]);
00154             scalar nbrMagSf = mag(half1Normals[facei]);
00155             scalar avSf = (magSf + nbrMagSf)/2.0;
00156 
00157             if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
00158             {
00159                 // Undetermined normal. Use dummy normal to force separation
00160                 // check. (note use of sqrt(VSMALL) since that is how mag
00161                 // scales)
00162                 half0Normals[facei] = point(1, 0, 0);
00163                 half1Normals[facei] = half0Normals[facei];
00164             }
00165             else
00166             {
00167                 scalar areaDiff = mag(magSf - nbrMagSf)/avSf;
00168 
00169                 if (areaDiff > maxAreaDiff)
00170                 {
00171                     maxAreaDiff = areaDiff;
00172                     maxAreaFacei = facei;
00173                 }
00174 
00175                 if (areaDiff > coupledPolyPatch::matchTol)
00176                 {
00177                     FatalErrorIn
00178                     (
00179                         "cyclicPolyPatch::calcTransforms()"
00180                     )   << "face " << facei << " area does not match neighbour "
00181                         << nbrFacei << " by "
00182                         << 100*areaDiff
00183                         << "% -- possible face ordering problem." << endl
00184                         << "patch:" << name()
00185                         << " my area:" << magSf
00186                         << " neighbour area:" << nbrMagSf
00187                         << " matching tolerance:"
00188                         << coupledPolyPatch::matchTol
00189                          << endl
00190                         << "Mesh face:" << start()+facei
00191                         << " vertices:"
00192                         << UIndirectList<point>(points, operator[](facei))()
00193                         << endl
00194                         << "Neighbour face:" << start()+nbrFacei
00195                         << " vertices:"
00196                         << UIndirectList<point>(points, operator[](nbrFacei))()
00197                         << endl
00198                         << "Rerun with cyclic debug flag set"
00199                         << " for more information." << exit(FatalError);
00200                 }
00201                 else
00202                 {
00203                     half0Normals[facei] /= magSf;
00204                     half1Normals[facei] /= nbrMagSf;
00205                 }
00206             }
00207         }
00208 
00209 
00210         // Print area match
00211         if (debug)
00212         {
00213             label nbrFacei = maxAreaFacei+size()/2;
00214             Pout<< "cyclicPolyPatch::calcTransforms :"
00215                 << " Max area error:" << 100*maxAreaDiff << "% at face:"
00216                 << maxAreaFacei << " at:" << half0Ctrs[maxAreaFacei]
00217                 << " coupled face:" << nbrFacei
00218                 << " at:" << half1Ctrs[maxAreaFacei]
00219                 << endl;
00220         }
00221 
00222 
00223 
00224         // See if transformation is prescribed
00225         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00226 
00227         switch (transform_)
00228         {
00229             case ROTATIONAL:
00230             {
00231                 // Specified single rotation tensor.
00232 
00233                 // Get best fitting face and its opposite number
00234                 label face0 = getConsistentRotationFace(half0Ctrs);
00235                 label face1 = face0;
00236 
00237                 vector n0 =
00238                     (
00239                         (half0Ctrs[face0] - rotationCentre_)
00240                       ^ rotationAxis_
00241                     );
00242                 vector n1 =
00243                     (
00244                         (half1Ctrs[face1] - rotationCentre_)
00245                       ^ -rotationAxis_
00246                     );
00247                 n0 /= mag(n0) + VSMALL;
00248                 n1 /= mag(n1) + VSMALL;
00249 
00250                 if (debug)
00251                 {
00252                     Pout<< "cyclicPolyPatch::calcTransforms :"
00253                         << " Specified rotation :"
00254                         << " n0:" << n0 << " from face " << half0Ctrs[face0]
00255                         << " and n1:" << n1 << " from face " << half1Ctrs[face1]
00256                         << endl;
00257                 }
00258 
00259                 // Calculate transformation tensors from face0,1 only.
00260                 // Note: can use tight tolerance now.
00261                 calcTransformTensors
00262                 (
00263                     pointField(1, half0Ctrs[face0]),
00264                     pointField(1, half1Ctrs[face1]),
00265                     vectorField(1, n0),
00266                     vectorField(1, n1),
00267                     scalarField(1, half0Tols[face0]),
00268                     1E-4,
00269                     ROTATIONAL
00270                 );
00271 
00272                 break;
00273             }
00274 
00275             case TRANSLATIONAL:
00276             {
00277                 // Calculate transformation tensors from all faces just to
00278                 // compare against user provided tolerance.
00279                 calcTransformTensors
00280                 (
00281                     half0Ctrs,
00282                     half1Ctrs,
00283                     half0Normals,
00284                     half1Normals,
00285                     half0Tols,
00286                     matchTol,
00287                     transform_
00288                 );
00289 
00290                 if (debug)
00291                 {
00292                     Pout<< "cyclicPolyPatch::calcTransforms :"
00293                         << " Specified separation vector : "
00294                         << separationVector_
00295                         << " . Calculated average separation : "
00296                         << average(coupledPolyPatch::separation())
00297                         << endl;
00298                 }
00299 
00300                 // Override computed transform with specified.
00301                 const scalar avgTol = average(half0Tols);
00302                 if
00303                 (
00304                     coupledPolyPatch::separation().size() != 1
00305                  || mag(separation(0) - separationVector_) > avgTol
00306                 )
00307                 {
00308                     WarningIn
00309                     (
00310                         "cyclicPolyPatch::calcTransforms()"
00311                     )   << "Specified separationVector " << separationVector_
00312                         << " differs from computed separation vector "
00313                         << coupledPolyPatch::separation() << endl
00314                         << "This probably means your geometry is not consistent"
00315                         << " with the specified separation and might lead"
00316                         << " to problems." << endl
00317                         << "Continuing with specified separation vector "
00318                         << separationVector_ << endl
00319                         << "patch:" << name()
00320                         << endl;
00321                 }
00322 
00323                 // Set transformation tensor.
00324                 const_cast<tensorField&>(coupledPolyPatch::forwardT()).clear();
00325                 const_cast<tensorField&>(coupledPolyPatch::reverseT()).clear();
00326                 const_cast<vectorField&>(coupledPolyPatch::separation()) =
00327                 vectorField
00328                 (
00329                     1,
00330                     separationVector_
00331                 );
00332 
00333                 break;
00334             }
00335 
00336             default:
00337             {
00338                 // Calculate transformation tensors from all faces.
00339                 calcTransformTensors
00340                 (
00341                     half0Ctrs,
00342                     half1Ctrs,
00343                     half0Normals,
00344                     half1Normals,
00345                     half0Tols,
00346                     matchTol,
00347                     transform_
00348                 );
00349 
00350                 break;
00351             }
00352         }
00353     }
00354 }
00355 
00356 
00357 // Get geometric zones of patch by looking at normals.
00358 // Method 1: any edge with sharpish angle is edge between two halves.
00359 //           (this will handle e.g. wedge geometries).
00360 //           Also two fully disconnected regions will be handled this way.
00361 // Method 2: sort faces into two halves based on face normal.
00362 bool Foam::cyclicPolyPatch::getGeometricHalves
00363 (
00364     const primitivePatch& pp,
00365     labelList& half0ToPatch,
00366     labelList& half1ToPatch
00367 ) const
00368 {
00369     // Calculate normals
00370     const vectorField& faceNormals = pp.faceNormals();
00371 
00372     // Find edges with sharp angles.
00373     boolList regionEdge(pp.nEdges(), false);
00374 
00375     const labelListList& edgeFaces = pp.edgeFaces();
00376 
00377     label nRegionEdges = 0;
00378 
00379     forAll(edgeFaces, edgeI)
00380     {
00381         const labelList& eFaces = edgeFaces[edgeI];
00382 
00383         // Check manifold edges for sharp angle.
00384         // (Non-manifold already handled by patchZones)
00385         if (eFaces.size() == 2)
00386         {
00387             if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
00388             {
00389                 regionEdge[edgeI] = true;
00390 
00391                 nRegionEdges++;
00392             }
00393         }
00394     }
00395 
00396 
00397     // For every face determine zone it is connected to (without crossing
00398     // any regionEdge)
00399     patchZones ppZones(pp, regionEdge);
00400 
00401     if (debug)
00402     {
00403         Pout<< "cyclicPolyPatch::getGeometricHalves : "
00404             << "Found " << nRegionEdges << " edges on patch " << name()
00405             << " where the cos of the angle between two connected faces"
00406             << " was less than " << featureCos_ << nl
00407             << "Patch divided by these and by single sides edges into "
00408             << ppZones.nZones() << " parts." << endl;
00409 
00410 
00411         // Dumping zones to obj files.
00412 
00413         labelList nZoneFaces(ppZones.nZones());
00414 
00415         for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
00416         {
00417             OFstream stream
00418             (
00419                 boundaryMesh().mesh().time().path()
00420                /name()+"_zone_"+Foam::name(zoneI)+".obj"
00421             );
00422             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone "
00423                 << zoneI << " face centres to OBJ file " << stream.name()
00424                 << endl;
00425 
00426             labelList zoneFaces(findIndices(ppZones, zoneI));
00427 
00428             forAll(zoneFaces, i)
00429             {
00430                 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
00431             }
00432 
00433             nZoneFaces[zoneI] = zoneFaces.size();
00434         }
00435     }
00436 
00437 
00438     if (ppZones.nZones() == 2)
00439     {
00440         half0ToPatch = findIndices(ppZones, 0);
00441         half1ToPatch = findIndices(ppZones, 1);
00442     }
00443     else
00444     {
00445         if (debug)
00446         {
00447             Pout<< "cyclicPolyPatch::getGeometricHalves :"
00448                 << " falling back to face-normal comparison" << endl;
00449         }
00450         label n0Faces = 0;
00451         half0ToPatch.setSize(pp.size());
00452 
00453         label n1Faces = 0;
00454         half1ToPatch.setSize(pp.size());
00455 
00456         // Compare to face 0 normal.
00457         forAll(faceNormals, faceI)
00458         {
00459             if ((faceNormals[faceI] & faceNormals[0]) > 0)
00460             {
00461                 half0ToPatch[n0Faces++] = faceI;
00462             }
00463             else
00464             {
00465                 half1ToPatch[n1Faces++] = faceI;
00466             }
00467         }
00468         half0ToPatch.setSize(n0Faces);
00469         half1ToPatch.setSize(n1Faces);
00470 
00471         if (debug)
00472         {
00473             Pout<< "cyclicPolyPatch::getGeometricHalves :"
00474                 << " Number of faces per zone:("
00475                 << n0Faces << ' ' << n1Faces << ')' << endl;
00476         }
00477     }
00478 
00479     if (half0ToPatch.size() != half1ToPatch.size())
00480     {
00481         fileName casePath(boundaryMesh().mesh().time().path());
00482 
00483         // Dump halves
00484         {
00485             fileName nm0(casePath/name()+"_half0_faces.obj");
00486             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
00487                 << " faces to OBJ file " << nm0 << endl;
00488             writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
00489 
00490             fileName nm1(casePath/name()+"_half1_faces.obj");
00491             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
00492                 << " faces to OBJ file " << nm1 << endl;
00493             writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
00494         }
00495 
00496         // Dump face centres
00497         {
00498             OFstream str0(casePath/name()+"_half0.obj");
00499             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
00500                 << " face centres to OBJ file " << str0.name() << endl;
00501 
00502             forAll(half0ToPatch, i)
00503             {
00504                 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
00505             }
00506 
00507             OFstream str1(casePath/name()+"_half1.obj");
00508             Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
00509                 << " face centres to OBJ file " << str1.name() << endl;
00510             forAll(half1ToPatch, i)
00511             {
00512                 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
00513             }
00514         }
00515 
00516         SeriousErrorIn
00517         (
00518             "cyclicPolyPatch::getGeometricHalves"
00519             "(const primitivePatch&, labelList&, labelList&) const"
00520         )   << "Patch " << name() << " gets decomposed in two zones of"
00521             << "inequal size: " << half0ToPatch.size()
00522             << " and " << half1ToPatch.size() << endl
00523             << "This means that the patch is either not two separate regions"
00524             << " or one region where the angle between the different regions"
00525             << " is not sufficiently sharp." << endl
00526             << "Please adapt the featureCos setting." << endl
00527             << "Continuing with incorrect face ordering from now on!" << endl;
00528 
00529         return false;
00530     }
00531     else
00532     {
00533         return true;
00534     }
00535 }
00536 
00537 
00538 // Given a split of faces into left and right half calculate the centres
00539 // and anchor points. Transform the left points so they align with the
00540 // right ones.
00541 void Foam::cyclicPolyPatch::getCentresAndAnchors
00542 (
00543     const primitivePatch& pp,
00544     const faceList& half0Faces,
00545     const faceList& half1Faces,
00546 
00547     pointField& ppPoints,
00548     pointField& half0Ctrs,
00549     pointField& half1Ctrs,
00550     pointField& anchors0,
00551     scalarField& tols
00552 ) const
00553 {
00554     // Get geometric data on both halves.
00555     half0Ctrs = calcFaceCentres(half0Faces, pp.points());
00556     anchors0 = getAnchorPoints(half0Faces, pp.points());
00557     half1Ctrs = calcFaceCentres(half1Faces, pp.points());
00558 
00559     switch (transform_)
00560     {
00561         case ROTATIONAL:
00562         {
00563             label face0 = getConsistentRotationFace(half0Ctrs);
00564             label face1 = getConsistentRotationFace(half1Ctrs);
00565 
00566             vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
00567             vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
00568             n0 /= mag(n0) + VSMALL;
00569             n1 /= mag(n1) + VSMALL;
00570 
00571             if (debug)
00572             {
00573                 Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
00574                     << " Specified rotation :"
00575                     << " n0:" << n0 << " n1:" << n1 << endl;
00576             }
00577 
00578             // Rotation (around origin)
00579             const tensor reverseT(rotationTensor(n0, -n1));
00580 
00581             // Rotation
00582             forAll(half0Ctrs, faceI)
00583             {
00584                 half0Ctrs[faceI] =
00585                     Foam::transform
00586                     (
00587                         reverseT,
00588                         half0Ctrs[faceI] - rotationCentre_
00589                     )
00590                   + rotationCentre_;
00591                 anchors0[faceI] =
00592                     Foam::transform
00593                     (
00594                         reverseT,
00595                         anchors0[faceI] - rotationCentre_
00596                     )
00597                   + rotationCentre_;
00598             }
00599 
00600             ppPoints =
00601                 Foam::transform
00602                 (
00603                     reverseT,
00604                     (pp.points() - rotationCentre_)()
00605                 )
00606               + rotationCentre_;
00607 
00608             break;
00609         }
00610 
00611         case TRANSLATIONAL:
00612         {
00613             // Transform 0 points.
00614 
00615             if (debug)
00616             {
00617                 Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
00618                     << "Specified translation : " << separationVector_
00619                     << endl;
00620             }
00621 
00622             half0Ctrs += separationVector_;
00623             anchors0 += separationVector_;
00624             ppPoints = pp.points() + separationVector_;
00625 
00626             break;
00627         }
00628 
00629         default:
00630         {
00631             // Assumes that cyclic is planar. This is also the initial
00632             // condition for patches without faces.
00633 
00634             // Determine the face with max area on both halves. These
00635             // two faces are used to determine the transformation tensors
00636             label max0I = findMaxArea(pp.points(), half0Faces);
00637             vector n0 = half0Faces[max0I].normal(pp.points());
00638             n0 /= mag(n0) + VSMALL;
00639 
00640             label max1I = findMaxArea(pp.points(), half1Faces);
00641             vector n1 = half1Faces[max1I].normal(pp.points());
00642             n1 /= mag(n1) + VSMALL;
00643 
00644             if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
00645             {
00646                 if (debug)
00647                 {
00648                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
00649                         << " Detected rotation :"
00650                         << " n0:" << n0 << " n1:" << n1 << endl;
00651                 }
00652 
00653                 // Rotation (around origin)
00654                 const tensor reverseT(rotationTensor(n0, -n1));
00655 
00656                 // Rotation
00657                 forAll(half0Ctrs, faceI)
00658                 {
00659                     half0Ctrs[faceI] = Foam::transform
00660                     (
00661                         reverseT,
00662                         half0Ctrs[faceI]
00663                     );
00664                     anchors0[faceI] = Foam::transform
00665                     (
00666                         reverseT,
00667                         anchors0[faceI]
00668                     );
00669                 }
00670                 ppPoints = Foam::transform(reverseT, pp.points());
00671             }
00672             else
00673             {
00674                 // Parallel translation. Get average of all used points.
00675 
00676                 primitiveFacePatch half0(half0Faces, pp.points());
00677                 const pointField& half0Pts = half0.localPoints();
00678                 const point ctr0(sum(half0Pts)/half0Pts.size());
00679 
00680                 primitiveFacePatch half1(half1Faces, pp.points());
00681                 const pointField& half1Pts = half1.localPoints();
00682                 const point ctr1(sum(half1Pts)/half1Pts.size());
00683 
00684                 if (debug)
00685                 {
00686                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
00687                         << " Detected translation :"
00688                         << " n0:" << n0 << " n1:" << n1
00689                         << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
00690                 }
00691 
00692                 half0Ctrs += ctr1 - ctr0;
00693                 anchors0 += ctr1 - ctr0;
00694                 ppPoints = pp.points() + ctr1 - ctr0;
00695             }
00696             break;
00697         }
00698     }
00699 
00700 
00701     // Calculate typical distance per face
00702     tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
00703 }
00704 
00705 
00706 // Calculates faceMap and rotation. Returns true if all ok.
00707 bool Foam::cyclicPolyPatch::matchAnchors
00708 (
00709     const bool report,
00710     const primitivePatch& pp,
00711     const labelList& half0ToPatch,
00712     const pointField& anchors0,
00713 
00714     const labelList& half1ToPatch,
00715     const faceList& half1Faces,
00716     const labelList& from1To0,
00717 
00718     const scalarField& tols,
00719 
00720     labelList& faceMap,
00721     labelList& rotation
00722 ) const
00723 {
00724     // Set faceMap such that half0 faces get first and corresponding half1
00725     // faces last.
00726 
00727     forAll(half0ToPatch, half0FaceI)
00728     {
00729         // Label in original patch
00730         label patchFaceI = half0ToPatch[half0FaceI];
00731 
00732         faceMap[patchFaceI] = half0FaceI;
00733 
00734         // No rotation
00735         rotation[patchFaceI] = 0;
00736     }
00737 
00738     bool fullMatch = true;
00739 
00740     forAll(from1To0, half1FaceI)
00741     {
00742         label patchFaceI = half1ToPatch[half1FaceI];
00743 
00744         // This face has to match the corresponding one on half0.
00745         label half0FaceI = from1To0[half1FaceI];
00746 
00747         label newFaceI = half0FaceI + pp.size()/2;
00748 
00749         faceMap[patchFaceI] = newFaceI;
00750 
00751         // Rotate patchFaceI such that its f[0] aligns with that of
00752         // the corresponding face
00753         // (which after shuffling will be at position half0FaceI)
00754 
00755         const point& wantedAnchor = anchors0[half0FaceI];
00756 
00757         rotation[newFaceI] = getRotation
00758         (
00759             pp.points(),
00760             half1Faces[half1FaceI],
00761             wantedAnchor,
00762             tols[half1FaceI]
00763         );
00764 
00765         if (rotation[newFaceI] == -1)
00766         {
00767             fullMatch = false;
00768 
00769             if (report)
00770             {
00771                 const face& f = half1Faces[half1FaceI];
00772                 SeriousErrorIn
00773                 (
00774                     "cyclicPolyPatch::matchAnchors(..)"
00775                 )   << "Patch:" << name() << " : "
00776                     << "Cannot find point on face " << f
00777                     << " with vertices:"
00778                     << UIndirectList<point>(pp.points(), f)()
00779                     << " that matches point " << wantedAnchor
00780                     << " when matching the halves of cyclic patch " << name()
00781                     << endl
00782                     << "Continuing with incorrect face ordering from now on!"
00783                     << endl;
00784             }
00785         }
00786     }
00787     return fullMatch;
00788 }
00789 
00790 
00791 Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
00792 (
00793     const pointField& faceCentres
00794 ) const
00795 {
00796     // Determine a face furthest away from the axis
00797 
00798     const scalarField magRadSqr =
00799         magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
00800 
00801     label rotFace = findMax(magRadSqr);
00802 
00803     if (debug)
00804     {
00805         Info<< "getConsistentRotationFace(const pointField&)" << nl
00806             << "    rotFace  = " << rotFace << nl
00807             << "    point    =  " << faceCentres[rotFace] << nl
00808             << "    distance = " << Foam::sqrt(magRadSqr[rotFace])
00809             << endl;
00810     }
00811 
00812     return rotFace;
00813 }
00814 
00815 
00816 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
00817 
00818 Foam::cyclicPolyPatch::cyclicPolyPatch
00819 (
00820     const word& name,
00821     const label size,
00822     const label start,
00823     const label index,
00824     const polyBoundaryMesh& bm
00825 )
00826 :
00827     coupledPolyPatch(name, size, start, index, bm),
00828     coupledPointsPtr_(NULL),
00829     coupledEdgesPtr_(NULL),
00830     featureCos_(0.9),
00831     transform_(UNKNOWN),
00832     rotationAxis_(vector::zero),
00833     rotationCentre_(point::zero),
00834     separationVector_(vector::zero)
00835 {}
00836 
00837 
00838 Foam::cyclicPolyPatch::cyclicPolyPatch
00839 (
00840     const word& name,
00841     const dictionary& dict,
00842     const label index,
00843     const polyBoundaryMesh& bm
00844 )
00845 :
00846     coupledPolyPatch(name, dict, index, bm),
00847     coupledPointsPtr_(NULL),
00848     coupledEdgesPtr_(NULL),
00849     featureCos_(0.9),
00850     transform_(UNKNOWN),
00851     rotationAxis_(vector::zero),
00852     rotationCentre_(point::zero),
00853     separationVector_(vector::zero)
00854 {
00855     if (dict.found("neighbourPatch"))
00856     {
00857         FatalIOErrorIn
00858         (
00859             "cyclicPolyPatch::cyclicPolyPatch\n"
00860             "(\n"
00861             "    const word& name,\n"
00862             "    const dictionary& dict,\n"
00863             "    const label index,\n"
00864             "    const polyBoundaryMesh& bm\n"
00865             ")",
00866             dict
00867         )   << "Found \"neighbourPatch\" entry when reading cyclic patch "
00868             << name << endl
00869             << "Is this mesh already with split cyclics?" << endl
00870             << "If so run a newer version that supports it"
00871             << ", if not comment out the \"neighbourPatch\" entry and re-run"
00872             << exit(FatalIOError);
00873     }
00874 
00875     dict.readIfPresent("featureCos", featureCos_);
00876 
00877     if (dict.found("transform"))
00878     {
00879         transform_ = transformTypeNames.read(dict.lookup("transform"));
00880         switch (transform_)
00881         {
00882             case ROTATIONAL:
00883             {
00884                 dict.lookup("rotationAxis") >> rotationAxis_;
00885                 dict.lookup("rotationCentre") >> rotationCentre_;
00886                 break;
00887             }
00888             case TRANSLATIONAL:
00889             {
00890                 dict.lookup("separationVector") >> separationVector_;
00891                 break;
00892             }
00893             default:
00894             {
00895                 // no additional info required
00896             }
00897         }
00898     }
00899 }
00900 
00901 
00902 Foam::cyclicPolyPatch::cyclicPolyPatch
00903 (
00904     const cyclicPolyPatch& pp,
00905     const polyBoundaryMesh& bm
00906 )
00907 :
00908     coupledPolyPatch(pp, bm),
00909     coupledPointsPtr_(NULL),
00910     coupledEdgesPtr_(NULL),
00911     featureCos_(pp.featureCos_),
00912     transform_(pp.transform_),
00913     rotationAxis_(pp.rotationAxis_),
00914     rotationCentre_(pp.rotationCentre_),
00915     separationVector_(pp.separationVector_)
00916 {}
00917 
00918 
00919 Foam::cyclicPolyPatch::cyclicPolyPatch
00920 (
00921     const cyclicPolyPatch& pp,
00922     const polyBoundaryMesh& bm,
00923     const label index,
00924     const label newSize,
00925     const label newStart
00926 )
00927 :
00928     coupledPolyPatch(pp, bm, index, newSize, newStart),
00929     coupledPointsPtr_(NULL),
00930     coupledEdgesPtr_(NULL),
00931     featureCos_(pp.featureCos_),
00932     transform_(pp.transform_),
00933     rotationAxis_(pp.rotationAxis_),
00934     rotationCentre_(pp.rotationCentre_),
00935     separationVector_(pp.separationVector_)
00936 {}
00937 
00938 
00939 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00940 
00941 Foam::cyclicPolyPatch::~cyclicPolyPatch()
00942 {
00943     deleteDemandDrivenData(coupledPointsPtr_);
00944     deleteDemandDrivenData(coupledEdgesPtr_);
00945 }
00946 
00947 
00948 
00949 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00950 
00951 void Foam::cyclicPolyPatch::initGeometry()
00952 {
00953     polyPatch::initGeometry();
00954 }
00955 
00956 void Foam::cyclicPolyPatch::calcGeometry()
00957 {
00958     polyPatch::calcGeometry();
00959     calcTransforms();
00960 }
00961 
00962 void Foam::cyclicPolyPatch::initMovePoints(const pointField& p)
00963 {
00964     polyPatch::initMovePoints(p);
00965 }
00966 
00967 void Foam::cyclicPolyPatch::movePoints(const pointField& p)
00968 {
00969     polyPatch::movePoints(p);
00970     calcTransforms();
00971 }
00972 
00973 void Foam::cyclicPolyPatch::initUpdateMesh()
00974 {
00975     polyPatch::initUpdateMesh();
00976 }
00977 
00978 void Foam::cyclicPolyPatch::updateMesh()
00979 {
00980     polyPatch::updateMesh();
00981     deleteDemandDrivenData(coupledPointsPtr_);
00982     deleteDemandDrivenData(coupledEdgesPtr_);
00983 }
00984 
00985 
00986 const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const
00987 {
00988     if (!coupledPointsPtr_)
00989     {
00990         // Look at cyclic patch as two halves, A and B.
00991         // Now all we know is that relative face index in halfA is same
00992         // as coupled face in halfB and also that the 0th vertex
00993         // corresponds.
00994 
00995         // From halfA point to halfB or -1.
00996         labelList coupledPoint(nPoints(), -1);
00997 
00998         for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
00999         {
01000             const face& fA = localFaces()[patchFaceA];
01001 
01002             forAll(fA, indexA)
01003             {
01004                 label patchPointA = fA[indexA];
01005 
01006                 if (coupledPoint[patchPointA] == -1)
01007                 {
01008                     const face& fB = localFaces()[patchFaceA + size()/2];
01009 
01010                     label indexB = (fB.size() - indexA) % fB.size();
01011 
01012                     // Filter out points on wedge axis
01013                     if (patchPointA != fB[indexB])
01014                     {
01015                         coupledPoint[patchPointA] = fB[indexB];
01016                     }
01017                 }
01018             }
01019         }
01020 
01021         coupledPointsPtr_ = new edgeList(nPoints());
01022         edgeList& connected = *coupledPointsPtr_;
01023 
01024         // Extract coupled points.
01025         label connectedI = 0;
01026 
01027         forAll(coupledPoint, i)
01028         {
01029             if (coupledPoint[i] != -1)
01030             {
01031                 connected[connectedI++] = edge(i, coupledPoint[i]);
01032             }
01033         }
01034 
01035         connected.setSize(connectedI);
01036 
01037         if (debug)
01038         {
01039             OFstream str
01040             (
01041                 boundaryMesh().mesh().time().path()
01042                /"coupledPoints.obj"
01043             );
01044             label vertI = 0;
01045 
01046             Pout<< "Writing file " << str.name() << " with coordinates of "
01047                 << "coupled points" << endl;
01048 
01049             forAll(connected, i)
01050             {
01051                 const point& a = points()[meshPoints()[connected[i][0]]];
01052                 const point& b = points()[meshPoints()[connected[i][1]]];
01053 
01054                 str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
01055                 str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
01056                 vertI += 2;
01057 
01058                 str<< "l " << vertI-1 << ' ' << vertI << nl;
01059             }
01060         }
01061 
01062         // Remove any addressing calculated for the coupled edges calculation
01063         const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
01064             .clearOut();
01065     }
01066     return *coupledPointsPtr_;
01067 }
01068 
01069 
01070 const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const
01071 {
01072     if (!coupledEdgesPtr_)
01073     {
01074         // Build map from points on halfA to points on halfB.
01075         const edgeList& pointCouples = coupledPoints();
01076 
01077         Map<label> aToB(2*pointCouples.size());
01078 
01079         forAll(pointCouples, i)
01080         {
01081             const edge& e = pointCouples[i];
01082 
01083             aToB.insert(e[0], e[1]);
01084         }
01085 
01086         // Map from edge on half A to points (in halfB indices)
01087         EdgeMap<label> edgeMap(nEdges());
01088 
01089         for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
01090         {
01091             const labelList& fEdges = faceEdges()[patchFaceA];
01092 
01093             forAll(fEdges, i)
01094             {
01095                 label edgeI = fEdges[i];
01096 
01097                 const edge& e = edges()[edgeI];
01098 
01099                 // Convert edge end points to corresponding points on halfB
01100                 // side.
01101                 Map<label>::const_iterator fnd0 = aToB.find(e[0]);
01102                 if (fnd0 != aToB.end())
01103                 {
01104                     Map<label>::const_iterator fnd1 = aToB.find(e[1]);
01105                     if (fnd1 != aToB.end())
01106                     {
01107                         edgeMap.insert(edge(fnd0(), fnd1()), edgeI);
01108                     }
01109                 }
01110             }
01111         }
01112 
01113         coupledEdgesPtr_ = new edgeList(nEdges()/2);
01114         edgeList& coupledEdges = *coupledEdgesPtr_;
01115         label coupleI = 0;
01116 
01117         for (label patchFaceB = size()/2; patchFaceB < size(); patchFaceB++)
01118         {
01119             const labelList& fEdges = faceEdges()[patchFaceB];
01120 
01121             forAll(fEdges, i)
01122             {
01123                 label edgeI = fEdges[i];
01124 
01125                 const edge& e = edges()[edgeI];
01126 
01127                 // Look up halfA edge from HashTable.
01128                 EdgeMap<label>::iterator iter = edgeMap.find(e);
01129 
01130                 if (iter != edgeMap.end())
01131                 {
01132                     label halfAEdgeI = iter();
01133 
01134                     // Store correspondence. Filter out edges on wedge axis.
01135                     if (halfAEdgeI != edgeI)
01136                     {
01137                         coupledEdges[coupleI++] = edge(halfAEdgeI, edgeI);
01138                     }
01139 
01140                     // Remove so we build unique list
01141                     edgeMap.erase(iter);
01142                 }
01143             }
01144         }
01145         coupledEdges.setSize(coupleI);
01146 
01147 
01148         // Some checks
01149 
01150         forAll(coupledEdges, i)
01151         {
01152             const edge& e = coupledEdges[i];
01153 
01154             if (e[0] == e[1] || e[0] < 0 || e[1] < 0)
01155             {
01156                 FatalErrorIn("cyclicPolyPatch::coupledEdges() const")
01157                     << "Problem : at position " << i
01158                     << " illegal couple:" << e
01159                     << abort(FatalError);
01160             }
01161         }
01162 
01163         if (debug)
01164         {
01165             OFstream str
01166             (
01167                 boundaryMesh().mesh().time().path()
01168                /"coupledEdges.obj"
01169             );
01170             label vertI = 0;
01171 
01172             Pout<< "Writing file " << str.name() << " with centres of "
01173                 << "coupled edges" << endl;
01174 
01175             forAll(coupledEdges, i)
01176             {
01177                 const edge& e = coupledEdges[i];
01178 
01179                 const point& a = edges()[e[0]].centre(localPoints());
01180                 const point& b = edges()[e[1]].centre(localPoints());
01181 
01182                 str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
01183                 str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
01184                 vertI += 2;
01185 
01186                 str<< "l " << vertI-1 << ' ' << vertI << nl;
01187             }
01188         }
01189 
01190         // Remove any addressing calculated for the coupled edges calculation
01191         const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
01192             .clearOut();
01193     }
01194     return *coupledEdgesPtr_;
01195 }
01196 
01197 
01198 void Foam::cyclicPolyPatch::initOrder(const primitivePatch& pp) const
01199 {}
01200 
01201 
01202 //  Return new ordering. Ordering is -faceMap: for every face index
01203 //  the new face -rotation:for every new face the clockwise shift
01204 //  of the original face. Return false if nothing changes (faceMap
01205 //  is identity, rotation is 0)
01206 bool Foam::cyclicPolyPatch::order
01207 (
01208     const primitivePatch& pp,
01209     labelList& faceMap,
01210     labelList& rotation
01211 ) const
01212 {
01213     faceMap.setSize(pp.size());
01214     faceMap = -1;
01215 
01216     rotation.setSize(pp.size());
01217     rotation = 0;
01218 
01219     if (pp.empty())
01220     {
01221         // No faces, nothing to change.
01222         return false;
01223     }
01224 
01225     if (pp.size()&1)
01226     {
01227         FatalErrorIn("cyclicPolyPatch::order(..)")
01228             << "Size of cyclic " << name() << " should be a multiple of 2"
01229             << ". It is " << pp.size() << abort(FatalError);
01230     }
01231 
01232     label halfSize = pp.size()/2;
01233 
01234     // Supplied primitivePatch already with new points.
01235     // Cyclics are limited to one transformation tensor
01236     // currently anyway (i.e. straight plane) so should not be too big a
01237     // problem.
01238 
01239 
01240     // Indices of faces on half0
01241     labelList half0ToPatch;
01242     // Indices of faces on half1
01243     labelList half1ToPatch;
01244 
01245 
01246     // 1. Test if already correctly oriented by starting from trivial ordering.
01247     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01248 
01249     half0ToPatch = identity(halfSize);
01250     half1ToPatch = half0ToPatch + halfSize;
01251 
01252     // Get faces
01253     faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
01254     faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
01255 
01256     // Get geometric quantities
01257     pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
01258     scalarField tols;
01259     getCentresAndAnchors
01260     (
01261         pp,
01262         half0Faces,
01263         half1Faces,
01264 
01265         ppPoints,
01266         half0Ctrs,
01267         half1Ctrs,
01268         anchors0,
01269         tols
01270     );
01271 
01272     if (debug)
01273     {
01274         Pout<< "half0 transformed faceCentres (avg)   : "
01275             << gAverage(half0Ctrs) << nl
01276             << "half1 untransformed faceCentres (avg) : "
01277             << gAverage(half1Ctrs) << endl;
01278     }
01279 
01280 
01281     // Geometric match of face centre vectors
01282     labelList from1To0;
01283     bool matchedAll = matchPoints
01284     (
01285         half1Ctrs,
01286         half0Ctrs,
01287         tols,
01288         false,
01289         from1To0
01290     );
01291 
01292     if (debug)
01293     {
01294         Pout<< "cyclicPolyPatch::order : test if already ordered:"
01295             << matchedAll << endl;
01296 
01297         // Dump halves
01298         fileName nm0("match1_"+name()+"_half0_faces.obj");
01299         Pout<< "cyclicPolyPatch::order : Writing half0"
01300             << " faces to OBJ file " << nm0 << endl;
01301         writeOBJ(nm0, half0Faces, ppPoints);
01302 
01303         fileName nm1("match1_"+name()+"_half1_faces.obj");
01304         Pout<< "cyclicPolyPatch::order : Writing half1"
01305             << " faces to OBJ file " << nm1 << endl;
01306         writeOBJ(nm1, half1Faces, pp.points());
01307 
01308         if (matchedAll)
01309         {
01310             OFstream ccStr
01311             (
01312                 boundaryMesh().mesh().time().path()
01313                /"match1_"+ name() + "_faceCentres.obj"
01314             );
01315             Pout<< "cyclicPolyPatch::order : "
01316                 << "Dumping currently found cyclic match as lines between"
01317                 << " corresponding face centres to file " << ccStr.name()
01318                 << endl;
01319 
01320             label vertI = 0;
01321             forAll(half1Ctrs, i)
01322             {
01323                 // Write edge between c1 and c0
01324                 const point& c0 = half0Ctrs[from1To0[i]];
01325                 const point& c1 = half1Ctrs[i];
01326                 writeOBJ(ccStr, c0, c1, vertI);
01327             }
01328         }
01329     }
01330 
01331 
01332     // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
01333     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01334 
01335     if (!matchedAll)
01336     {
01337         label faceI = 0;
01338         for (label i = 0; i < halfSize; i++)
01339         {
01340             half0ToPatch[i] = faceI++;
01341             half1ToPatch[i] = faceI++;
01342         }
01343 
01344         // And redo all matching
01345         half0Faces = UIndirectList<face>(pp, half0ToPatch);
01346         half1Faces = UIndirectList<face>(pp, half1ToPatch);
01347 
01348         getCentresAndAnchors
01349         (
01350             pp,
01351             half0Faces,
01352             half1Faces,
01353 
01354             ppPoints,
01355             half0Ctrs,
01356             half1Ctrs,
01357             anchors0,
01358             tols
01359         );
01360 
01361         // Geometric match of face centre vectors
01362         matchedAll = matchPoints
01363         (
01364             half1Ctrs,
01365             half0Ctrs,
01366             tols,
01367             false,
01368             from1To0
01369         );
01370 
01371         if (debug)
01372         {
01373             Pout<< "cyclicPolyPatch::order : test if pairwise ordered:"
01374                 << matchedAll << endl;
01375             // Dump halves
01376             fileName nm0("match2_"+name()+"_half0_faces.obj");
01377             Pout<< "cyclicPolyPatch::order : Writing half0"
01378                 << " faces to OBJ file " << nm0 << endl;
01379             writeOBJ(nm0, half0Faces, ppPoints);
01380 
01381             fileName nm1("match2_"+name()+"_half1_faces.obj");
01382             Pout<< "cyclicPolyPatch::order : Writing half1"
01383                 << " faces to OBJ file " << nm1 << endl;
01384             writeOBJ(nm1, half1Faces, pp.points());
01385 
01386             OFstream ccStr
01387             (
01388                 boundaryMesh().mesh().time().path()
01389                /"match2_"+name()+"_faceCentres.obj"
01390             );
01391             Pout<< "cyclicPolyPatch::order : "
01392                 << "Dumping currently found cyclic match as lines between"
01393                 << " corresponding face centres to file " << ccStr.name()
01394                 << endl;
01395 
01396             // Recalculate untransformed face centres
01397             //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
01398             label vertI = 0;
01399 
01400             forAll(half1Ctrs, i)
01401             {
01402                 if (from1To0[i] != -1)
01403                 {
01404                     // Write edge between c1 and c0
01405                     //const point& c0 = rawHalf0Ctrs[from1To0[i]];
01406                     const point& c0 = half0Ctrs[from1To0[i]];
01407                     const point& c1 = half1Ctrs[i];
01408                     writeOBJ(ccStr, c0, c1, vertI);
01409                 }
01410             }
01411         }
01412     }
01413 
01414 
01415     // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
01416     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01417 
01418     if (!matchedAll)
01419     {
01420         label baffleI = 0;
01421 
01422         forAll(pp, faceI)
01423         {
01424             const face& f = pp.localFaces()[faceI];
01425             const labelList& pFaces = pp.pointFaces()[f[0]];
01426 
01427             label matchedFaceI = -1;
01428 
01429             forAll(pFaces, i)
01430             {
01431                 label otherFaceI = pFaces[i];
01432 
01433                 if (otherFaceI > faceI)
01434                 {
01435                     const face& otherF = pp.localFaces()[otherFaceI];
01436 
01437                     // Note: might pick up two similar oriented faces
01438                     //       (but that is illegal anyway)
01439                     if (f == otherF)
01440                     {
01441                         matchedFaceI = otherFaceI;
01442                         break;
01443                     }
01444                 }
01445             }
01446 
01447             if (matchedFaceI != -1)
01448             {
01449                 half0ToPatch[baffleI] = faceI;
01450                 half1ToPatch[baffleI] = matchedFaceI;
01451                 baffleI++;
01452             }
01453         }
01454 
01455         if (baffleI == halfSize)
01456         {
01457             // And redo all matching
01458             half0Faces = UIndirectList<face>(pp, half0ToPatch);
01459             half1Faces = UIndirectList<face>(pp, half1ToPatch);
01460 
01461             getCentresAndAnchors
01462             (
01463                 pp,
01464                 half0Faces,
01465                 half1Faces,
01466 
01467                 ppPoints,
01468                 half0Ctrs,
01469                 half1Ctrs,
01470                 anchors0,
01471                 tols
01472             );
01473 
01474             // Geometric match of face centre vectors
01475             matchedAll = matchPoints
01476             (
01477                 half1Ctrs,
01478                 half0Ctrs,
01479                 tols,
01480                 false,
01481                 from1To0
01482             );
01483 
01484             if (debug)
01485             {
01486                 Pout<< "cyclicPolyPatch::order : test if baffles:"
01487                     << matchedAll << endl;
01488                 // Dump halves
01489                 fileName nm0("match3_"+name()+"_half0_faces.obj");
01490                 Pout<< "cyclicPolyPatch::order : Writing half0"
01491                     << " faces to OBJ file " << nm0 << endl;
01492                 writeOBJ(nm0, half0Faces, ppPoints);
01493 
01494                 fileName nm1("match3_"+name()+"_half1_faces.obj");
01495                 Pout<< "cyclicPolyPatch::order : Writing half1"
01496                     << " faces to OBJ file " << nm1 << endl;
01497                 writeOBJ(nm1, half1Faces, pp.points());
01498 
01499                 OFstream ccStr
01500                 (
01501                     boundaryMesh().mesh().time().path()
01502                    /"match3_"+ name() + "_faceCentres.obj"
01503                 );
01504                 Pout<< "cyclicPolyPatch::order : "
01505                     << "Dumping currently found cyclic match as lines between"
01506                     << " corresponding face centres to file " << ccStr.name()
01507                     << endl;
01508 
01509                 // Recalculate untransformed face centres
01510                 //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
01511                 label vertI = 0;
01512 
01513                 forAll(half1Ctrs, i)
01514                 {
01515                     if (from1To0[i] != -1)
01516                     {
01517                         // Write edge between c1 and c0
01518                         //const point& c0 = rawHalf0Ctrs[from1To0[i]];
01519                         const point& c0 = half0Ctrs[from1To0[i]];
01520                         const point& c1 = half1Ctrs[i];
01521                         writeOBJ(ccStr, c0, c1, vertI);
01522                     }
01523                 }
01524             }
01525         }
01526     }
01527 
01528 
01529     // 4. Automatic geometric ordering
01530     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01531 
01532     if (!matchedAll)
01533     {
01534         // Split faces according to feature angle or topology
01535         bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
01536 
01537         if (!okSplit)
01538         {
01539             // Did not split into two equal parts.
01540             return false;
01541         }
01542 
01543         // And redo all matching
01544         half0Faces = UIndirectList<face>(pp, half0ToPatch);
01545         half1Faces = UIndirectList<face>(pp, half1ToPatch);
01546 
01547         getCentresAndAnchors
01548         (
01549             pp,
01550             half0Faces,
01551             half1Faces,
01552 
01553             ppPoints,
01554             half0Ctrs,
01555             half1Ctrs,
01556             anchors0,
01557             tols
01558         );
01559 
01560         // Geometric match of face centre vectors
01561         matchedAll = matchPoints
01562         (
01563             half1Ctrs,
01564             half0Ctrs,
01565             tols,
01566             false,
01567             from1To0
01568         );
01569 
01570         if (debug)
01571         {
01572             Pout<< "cyclicPolyPatch::order : automatic ordering result:"
01573                 << matchedAll << endl;
01574             // Dump halves
01575             fileName nm0("match4_"+name()+"_half0_faces.obj");
01576             Pout<< "cyclicPolyPatch::order : Writing half0"
01577                 << " faces to OBJ file " << nm0 << endl;
01578             writeOBJ(nm0, half0Faces, ppPoints);
01579 
01580             fileName nm1("match4_"+name()+"_half1_faces.obj");
01581             Pout<< "cyclicPolyPatch::order : Writing half1"
01582                 << " faces to OBJ file " << nm1 << endl;
01583             writeOBJ(nm1, half1Faces, pp.points());
01584 
01585             OFstream ccStr
01586             (
01587                 boundaryMesh().mesh().time().path()
01588                /"match4_"+ name() + "_faceCentres.obj"
01589             );
01590             Pout<< "cyclicPolyPatch::order : "
01591                 << "Dumping currently found cyclic match as lines between"
01592                 << " corresponding face centres to file " << ccStr.name()
01593                 << endl;
01594 
01595             // Recalculate untransformed face centres
01596             //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
01597             label vertI = 0;
01598 
01599             forAll(half1Ctrs, i)
01600             {
01601                 if (from1To0[i] != -1)
01602                 {
01603                     // Write edge between c1 and c0
01604                     //const point& c0 = rawHalf0Ctrs[from1To0[i]];
01605                     const point& c0 = half0Ctrs[from1To0[i]];
01606                     const point& c1 = half1Ctrs[i];
01607                     writeOBJ(ccStr, c0, c1, vertI);
01608                 }
01609             }
01610         }
01611     }
01612 
01613 
01614     if (!matchedAll || debug)
01615     {
01616         // Dump halves
01617         fileName nm0(name()+"_half0_faces.obj");
01618         Pout<< "cyclicPolyPatch::order : Writing half0"
01619             << " faces to OBJ file " << nm0 << endl;
01620         writeOBJ(nm0, half0Faces, pp.points());
01621 
01622         fileName nm1(name()+"_half1_faces.obj");
01623         Pout<< "cyclicPolyPatch::order : Writing half1"
01624             << " faces to OBJ file " << nm1 << endl;
01625         writeOBJ(nm1, half1Faces, pp.points());
01626 
01627         OFstream ccStr
01628         (
01629             boundaryMesh().mesh().time().path()
01630            /name() + "_faceCentres.obj"
01631         );
01632         Pout<< "cyclicPolyPatch::order : "
01633             << "Dumping currently found cyclic match as lines between"
01634             << " corresponding face centres to file " << ccStr.name()
01635             << endl;
01636 
01637         // Recalculate untransformed face centres
01638         //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
01639         label vertI = 0;
01640 
01641         forAll(half1Ctrs, i)
01642         {
01643             if (from1To0[i] != -1)
01644             {
01645                 // Write edge between c1 and c0
01646                 //const point& c0 = rawHalf0Ctrs[from1To0[i]];
01647                 const point& c0 = half0Ctrs[from1To0[i]];
01648                 const point& c1 = half1Ctrs[i];
01649                 writeOBJ(ccStr, c0, c1, vertI);
01650             }
01651         }
01652     }
01653 
01654 
01655     if (!matchedAll)
01656     {
01657         SeriousErrorIn
01658         (
01659             "cyclicPolyPatch::order"
01660             "(const primitivePatch&, labelList&, labelList&) const"
01661         )   << "Patch:" << name() << " : "
01662             << "Cannot match vectors to faces on both sides of patch" << endl
01663             << "    Perhaps your faces do not match?"
01664             << " The obj files written contain the current match." << endl
01665             << "    Continuing with incorrect face ordering from now on!"
01666             << endl;
01667 
01668             return false;
01669     }
01670 
01671 
01672     // Set faceMap such that half0 faces get first and corresponding half1
01673     // faces last.
01674     matchAnchors
01675     (
01676         true,                   // report if anchor matching error
01677         pp,
01678         half0ToPatch,
01679         anchors0,
01680         half1ToPatch,
01681         half1Faces,
01682         from1To0,
01683         tols,
01684         faceMap,
01685         rotation
01686     );
01687 
01688     // Return false if no change necessary, true otherwise.
01689 
01690     forAll(faceMap, faceI)
01691     {
01692         if (faceMap[faceI] != faceI || rotation[faceI] != 0)
01693         {
01694             return true;
01695         }
01696     }
01697 
01698     return false;
01699 }
01700 
01701 
01702 void Foam::cyclicPolyPatch::write(Ostream& os) const
01703 {
01704     polyPatch::write(os);
01705     os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
01706     switch (transform_)
01707     {
01708         case ROTATIONAL:
01709         {
01710             os.writeKeyword("transform") << transformTypeNames[ROTATIONAL]
01711                 << token::END_STATEMENT << nl;
01712             os.writeKeyword("rotationAxis") << rotationAxis_
01713                 << token::END_STATEMENT << nl;
01714             os.writeKeyword("rotationCentre") << rotationCentre_
01715                 << token::END_STATEMENT << nl;
01716             break;
01717         }
01718         case TRANSLATIONAL:
01719         {
01720             os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
01721                 << token::END_STATEMENT << nl;
01722             os.writeKeyword("separationVector") << separationVector_
01723                 << token::END_STATEMENT << nl;
01724             break;
01725         }
01726         default:
01727         {
01728             // no additional info to write
01729         }
01730     }
01731 }
01732 
01733 
01734 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines