diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index 21b03fdaa68b5f088d93c63bb612d07f697f8e6c..6e8eff98c2757d80540ebdefb8838cc2242be6a3 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -43,6 +43,18 @@ namespace Foam addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word); addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary); + + +template<> +const char* NamedEnum<cyclicPolyPatch::transformType, 3>::names[] = +{ + "unknown", + "rotational", + "translational" +}; + +const NamedEnum<cyclicPolyPatch::transformType, 3> + cyclicPolyPatch::transformTypeNames; } @@ -52,6 +64,23 @@ void Foam::cyclicPolyPatch::calcTransforms() { if (size() > 0) { + const pointField& points = this->points(); + + maxI_ = -1; + scalar maxAreaSqr = -GREAT; + + for (label faceI = 0; faceI < size()/2; faceI++) + { + const face& f = operator[](faceI); + scalar areaSqr = magSqr(f.normal(points)); + + if (areaSqr > maxAreaSqr) + { + maxAreaSqr = areaSqr; + maxI_ = faceI; + } + } + primitivePatch half0 ( SubList<face> @@ -59,7 +88,7 @@ void Foam::cyclicPolyPatch::calcTransforms() *this, size()/2 ), - points() + points ); pointField half0Ctrs(calcFaceCentres(half0, half0.points())); scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs)); @@ -72,7 +101,7 @@ void Foam::cyclicPolyPatch::calcTransforms() size()/2, size()/2 ), - points() + points ); pointField half1Ctrs(calcFaceCentres(half1, half1.points())); @@ -97,9 +126,9 @@ void Foam::cyclicPolyPatch::calcTransforms() for (label facei = 0; facei < size()/2; facei++) { - half0Normals[facei] = operator[](facei).normal(points()); + half0Normals[facei] = operator[](facei).normal(points); label nbrFacei = facei+size()/2; - half1Normals[facei] = operator[](nbrFacei).normal(points()); + half1Normals[facei] = operator[](nbrFacei).normal(points); scalar magSf = mag(half0Normals[facei]); scalar nbrMagSf = mag(half1Normals[facei]); @@ -129,11 +158,11 @@ void Foam::cyclicPolyPatch::calcTransforms() << endl << "Mesh face:" << start()+facei << " vertices:" - << IndirectList<point>(points(), operator[](facei))() + << IndirectList<point>(points, operator[](facei))() << endl << "Neighbour face:" << start()+nbrFacei << " vertices:" - << IndirectList<point>(points(), operator[](nbrFacei))() + << IndirectList<point>(points, operator[](nbrFacei))() << endl << "Rerun with cyclic debug flag set" << " for more information." << exit(FatalError); @@ -185,7 +214,7 @@ bool Foam::cyclicPolyPatch::getGeometricHalves { const labelList& eFaces = edgeFaces[edgeI]; - // Check manifold edges for sharp angle. + // Check manifold edges for sharp angle. // (Non-manifold already handled by patchZones) if (eFaces.size() == 2) { @@ -353,16 +382,36 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors ) const { // Get geometric data on both halves. - - vector n0 = half0Faces[0].normal(pp.points()); - n0 /= mag(n0)+VSMALL; - vector n1 = half1Faces[0].normal(pp.points()); - n1 /= mag(n1)+VSMALL; - half0Ctrs = calcFaceCentres(half0Faces, pp.points()); anchors0 = getAnchorPoints(half0Faces, pp.points()); half1Ctrs = calcFaceCentres(half1Faces, pp.points()); + vector n0 = vector::zero; + vector n1 = vector::zero; + switch (transform_) + { + case ROTATIONAL: + { + label face0 = getConsistentRotationFace(half0Ctrs); + label face1 = getConsistentRotationFace(half1Ctrs); + + n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_); + n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_); + n0 /= mag(n0) + VSMALL; + n1 /= mag(n1) + VSMALL; + break; + } + default: + { + // Assumes that cyclic is correctly ordered, so that face[maxI_] + // on each side is equivalent. + n0 = half0Faces[maxI_].normal(pp.points()); + n0 /= mag(n0) + VSMALL; + n1 = half1Faces[maxI_].normal(pp.points()); + n1 /= mag(n1) + VSMALL; + } + } + if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol) { if (debug) @@ -373,7 +422,6 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors // Rotation (around origin) const tensor reverseT(rotationTensor(n0, -n1)); - const tensor forwardT(rotationTensor(-n1, n0)); // Rotation forAll(half0Ctrs, faceI) @@ -495,6 +543,44 @@ bool Foam::cyclicPolyPatch::matchAnchors } +Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace +( + const pointField& faceCentres +) const +{ + const scalarField magRadSqr = + magSqr((faceCentres - rotationCentre_) ^ rotationAxis_); + scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_; + axisLen = axisLen - min(axisLen); + const scalarField magLenSqr = magRadSqr + axisLen*axisLen; + + label rotFace = -1; + scalar maxMagLenSqr = -GREAT; + scalar maxMagRadSqr = -GREAT; + forAll(faceCentres, i) + { + if (magLenSqr[i] >= maxMagLenSqr) + { + if (magRadSqr[i] > maxMagRadSqr) + { + rotFace = i; + maxMagLenSqr = magLenSqr[i]; + maxMagRadSqr = magRadSqr[i]; + } + } + } + + if (debug) + { + Info<< "getConsistentRotationFace(const pointField&)" << nl + << " rotFace = " << rotFace << nl + << " point = " << faceCentres[rotFace] << endl; + } + + return rotFace; +} + + // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * // Foam::cyclicPolyPatch::cyclicPolyPatch @@ -509,7 +595,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(name, size, start, index, bm), coupledPointsPtr_(NULL), coupledEdgesPtr_(NULL), - featureCos_(0.9) + featureCos_(0.9), + maxI_(-1), + transform_(UNKNOWN), + rotationAxis_(vector::zero), + rotationCentre_(point::zero) { calcTransforms(); } @@ -526,12 +616,33 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(name, dict, index, bm), coupledPointsPtr_(NULL), coupledEdgesPtr_(NULL), - featureCos_(0.9) + featureCos_(0.9), + maxI_(-1), + transform_(UNKNOWN), + rotationAxis_(vector::zero), + rotationCentre_(point::zero) { if (dict.found("featureCos")) { dict.lookup("featureCos") >> featureCos_; } + if (dict.found("transform")) + { + transform_ = transformTypeNames.read(dict.lookup("transform")); + switch (transform_) + { + case ROTATIONAL: + { + dict.lookup("rotationAxis") >> rotationAxis_; + dict.lookup("rotationCentre") >> rotationCentre_; + break; + } + default: + { + // no additioanl info required + } + } + } calcTransforms(); } @@ -546,7 +657,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(pp, bm), coupledPointsPtr_(NULL), coupledEdgesPtr_(NULL), - featureCos_(0.9) + featureCos_(pp.featureCos_), + maxI_(pp.maxI_), + transform_(pp.transform_), + rotationAxis_(pp.rotationAxis_), + rotationCentre_(pp.rotationCentre_) { calcTransforms(); } @@ -564,7 +679,11 @@ Foam::cyclicPolyPatch::cyclicPolyPatch coupledPolyPatch(pp, bm, index, newSize, newStart), coupledPointsPtr_(NULL), coupledEdgesPtr_(NULL), - featureCos_(0.9) + featureCos_(pp.featureCos_), + maxI_(pp.maxI_), + transform_(pp.transform_), + rotationAxis_(pp.rotationAxis_), + rotationCentre_(pp.rotationCentre_) { calcTransforms(); } @@ -1031,7 +1150,7 @@ bool Foam::cyclicPolyPatch::order Pout<< "cyclicPolyPatch::order : test if baffles:" << matchedAll << endl; } - + } } @@ -1176,6 +1295,29 @@ void Foam::cyclicPolyPatch::write(Ostream& os) const { polyPatch::write(os); os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl; + switch (transform_) + { + case ROTATIONAL: + { + os.writeKeyword("transform") << transformTypeNames[ROTATIONAL] + << token::END_STATEMENT << nl; + os.writeKeyword("rotationAxis") << rotationAxis_ + << token::END_STATEMENT << nl; + os.writeKeyword("rotationCentre") << rotationCentre_ + << token::END_STATEMENT << nl; + break; + } + case TRANSLATIONAL: + { + os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL] + << token::END_STATEMENT << nl; + break; + } + default: + { + // no additional info to write + } + } } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H index b9bf48aa09210be466df9f8888f768a039693428..186f8d77451c1296d21e629938866761a48263a4 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H @@ -70,6 +70,19 @@ class cyclicPolyPatch : public coupledPolyPatch { +public: + + enum transformType + { + UNKNOWN, + ROTATIONAL, + TRANSLATIONAL + }; + static const NamedEnum<transformType, 3> transformTypeNames; + + +private: + // Private data //- List of edges formed from connected points. e[0] is the point on @@ -85,6 +98,18 @@ class cyclicPolyPatch // Used to split cyclic into halves. scalar featureCos_; + //- Index of largest cell face + label maxI_; + + //- Type of transformation - rotational or translational + transformType transform_; + + //- Axis of rotation for rotational cyclics + vector rotationAxis_; + + //- point on axis of rotation for rotational cyclics + point rotationCentre_; + // Private member functions @@ -131,6 +156,14 @@ class cyclicPolyPatch labelList& rotation ) const; + //- For rotational cases, try to find a unique face on each side + // of the cyclic. + label getConsistentRotationFace + ( + const pointField& faceCentres + ) const; + + protected: // Protected Member functions