diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C index ac7eac71c1..8bc34d2026 100644 --- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C +++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C @@ -300,7 +300,11 @@ label mergePatchFaces const faceZone& fZone = mesh.faceZones()[zoneID]; zoneFlip = fZone.flipMap()[fZone.whichFace(newMasterI)]; } - label patchI = mesh.boundaryMesh().whichPatch(newMasterI); + labelPair patchIDs = polyTopoChange::whichPatch + ( + mesh.boundaryMesh(), + newMasterI + ); Pout<< "Restoring new master face " << newMasterI @@ -316,10 +320,11 @@ label mergePatchFaces own, // owner -1, // neighbour false, // face flip - patchI, // patch for face + patchIDs[0], // patch for face false, // remove from zone zoneID, // zone for face - zoneFlip // face flip in zone + zoneFlip, // face flip in zone + patchIDs[1] // subPatch ) ); @@ -341,9 +346,10 @@ label mergePatchFaces -1, // masterEdgeID, newMasterI, // masterFaceID, false, // flipFaceFlux, - patchI, // patchID, + patchIDs[0], // patchID, zoneID, // zoneID, - zoneFlip // zoneFlip + zoneFlip, // zoneFlip + patchIDs[1] // subPatch ) ); } diff --git a/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.C b/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.C index b54cc404f3..b261ad164a 100644 --- a/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.C +++ b/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "cellSplitter.H" @@ -51,17 +49,12 @@ defineTypeNameAndDebug(cellSplitter, 0); void Foam::cellSplitter::getFaceInfo ( const label faceI, - label& patchID, + labelPair& patchIDs, label& zoneID, label& zoneFlip ) const { - patchID = -1; - - if (!mesh_.isInternalFace(faceI)) - { - patchID = mesh_.boundaryMesh().whichPatch(faceI); - } + patchIDs = polyTopoChange::whichPatch(mesh_.boundaryMesh(), faceI); zoneID = mesh_.faceZones().whichZone(faceI); @@ -172,17 +165,16 @@ void Foam::cellSplitter::setRefinement label anchorPoint = mesh_.cellPoints()[cellI][0]; - label addedPointI = - meshMod.setAction + label addedPointI = meshMod.setAction + ( + polyAddPoint ( - polyAddPoint - ( - iter(), // point - anchorPoint, // master point - -1, // zone for point - true // supports a cell - ) - ); + iter(), // point + anchorPoint, // master point + -1, // zone for point + true // supports a cell + ) + ); addedPoints_.insert(cellI, addedPointI); //Pout<< "Added point " << addedPointI @@ -212,18 +204,17 @@ void Foam::cellSplitter::setRefinement // Add other pyramids for (label i = 1; i < cFaces.size(); i++) { - label addedCellI = - meshMod.setAction + label addedCellI = meshMod.setAction + ( + polyAddCell ( - polyAddCell - ( - -1, // master point - -1, // master edge - -1, // master face - cellI, // master cell - -1 // zone - ) - ); + -1, // master point + -1, // master edge + -1, // master face + cellI, // master cell + -1 // zone + ) + ); newCells[i] = addedCellI; } @@ -310,7 +301,8 @@ void Foam::cellSplitter::setRefinement false, // flux flip -1, // patch for face -1, // zone for face - false // face zone flip + false, // face zone flip + -1 // subPatch ) ); } @@ -356,7 +348,8 @@ void Foam::cellSplitter::setRefinement false, // flux flip -1, // patch for face -1, // zone for face - false // face zone flip + false, // face zone flip + -1 // subPatch ) ); } @@ -411,7 +404,8 @@ void Foam::cellSplitter::setRefinement -1, // patch for face false, // remove from zone -1, // zone for face - false // face zone flip + false, // face zone flip + -1 // subPatch ) ); } @@ -429,7 +423,8 @@ void Foam::cellSplitter::setRefinement -1, // patch for face false, // remove from zone -1, // zone for face - false // face zone flip + false, // face zone flip + -1 // subPatch ) ); } @@ -439,8 +434,9 @@ void Foam::cellSplitter::setRefinement { label newOwn = newOwner(faceI, cellToCells); - label patchID, zoneID, zoneFlip; - getFaceInfo(faceI, patchID, zoneID, zoneFlip); + labelPair patchIDs; + label zoneID, zoneFlip; + getFaceInfo(faceI, patchIDs, zoneID, zoneFlip); meshMod.setAction ( @@ -451,10 +447,11 @@ void Foam::cellSplitter::setRefinement newOwn, // owner -1, // neighbour false, // flux flip - patchID, // patch for face + patchIDs[0], // patch for face false, // remove from zone zoneID, // zone for face - zoneFlip // face zone flip + zoneFlip, // face zone flip + patchIDs[1] ) ); } diff --git a/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.H b/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.H index 2dc6b682a1..8226c458be 100644 --- a/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.H +++ b/applications/utilities/mesh/advanced/modifyMesh/cellSplitter.H @@ -39,6 +39,7 @@ SourceFiles #include "Map.H" #include "edge.H" +#include "labelPair.H" #include "typeInfo.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -73,7 +74,7 @@ class cellSplitter void getFaceInfo ( const label faceI, - label& patchID, + labelPair& patchIDs, label& zoneID, label& zoneFlip ) const; diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L index 46eb810e05..feb48b741c 100644 --- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L +++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L @@ -1124,7 +1124,8 @@ int main(int argc, char *argv[]) false, // flipFaceFlux -1, // patchID faceZonei, // zoneID - false // zoneFlip + false, // zoneFlip + -1 // subPatchID ); // Mark face as being done @@ -1158,7 +1159,8 @@ int main(int argc, char *argv[]) false, // flipFaceFlux patchi, // patchID -1, // zoneID - false // zoneFlip + false, // zoneFlip + -1 // subPatchID ); // For baffles create the opposite face @@ -1175,7 +1177,8 @@ int main(int argc, char *argv[]) false, // flipFaceFlux patchi, // patchID -1, // zoneID - false // zoneFlip + false, // zoneFlip + -1 // subPatchID ); } @@ -1203,13 +1206,14 @@ int main(int argc, char *argv[]) faces[facei], owner[facei], neighbour[facei], - -1, //masterPointID - -1, //masterEdgeID - facei, //masterFace - false, //flipFaceFlux - -1, //patchID - -1, //zoneID - false //zoneFlip + -1, // masterPointID + -1, // masterEdgeID + facei, // masterFace + false, // flipFaceFlux + -1, // patchID + -1, // zoneID + false, // zoneFlip + -1 // subPatchID ); } } diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C index 4552b5eaa3..17fbe68016 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh.C @@ -136,7 +136,8 @@ void Foam::extrude2DMesh::setRefinement false, // flipFaceFlux -1, // patchID zoneID, // zoneID - zoneFlip // zoneFlip + zoneFlip, // zoneFlip + -1 // subPatchID ); } @@ -173,7 +174,8 @@ void Foam::extrude2DMesh::setRefinement false, // flipFaceFlux patchI, // patchID zoneID, // zoneID - zoneFlip // zoneFlip + zoneFlip, // zoneFlip + -1 //?TBD subPatchID ); } } @@ -236,7 +238,8 @@ void Foam::extrude2DMesh::setRefinement false, // flipFaceFlux frontPatchI, // patchID -1, // zoneID - false // zoneFlip + false, // zoneFlip + -1 //?TDB subPatchID ); // Offset to create front face. @@ -255,7 +258,8 @@ void Foam::extrude2DMesh::setRefinement false, // flipFaceFlux frontPatchI, // patchID -1, // zoneID - false // zoneFlip + false, // zoneFlip + -1 //?TDB subPatchID ); } } diff --git a/applications/utilities/mesh/generation/snappyHexMesh/Make/options b/applications/utilities/mesh/generation/snappyHexMesh/Make/options index d5262d8d39..fa3975b68f 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/Make/options +++ b/applications/utilities/mesh/generation/snappyHexMesh/Make/options @@ -9,5 +9,4 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ - -L$(FOAM_MPI_LIBBIN) -lparMetisDecompositionMethod \ -lautoMesh diff --git a/src/Allwmake b/src/Allwmake index 5885810667..341ee2ed6f 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -18,16 +18,16 @@ wmake libso finiteVolume ( cd decompositionAgglomeration && ./Allwmake ) -#wmake libso sampling +wmake libso sampling wmake libso dynamicMesh -#wmake libso dynamicFvMesh -#wmake libso topoChangerFvMesh -#wmake libso fvMotionSolver -#wmake libso engine -# -#wmake libso ODE -#wmake libso randomProcesses +wmake libso dynamicFvMesh +wmake libso topoChangerFvMesh +wmake libso fvMotionSolver +wmake libso engine + +wmake libso ODE +wmake libso randomProcesses ( cd thermophysicalModels && ./Allwmake ) ( cd transportModels && ./Allwmake ) @@ -36,7 +36,7 @@ wmake libso dynamicMesh ( cd postProcessing && ./Allwmake ) ( cd conversion && ./Allwmake ) -#wmake libso autoMesh -#wmake libso errorEstimation +wmake libso autoMesh +wmake libso errorEstimation # ----------------------------------------------------------------- end-of-file diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C index 1a5273aa2d..e2a7e11c2c 100644 --- a/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/constraint/cyclic/cyclicPointPatchField.C @@ -26,6 +26,7 @@ License #include "cyclicPointPatchField.H" #include "Swap.H" +#include "transformField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -124,6 +125,13 @@ cyclicPointPatchField::cyclicPointPatchField template void cyclicPointPatchField::swapAdd(Field& pField) const { + // Get neighbouring pointPatch + const cyclicPointPatch& nbrPatch = cyclicPatch_.neighbPatch(); + + // Get neighbouring patch internal field. Written out since cannot get + // access to neighbouring patch field. + Field nbrPf(pField, nbrPatch.meshPoints()); + Field pf(this->patchInternalField(pField)); const edgeList& pairs = cyclicPatch_.transformPairs(); @@ -132,16 +140,14 @@ void cyclicPointPatchField::swapAdd(Field& pField) const { forAll(pairs, pairi) { - Type tmp = pf[pairs[pairi][0]]; - pf[pairs[pairi][0]] = transform(forwardT(), pf[pairs[pairi][1]]); - pf[pairs[pairi][1]] = transform(reverseT(), tmp); + pf[pairs[pairi][0]] = transform(forwardT(), nbrPf[pairs[pairi][1]]); } } else { forAll(pairs, pairi) { - Swap(pf[pairs[pairi][0]], pf[pairs[pairi][1]]); + pf[pairs[pairi][0]] = nbrPf[pairs[pairi][1]]; } } diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C index 102ace3caa..01cb10c749 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C @@ -144,15 +144,7 @@ tmp > pointPatchField::patchInternalField // get addressing const labelList& meshPoints = patch().meshPoints(); - tmp > tvalues(new Field(meshPoints.size())); - Field& values = tvalues(); - - forAll (meshPoints, pointI) - { - values[pointI] = iF[meshPoints[pointI]]; - } - - return tvalues; + return tmp >(new Field(iF, meshPoints)); } diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H index e03380597a..a00f57b65b 100644 --- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H +++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H @@ -38,6 +38,7 @@ SourceFiles #include "coupledFacePointPatch.H" #include "cyclicPolyPatch.H" +#include "pointBoundaryMesh.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -123,6 +124,14 @@ public: return cyclicPolyPatch_; } + //- Return neighbour point patch + const cyclicPointPatch& neighbPatch() const + { + label patchI = cyclicPolyPatch_.neighbPatchID(); + const pointPatch& pp = this->boundaryMesh()[patchI]; + return refCast(pp); + } + //- Are the cyclic planes parallel bool parallel() const { @@ -145,7 +154,8 @@ public: // Access functions for demand driven data //- Return the set of pairs of points that require transformation - // and/or mapping + // and/or mapping. First index is on this patch, second on the + // neighbour patch. virtual const edgeList& transformPairs() const; }; diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H index 4156682818..f1c3d390b3 100644 --- a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H +++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H @@ -132,7 +132,7 @@ public: // associated with any faces virtual const labelList& loneMeshPoints() const; - //- Return point unit normals. Not impelemented. + //- Return point unit normals. Not implemented. virtual const vectorField& pointNormals() const; }; diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H index a280497bf0..968d34a899 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H @@ -39,6 +39,7 @@ SourceFiles #define coupledPolyPatch_H #include "polyPatch.H" +#include "diagTensorField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -231,20 +232,34 @@ public: return true; } - //- Are the planes separated. - virtual bool separated() const = 0; + //- Transform a patch-based field. + //!! TDB with macros? + virtual void transform(scalarField& l) const = 0; + virtual void transform(vectorField& l) const = 0; + virtual void transform(sphericalTensorField& l) const = 0; + virtual void transform(diagTensorField& l) const = 0; + virtual void transform(symmTensorField& l) const = 0; + virtual void transform(tensorField& l) const = 0; - //- If the planes are separated the separation vector. - virtual const vector& separation() const = 0; + //- Transform a patch-based position + virtual void transformPosition(pointField& l) const = 0; - //- Are the cyclic planes parallel. - virtual bool parallel() const = 0; + // Low level geometric information - //- Return face transformation tensor. - virtual const tensor& forwardT() const = 0; + //- Are the planes separated. + virtual bool separated() const = 0; - //- Return neighbour-cell transformation tensor. - virtual const tensor& reverseT() const = 0; + //- If the planes are separated the separation vector. + virtual const vector& separation() const = 0; + + //- Are the cyclic planes parallel. + virtual bool parallel() const = 0; + + //- Return face transformation tensor. + virtual const tensor& forwardT() const = 0; + + //- Return neighbour-cell transformation tensor. + virtual const tensor& reverseT() const = 0; //- Initialise the calculation of the patch geometry diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index 577bec90ff..c1352c9083 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -35,6 +35,7 @@ License #include "EdgeMap.H" #include "Time.H" #include "diagTensor.H" +#include "transformField.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -250,195 +251,13 @@ Pout<< "cyclicPolyPatch::calcTransforms : calculated transforms for:" } -// Get geometric zones of patch by looking at normals. -// Method 1: any edge with sharpish angle is edge between two halves. -// (this will handle e.g. wedge geometries). -// Also two fully disconnected regions will be handled this way. -// Method 2: sort faces into two halves based on face normal. -bool Foam::cyclicPolyPatch::getGeometricHalves -( - const primitivePatch& pp, - labelList& half0ToPatch, - labelList& half1ToPatch -) const -{ - // Calculate normals - const vectorField& faceNormals = pp.faceNormals(); - - // Find edges with sharp angles. - boolList regionEdge(pp.nEdges(), false); - - const labelListList& edgeFaces = pp.edgeFaces(); - - label nRegionEdges = 0; - - forAll(edgeFaces, edgeI) - { - const labelList& eFaces = edgeFaces[edgeI]; - - // Check manifold edges for sharp angle. - // (Non-manifold already handled by patchZones) - if (eFaces.size() == 2) - { - if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_) - { - regionEdge[edgeI] = true; - - nRegionEdges++; - } - } - } - - - // For every face determine zone it is connected to (without crossing - // any regionEdge) - patchZones ppZones(pp, regionEdge); - - if (debug) - { - Pout<< "cyclicPolyPatch::getGeometricHalves : " - << "Found " << nRegionEdges << " edges on patch " << name() - << " where the cos of the angle between two connected faces" - << " was less than " << featureCos_ << nl - << "Patch divided by these and by single sides edges into " - << ppZones.nZones() << " parts." << endl; - - - // Dumping zones to obj files. - - labelList nZoneFaces(ppZones.nZones()); - - for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++) - { - OFstream stream - ( - boundaryMesh().mesh().time().path() - /name()+"_zone_"+Foam::name(zoneI)+".obj" - ); - Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone " - << zoneI << " face centres to OBJ file " << stream.name() - << endl; - - labelList zoneFaces(findIndices(ppZones, zoneI)); - - forAll(zoneFaces, i) - { - writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points())); - } - - nZoneFaces[zoneI] = zoneFaces.size(); - } - } - - - if (ppZones.nZones() == 2) - { - half0ToPatch = findIndices(ppZones, 0); - half1ToPatch = findIndices(ppZones, 1); - } - else - { - if (debug) - { - Pout<< "cyclicPolyPatch::getGeometricHalves :" - << " falling back to face-normal comparison" << endl; - } - label n0Faces = 0; - half0ToPatch.setSize(pp.size()); - - label n1Faces = 0; - half1ToPatch.setSize(pp.size()); - - // Compare to face 0 normal. - forAll(faceNormals, faceI) - { - if ((faceNormals[faceI] & faceNormals[0]) > 0) - { - half0ToPatch[n0Faces++] = faceI; - } - else - { - half1ToPatch[n1Faces++] = faceI; - } - } - half0ToPatch.setSize(n0Faces); - half1ToPatch.setSize(n1Faces); - - if (debug) - { - Pout<< "cyclicPolyPatch::getGeometricHalves :" - << " Number of faces per zone:(" - << n0Faces << ' ' << n1Faces << ')' << endl; - } - } - - if (half0ToPatch.size() != half1ToPatch.size()) - { - fileName casePath(boundaryMesh().mesh().time().path()); - - // Dump halves - { - fileName nm0(casePath/name()+"_half0_faces.obj"); - Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0" - << " faces to OBJ file " << nm0 << endl; - writeOBJ(nm0, IndirectList(pp, half0ToPatch)(), pp.points()); - - fileName nm1(casePath/name()+"_half1_faces.obj"); - Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1" - << " faces to OBJ file " << nm1 << endl; - writeOBJ(nm1, IndirectList(pp, half1ToPatch)(), pp.points()); - } - - // Dump face centres - { - OFstream str0(casePath/name()+"_half0.obj"); - Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0" - << " face centres to OBJ file " << str0.name() << endl; - - forAll(half0ToPatch, i) - { - writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points())); - } - - OFstream str1(casePath/name()+"_half1.obj"); - Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1" - << " face centres to OBJ file " << str1.name() << endl; - forAll(half1ToPatch, i) - { - writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points())); - } - } - - SeriousErrorIn - ( - "cyclicPolyPatch::getGeometricHalves" - "(const primitivePatch&, labelList&, labelList&) const" - ) << "Patch " << name() << " gets decomposed in two zones of" - << "inequal size: " << half0ToPatch.size() - << " and " << half1ToPatch.size() << endl - << "This means that the patch is either not two separate regions" - << " or one region where the angle between the different regions" - << " is not sufficiently sharp." << endl - << "Please adapt the featureCos setting." << endl - << "Continuing with incorrect face ordering from now on!" << endl; - - return false; - } - else - { - return true; - } -} - - // Given a split of faces into left and right half calculate the centres // and anchor points. Transform the left points so they align with the // right ones. void Foam::cyclicPolyPatch::getCentresAndAnchors ( - const primitivePatch& pp, - const faceList& half0Faces, - const faceList& half1Faces, + const primitivePatch& pp0, + const primitivePatch& pp1, pointField& half0Ctrs, pointField& half1Ctrs, @@ -447,9 +266,9 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors ) const { // Get geometric data on both halves. - half0Ctrs = calcFaceCentres(half0Faces, pp.points()); - anchors0 = getAnchorPoints(half0Faces, pp.points()); - half1Ctrs = calcFaceCentres(half1Faces, pp.points()); + half0Ctrs = calcFaceCentres(pp0, pp0.points()); + anchors0 = getAnchorPoints(pp0, pp0.points()); + half1Ctrs = calcFaceCentres(pp1, pp1.points()); vector n0 = vector::zero; vector n1 = vector::zero; @@ -473,12 +292,12 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors // Determine the face with max area on both halves. These // two faces are used to determine the transformation tensors - label max0I = findMaxArea(pp.points(), half0Faces); - n0 = half0Faces[max0I].normal(pp.points()); + label max0I = findMaxArea(pp0.points(), pp0); + n0 = pp0[max0I].normal(pp0.points()); n0 /= mag(n0) + VSMALL; - label max1I = findMaxArea(pp.points(), half1Faces); - n1 = half1Faces[max1I].normal(pp.points()); + label max1I = findMaxArea(pp1.points(), pp1); + n1 = pp1[max1I].normal(pp1.points()); n1 /= mag(n1) + VSMALL; } } @@ -505,13 +324,8 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors { // Parallel translation. Get average of all used points. - primitiveFacePatch half0(half0Faces, pp.points()); - const pointField& half0Pts = half0.localPoints(); - const point ctr0(sum(half0Pts)/half0Pts.size()); - - primitiveFacePatch half1(half1Faces, pp.points()); - const pointField& half1Pts = half1.localPoints(); - const point ctr1(sum(half1Pts)/half1Pts.size()); + const point ctr0(sum(pp0.localPoints())/pp0.nPoints()); + const point ctr1(sum(pp1.localPoints())/pp1.nPoints()); if (debug) { @@ -525,92 +339,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors } // Calculate typical distance per face - tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs); -} - - -// Calculates faceMap and rotation. Returns true if all ok. -bool Foam::cyclicPolyPatch::matchAnchors -( - const bool report, - const primitivePatch& pp, - const labelList& half0ToPatch, - const pointField& anchors0, - - const labelList& half1ToPatch, - const faceList& half1Faces, - const labelList& from1To0, - - const scalarField& tols, - - labelList& faceMap, - labelList& rotation -) const -{ - // Set faceMap such that half0 faces get first and corresponding half1 - // faces last. - - forAll(half0ToPatch, half0FaceI) - { - // Label in original patch - label patchFaceI = half0ToPatch[half0FaceI]; - - faceMap[patchFaceI] = half0FaceI; - - // No rotation - rotation[patchFaceI] = 0; - } - - bool fullMatch = true; - - forAll(from1To0, half1FaceI) - { - label patchFaceI = half1ToPatch[half1FaceI]; - - // This face has to match the corresponding one on half0. - label half0FaceI = from1To0[half1FaceI]; - - label newFaceI = half0FaceI + pp.size()/2; - - faceMap[patchFaceI] = newFaceI; - - // Rotate patchFaceI such that its f[0] aligns with that of - // the corresponding face - // (which after shuffling will be at position half0FaceI) - - const point& wantedAnchor = anchors0[half0FaceI]; - - rotation[newFaceI] = getRotation - ( - pp.points(), - half1Faces[half1FaceI], - wantedAnchor, - tols[half1FaceI] - ); - - if (rotation[newFaceI] == -1) - { - fullMatch = false; - - if (report) - { - const face& f = half1Faces[half1FaceI]; - SeriousErrorIn - ( - "cyclicPolyPatch::matchAnchors(..)" - ) << "Patch:" << name() << " : " - << "Cannot find point on face " << f - << " with vertices:" - << IndirectList(pp.points(), f)() - << " that matches point " << wantedAnchor - << " when matching the halves of cyclic patch " << name() - << endl - << "Continuing with incorrect face ordering from now on!" - << endl; - } - } - } - return fullMatch; + tols = calcFaceTol(pp1, pp1.points(), half1Ctrs); } @@ -666,12 +395,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(name, size, start, index, bm), neighbPatchName_(""), neighbPatchID_(-1), - coupledPointsPtr_(NULL), - coupledEdgesPtr_(NULL), - featureCos_(0.9), transform_(UNKNOWN), rotationAxis_(vector::zero), rotationCentre_(point::zero), + coupledPointsPtr_(NULL), + coupledEdgesPtr_(NULL), separated_(false), separation_(vector::zero), parallel_(true), @@ -694,20 +422,17 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(name, dict, index, bm), neighbPatchName_(dict.lookup("neighbourPatch")), neighbPatchID_(-1), - coupledPointsPtr_(NULL), - coupledEdgesPtr_(NULL), - featureCos_(0.9), transform_(UNKNOWN), rotationAxis_(vector::zero), rotationCentre_(point::zero), + coupledPointsPtr_(NULL), + coupledEdgesPtr_(NULL), separated_(false), separation_(vector::zero), parallel_(true), forwardT_(I), reverseT_(I) { - dict.readIfPresent("featureCos", featureCos_); - if (dict.found("transform")) { transform_ = transformTypeNames.read(dict.lookup("transform")); @@ -740,12 +465,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(pp, bm), neighbPatchName_(pp.neighbPatchName()), neighbPatchID_(-1), - coupledPointsPtr_(NULL), - coupledEdgesPtr_(NULL), - featureCos_(pp.featureCos_), transform_(pp.transform_), rotationAxis_(pp.rotationAxis_), rotationCentre_(pp.rotationCentre_), + coupledPointsPtr_(NULL), + coupledEdgesPtr_(NULL), separated_(false), separation_(vector::zero), parallel_(true), @@ -770,12 +494,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(pp, bm, index, newSize, newStart), neighbPatchName_(neighbPatchName), neighbPatchID_(-1), - coupledPointsPtr_(NULL), - coupledEdgesPtr_(NULL), - featureCos_(pp.featureCos_), transform_(pp.transform_), rotationAxis_(pp.rotationAxis_), - rotationCentre_(pp.rotationCentre_) + rotationCentre_(pp.rotationCentre_), + coupledPointsPtr_(NULL), + coupledEdgesPtr_(NULL) { // Neighbour patch might not be valid yet so no transformation // calculation possible. @@ -794,6 +517,19 @@ Foam::cyclicPolyPatch::~cyclicPolyPatch() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +void Foam::cyclicPolyPatch::transformPosition(pointField& l) const +{ + if (!parallel()) + { + Foam::transform(reverseT_, l); + } + else if (separated()) + { + l += separation_; + } +} + + void Foam::cyclicPolyPatch::initGeometry() { polyPatch::initGeometry(); @@ -872,89 +608,96 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const { if (!coupledPointsPtr_) { - WarningIn("cyclicPolyPatch::coupledPoints() const") - << "to be implemented." << endl; + const faceList& nbrLocalFaces = neighbPatch().localFaces(); + const labelList& nbrMeshPoints = neighbPatch().meshPoints(); - coupledPointsPtr_ = new edgeList(); + // Now all we know is that relative face index in *this is same + // as coupled face in nbrPatch and also that the 0th vertex + // corresponds. + + // From local point to nbrPatch or -1. + labelList coupledPoint(nPoints(), -1); + + forAll(*this, patchFaceI) + { + const face& fA = localFaces()[patchFaceI]; + + forAll(fA, indexA) + { + label patchPointA = fA[indexA]; + + if (coupledPoint[patchPointA] == -1) + { + const face& fB = nbrLocalFaces[patchFaceI]; + + label indexB = (fB.size() - indexA) % fB.size(); + + // Filter out points on wedge axis + if (patchPointA != fB[indexB]) + { + coupledPoint[patchPointA] = fB[indexB]; + } + } + } + } + + coupledPointsPtr_ = new edgeList(nPoints()); + edgeList& connected = *coupledPointsPtr_; + + // Extract coupled points. + label connectedI = 0; + + forAll(coupledPoint, i) + { + if (coupledPoint[i] != -1) + { + connected[connectedI++] = edge(i, coupledPoint[i]); + } + } + + connected.setSize(connectedI); + + if (debug) + { + OFstream str + ( + boundaryMesh().mesh().time().path() + /name() + "_coupledPoints.obj" + ); + label vertI = 0; + + Pout<< "Writing file " << str.name() << " with coordinates of " + << "coupled points" << endl; + + forAll(connected, i) + { + const point& a = points()[meshPoints()[connected[i][0]]]; + const point& b = points()[nbrMeshPoints[connected[i][1]]]; + + str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl; + str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl; + vertI += 2; + + str<< "l " << vertI-1 << ' ' << vertI << nl; + } + } + + // Remove any addressing calculated for the coupled edges calculation + const_cast + ( + static_cast + ( + *this + ) + ).clearOut(); + const_cast + ( + static_cast + ( + neighbPatch() + ) + ).clearOut(); } - - -// // Look at cyclic patch as two halves, A and B. -// // Now all we know is that relative face index in halfA is same -// // as coupled face in halfB and also that the 0th vertex -// // corresponds. -// -// // From halfA point to halfB or -1. -// labelList coupledPoint(nPoints(), -1); -// -// for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++) -// { -// const face& fA = localFaces()[patchFaceA]; -// -// forAll(fA, indexA) -// { -// label patchPointA = fA[indexA]; -// -// if (coupledPoint[patchPointA] == -1) -// { -// const face& fB = localFaces()[patchFaceA + size()/2]; -// -// label indexB = (fB.size() - indexA) % fB.size(); -// -// // Filter out points on wedge axis -// if (patchPointA != fB[indexB]) -// { -// coupledPoint[patchPointA] = fB[indexB]; -// } -// } -// } -// } -// -// coupledPointsPtr_ = new edgeList(nPoints()); -// edgeList& connected = *coupledPointsPtr_; -// -// // Extract coupled points. -// label connectedI = 0; -// -// forAll(coupledPoint, i) -// { -// if (coupledPoint[i] != -1) -// { -// connected[connectedI++] = edge(i, coupledPoint[i]); -// } -// } -// -// connected.setSize(connectedI); -// -// if (debug) -// { -// OFstream str -// ( -// boundaryMesh().mesh().time().path() -// /"coupledPoints.obj" -// ); -// label vertI = 0; -// -// Pout<< "Writing file " << str.name() << " with coordinates of " -// << "coupled points" << endl; -// -// forAll(connected, i) -// { -// const point& a = points()[meshPoints()[connected[i][0]]]; -// const point& b = points()[meshPoints()[connected[i][1]]]; -// -// str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl; -// str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl; -// vertI += 2; -// -// str<< "l " << vertI-1 << ' ' << vertI << nl; -// } -// } -// -// // Remove any addressing calculated for the coupled edges calculation -// const_cast(static_cast(*this)) -// .clearOut(); -// } return *coupledPointsPtr_; } @@ -963,139 +706,164 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const { if (!coupledEdgesPtr_) { - WarningIn("cyclicPolyPatch::coupledEdges() const") - << "to be implemented." << endl; + const edgeList& pointCouples = coupledPoints(); - coupledEdgesPtr_ = new edgeList(); + // Build map from points on *this (A) to points on neighbourpatch (B) + Map