51 Foam::label Foam::cyclicPolyPatch::findMaxArea
58 scalar maxAreaSqr = -GREAT;
62 scalar areaSqr =
magSqr(faces[faceI].normal(points));
64 if (areaSqr > maxAreaSqr)
74 void Foam::cyclicPolyPatch::calcTransforms()
93 pointField half0Ctrs(calcFaceCentres(half0, half0.points()));
95 scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs));
107 pointField half1Ctrs(calcFaceCentres(half1, half1.points()));
115 Pout<<
"cyclicPolyPatch::calcTransforms : Writing half0"
116 <<
" faces to OBJ file " << nm0 <<
endl;
117 writeOBJ(nm0, half0, half0.points());
120 Pout<<
"cyclicPolyPatch::calcTransforms : Writing half1"
121 <<
" faces to OBJ file " << nm1 <<
endl;
122 writeOBJ(nm1, half1, half1.points());
126 Pout<<
"cyclicPolyPatch::calcTransforms :"
127 <<
" Writing coupled face centres as lines to " << str.
name()
131 const point& p0 = half0Ctrs[i];
132 str <<
"v " << p0.
x() <<
' ' << p0.y() <<
' ' << p0.z() <<
nl;
134 const point& p1 = half1Ctrs[i];
135 str <<
"v " << p1.
x() <<
' ' << p1.y() <<
' ' << p1.z() <<
nl;
137 str <<
"l " << vertI-1 <<
' ' << vertI <<
nl;
144 scalar maxAreaDiff = -GREAT;
145 label maxAreaFacei = -1;
147 for (label facei = 0; facei < size()/2; facei++)
149 half0Normals[facei] = operator[](facei).normal(points);
150 label nbrFacei = facei+size()/2;
151 half1Normals[facei] = operator[](nbrFacei).normal(points);
153 scalar magSf =
mag(half0Normals[facei]);
154 scalar nbrMagSf =
mag(half1Normals[facei]);
155 scalar avSf = (magSf + nbrMagSf)/2.0;
157 if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
162 half0Normals[facei] =
point(1, 0, 0);
163 half1Normals[facei] = half0Normals[facei];
167 scalar areaDiff =
mag(magSf - nbrMagSf)/avSf;
169 if (areaDiff > maxAreaDiff)
171 maxAreaDiff = areaDiff;
172 maxAreaFacei = facei;
179 "cyclicPolyPatch::calcTransforms()"
180 ) <<
"face " << facei <<
" area does not match neighbour "
181 << nbrFacei <<
" by "
183 <<
"% -- possible face ordering problem." << endl
184 <<
"patch:" <<
name()
185 <<
" my area:" << magSf
186 <<
" neighbour area:" << nbrMagSf
187 <<
" matching tolerance:"
190 <<
"Mesh face:" << start()+facei
192 << UIndirectList<point>(
points, operator[](facei))()
194 <<
"Neighbour face:" << start()+nbrFacei
196 << UIndirectList<point>(
points, operator[](nbrFacei))()
198 <<
"Rerun with cyclic debug flag set"
203 half0Normals[facei] /= magSf;
204 half1Normals[facei] /= nbrMagSf;
213 label nbrFacei = maxAreaFacei+size()/2;
214 Pout<<
"cyclicPolyPatch::calcTransforms :"
215 <<
" Max area error:" << 100*maxAreaDiff <<
"% at face:"
216 << maxAreaFacei <<
" at:" << half0Ctrs[maxAreaFacei]
217 <<
" coupled face:" << nbrFacei
218 <<
" at:" << half1Ctrs[maxAreaFacei]
234 label face0 = getConsistentRotationFace(half0Ctrs);
239 (half0Ctrs[face0] - rotationCentre_)
244 (half1Ctrs[face1] - rotationCentre_)
247 n0 /=
mag(n0) + VSMALL;
248 n1 /=
mag(n1) + VSMALL;
252 Pout<<
"cyclicPolyPatch::calcTransforms :"
253 <<
" Specified rotation :"
254 <<
" n0:" << n0 <<
" from face " << half0Ctrs[face0]
255 <<
" and n1:" << n1 <<
" from face " << half1Ctrs[face1]
292 Pout<<
"cyclicPolyPatch::calcTransforms :"
293 <<
" Specified separation vector : "
295 <<
" . Calculated average separation : "
301 const scalar avgTol =
average(half0Tols);
305 ||
mag(separation(0) - separationVector_) > avgTol
310 "cyclicPolyPatch::calcTransforms()"
311 ) <<
"Specified separationVector " << separationVector_
312 <<
" differs from computed separation vector "
314 <<
"This probably means your geometry is not consistent"
315 <<
" with the specified separation and might lead"
316 <<
" to problems." << endl
317 <<
"Continuing with specified separation vector "
318 << separationVector_ << endl
319 <<
"patch:" <<
name()
362 bool Foam::cyclicPolyPatch::getGeometricHalves
373 boolList regionEdge(pp.nEdges(),
false);
377 label nRegionEdges = 0;
381 const labelList& eFaces = edgeFaces[edgeI];
385 if (eFaces.size() == 2)
387 if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
389 regionEdge[edgeI] =
true;
399 patchZones ppZones(pp, regionEdge);
403 Pout<<
"cyclicPolyPatch::getGeometricHalves : "
404 <<
"Found " << nRegionEdges <<
" edges on patch " <<
name()
405 <<
" where the cos of the angle between two connected faces"
406 <<
" was less than " << featureCos_ <<
nl
407 <<
"Patch divided by these and by single sides edges into "
408 << ppZones.nZones() <<
" parts." <<
endl;
415 for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
419 boundaryMesh().
mesh().time().path()
422 Pout<<
"cyclicPolyPatch::getGeometricHalves : Writing zone "
423 << zoneI <<
" face centres to OBJ file " << stream.
name()
430 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
433 nZoneFaces[zoneI] = zoneFaces.size();
438 if (ppZones.nZones() == 2)
447 Pout<<
"cyclicPolyPatch::getGeometricHalves :"
448 <<
" falling back to face-normal comparison" <<
endl;
451 half0ToPatch.setSize(pp.size());
454 half1ToPatch.setSize(pp.size());
457 forAll(faceNormals, faceI)
459 if ((faceNormals[faceI] & faceNormals[0]) > 0)
461 half0ToPatch[n0Faces++] = faceI;
465 half1ToPatch[n1Faces++] = faceI;
468 half0ToPatch.setSize(n0Faces);
469 half1ToPatch.setSize(n1Faces);
473 Pout<<
"cyclicPolyPatch::getGeometricHalves :"
474 <<
" Number of faces per zone:("
475 << n0Faces <<
' ' << n1Faces <<
')' <<
endl;
479 if (half0ToPatch.size() != half1ToPatch.size())
486 Pout<<
"cyclicPolyPatch::getGeometricHalves : Writing half0"
487 <<
" faces to OBJ file " << nm0 <<
endl;
488 writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
491 Pout<<
"cyclicPolyPatch::getGeometricHalves : Writing half1"
492 <<
" faces to OBJ file " << nm1 <<
endl;
493 writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
499 Pout<<
"cyclicPolyPatch::getGeometricHalves : Writing half0"
500 <<
" face centres to OBJ file " << str0.
name() <<
endl;
504 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
508 Pout<<
"cyclicPolyPatch::getGeometricHalves : Writing half1"
509 <<
" face centres to OBJ file " << str1.
name() <<
endl;
512 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
518 "cyclicPolyPatch::getGeometricHalves"
519 "(const primitivePatch&, labelList&, labelList&) const"
520 ) <<
"Patch " <<
name() <<
" gets decomposed in two zones of"
521 <<
"inequal size: " << half0ToPatch.size()
522 <<
" and " << half1ToPatch.size() << endl
523 <<
"This means that the patch is either not two separate regions"
524 <<
" or one region where the angle between the different regions"
525 <<
" is not sufficiently sharp." << endl
526 <<
"Please adapt the featureCos setting." << endl
527 <<
"Continuing with incorrect face ordering from now on!" <<
endl;
541 void Foam::cyclicPolyPatch::getCentresAndAnchors
555 half0Ctrs = calcFaceCentres(half0Faces, pp.points());
556 anchors0 = getAnchorPoints(half0Faces, pp.points());
557 half1Ctrs = calcFaceCentres(half1Faces, pp.points());
563 label face0 = getConsistentRotationFace(half0Ctrs);
564 label face1 = getConsistentRotationFace(half1Ctrs);
566 vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
567 vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
568 n0 /=
mag(n0) + VSMALL;
569 n1 /=
mag(n1) + VSMALL;
573 Pout<<
"cyclicPolyPatch::getCentresAndAnchors :"
574 <<
" Specified rotation :"
575 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
588 half0Ctrs[faceI] - rotationCentre_
595 anchors0[faceI] - rotationCentre_
604 (pp.points() - rotationCentre_)()
617 Pout<<
"cyclicPolyPatch::getCentresAndAnchors :"
618 <<
"Specified translation : " << separationVector_
622 half0Ctrs += separationVector_;
623 anchors0 += separationVector_;
624 ppPoints = pp.points() + separationVector_;
636 label max0I = findMaxArea(pp.points(), half0Faces);
637 vector n0 = half0Faces[max0I].normal(pp.points());
638 n0 /=
mag(n0) + VSMALL;
640 label max1I = findMaxArea(pp.points(), half1Faces);
641 vector n1 = half1Faces[max1I].normal(pp.points());
642 n1 /=
mag(n1) + VSMALL;
648 Pout<<
"cyclicPolyPatch::getCentresAndAnchors :"
649 <<
" Detected rotation :"
650 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
677 const pointField& half0Pts = half0.localPoints();
678 const point ctr0(
sum(half0Pts)/half0Pts.size());
681 const pointField& half1Pts = half1.localPoints();
682 const point ctr1(
sum(half1Pts)/half1Pts.size());
686 Pout<<
"cyclicPolyPatch::getCentresAndAnchors :"
687 <<
" Detected translation :"
688 <<
" n0:" << n0 <<
" n1:" << n1
689 <<
" ctr0:" << ctr0 <<
" ctr1:" << ctr1 <<
endl;
692 half0Ctrs += ctr1 - ctr0;
693 anchors0 += ctr1 - ctr0;
694 ppPoints = pp.points() + ctr1 - ctr0;
702 tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
707 bool Foam::cyclicPolyPatch::matchAnchors
727 forAll(half0ToPatch, half0FaceI)
730 label patchFaceI = half0ToPatch[half0FaceI];
732 faceMap[patchFaceI] = half0FaceI;
735 rotation[patchFaceI] = 0;
738 bool fullMatch =
true;
740 forAll(from1To0, half1FaceI)
742 label patchFaceI = half1ToPatch[half1FaceI];
745 label half0FaceI = from1To0[half1FaceI];
747 label newFaceI = half0FaceI + pp.size()/2;
749 faceMap[patchFaceI] = newFaceI;
755 const point& wantedAnchor = anchors0[half0FaceI];
757 rotation[newFaceI] = getRotation
760 half1Faces[half1FaceI],
765 if (rotation[newFaceI] == -1)
771 const face&
f = half1Faces[half1FaceI];
774 "cyclicPolyPatch::matchAnchors(..)"
775 ) <<
"Patch:" <<
name() <<
" : "
776 <<
"Cannot find point on face " << f
778 << UIndirectList<point>(pp.points(),
f)()
779 <<
" that matches point " << wantedAnchor
780 <<
" when matching the halves of cyclic patch " <<
name()
782 <<
"Continuing with incorrect face ordering from now on!"
791 Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
799 magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
801 label rotFace =
findMax(magRadSqr);
805 Info<<
"getConsistentRotationFace(const pointField&)" <<
nl
806 <<
" rotFace = " << rotFace <<
nl
807 <<
" point = " << faceCentres[rotFace] <<
nl
808 <<
" distance = " <<
Foam::sqrt(magRadSqr[rotFace])
828 coupledPointsPtr_(NULL),
829 coupledEdgesPtr_(NULL),
847 coupledPointsPtr_(NULL),
848 coupledEdgesPtr_(NULL),
855 if (dict.
found(
"neighbourPatch"))
859 "cyclicPolyPatch::cyclicPolyPatch\n"
861 " const word& name,\n"
862 " const dictionary& dict,\n"
863 " const label index,\n"
864 " const polyBoundaryMesh& bm\n"
867 ) <<
"Found \"neighbourPatch\" entry when reading cyclic patch "
869 <<
"Is this mesh already with split cyclics?" << endl
870 <<
"If so run a newer version that supports it"
871 <<
", if not comment out the \"neighbourPatch\" entry and re-run"
877 if (dict.
found(
"transform"))
879 transform_ = transformTypeNames.read(dict.
lookup(
"transform"));
884 dict.
lookup(
"rotationAxis") >> rotationAxis_;
885 dict.
lookup(
"rotationCentre") >> rotationCentre_;
890 dict.
lookup(
"separationVector") >> separationVector_;
909 coupledPointsPtr_(NULL),
910 coupledEdgesPtr_(NULL),
911 featureCos_(pp.featureCos_),
912 transform_(pp.transform_),
913 rotationAxis_(pp.rotationAxis_),
914 rotationCentre_(pp.rotationCentre_),
915 separationVector_(pp.separationVector_)
929 coupledPointsPtr_(NULL),
930 coupledEdgesPtr_(NULL),
931 featureCos_(pp.featureCos_),
932 transform_(pp.transform_),
933 rotationAxis_(pp.rotationAxis_),
934 rotationCentre_(pp.rotationCentre_),
935 separationVector_(pp.separationVector_)
988 if (!coupledPointsPtr_)
998 for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
1000 const face& fA = localFaces()[patchFaceA];
1004 label patchPointA = fA[indexA];
1006 if (coupledPoint[patchPointA] == -1)
1008 const face& fB = localFaces()[patchFaceA + size()/2];
1010 label indexB = (fB.
size() - indexA) % fB.
size();
1013 if (patchPointA != fB[indexB])
1015 coupledPoint[patchPointA] = fB[indexB];
1022 edgeList& connected = *coupledPointsPtr_;
1025 label connectedI = 0;
1029 if (coupledPoint[i] != -1)
1031 connected[connectedI++] =
edge(i, coupledPoint[i]);
1035 connected.
setSize(connectedI);
1042 /
"coupledPoints.obj"
1046 Pout<<
"Writing file " << str.
name() <<
" with coordinates of "
1047 <<
"coupled points" <<
endl;
1051 const point& a =
points()[meshPoints()[connected[i][0]]];
1052 const point&
b =
points()[meshPoints()[connected[i][1]]];
1054 str<<
"v " << a.
x() <<
' ' << a.
y() <<
' ' << a.
z() <<
nl;
1055 str<<
"v " << b.
x() <<
' ' << b.
y() <<
' ' << b.
z() <<
nl;
1058 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1066 return *coupledPointsPtr_;
1072 if (!coupledEdgesPtr_)
1075 const edgeList& pointCouples = coupledPoints();
1081 const edge&
e = pointCouples[i];
1083 aToB.insert(e[0], e[1]);
1089 for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
1091 const labelList& fEdges = faceEdges()[patchFaceA];
1095 label edgeI = fEdges[i];
1097 const edge&
e = edges()[edgeI];
1102 if (fnd0 != aToB.
end())
1105 if (fnd1 != aToB.
end())
1113 coupledEdgesPtr_ =
new edgeList(nEdges()/2);
1114 edgeList& coupledEdges = *coupledEdgesPtr_;
1117 for (label patchFaceB = size()/2; patchFaceB < size(); patchFaceB++)
1119 const labelList& fEdges = faceEdges()[patchFaceB];
1123 label edgeI = fEdges[i];
1125 const edge&
e = edges()[edgeI];
1130 if (iter != edgeMap.
end())
1132 label halfAEdgeI = iter();
1135 if (halfAEdgeI != edgeI)
1137 coupledEdges[coupleI++] =
edge(halfAEdgeI, edgeI);
1141 edgeMap.
erase(iter);
1145 coupledEdges.
setSize(coupleI);
1152 const edge&
e = coupledEdges[i];
1154 if (e[0] == e[1] || e[0] < 0 || e[1] < 0)
1157 <<
"Problem : at position " << i
1158 <<
" illegal couple:" << e
1172 Pout<<
"Writing file " << str.
name() <<
" with centres of "
1173 <<
"coupled edges" <<
endl;
1177 const edge&
e = coupledEdges[i];
1179 const point& a = edges()[e[0]].
centre(localPoints());
1180 const point&
b = edges()[e[1]].centre(localPoints());
1182 str<<
"v " << a.
x() <<
' ' << a.
y() <<
' ' << a.
z() <<
nl;
1183 str<<
"v " << b.
x() <<
' ' << b.
y() <<
' ' << b.
z() <<
nl;
1186 str<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
1194 return *coupledEdgesPtr_;
1228 <<
"Size of cyclic " <<
name() <<
" should be a multiple of 2"
1232 label halfSize = pp.size()/2;
1250 half1ToPatch = half0ToPatch + halfSize;
1257 pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
1259 getCentresAndAnchors
1274 Pout<<
"half0 transformed faceCentres (avg) : "
1276 <<
"half1 untransformed faceCentres (avg) : "
1294 Pout<<
"cyclicPolyPatch::order : test if already ordered:"
1295 << matchedAll <<
endl;
1299 Pout<<
"cyclicPolyPatch::order : Writing half0"
1300 <<
" faces to OBJ file " << nm0 <<
endl;
1301 writeOBJ(nm0, half0Faces, ppPoints);
1304 Pout<<
"cyclicPolyPatch::order : Writing half1"
1305 <<
" faces to OBJ file " << nm1 <<
endl;
1313 /
"match1_"+
name() +
"_faceCentres.obj"
1315 Pout<<
"cyclicPolyPatch::order : "
1316 <<
"Dumping currently found cyclic match as lines between"
1317 <<
" corresponding face centres to file " << ccStr.
name()
1324 const point& c0 = half0Ctrs[from1To0[i]];
1325 const point& c1 = half1Ctrs[i];
1338 for (label i = 0; i < halfSize; i++)
1340 half0ToPatch[i] = faceI++;
1341 half1ToPatch[i] = faceI++;
1348 getCentresAndAnchors
1373 Pout<<
"cyclicPolyPatch::order : test if pairwise ordered:"
1374 << matchedAll <<
endl;
1377 Pout<<
"cyclicPolyPatch::order : Writing half0"
1378 <<
" faces to OBJ file " << nm0 <<
endl;
1379 writeOBJ(nm0, half0Faces, ppPoints);
1382 Pout<<
"cyclicPolyPatch::order : Writing half1"
1383 <<
" faces to OBJ file " << nm1 <<
endl;
1384 writeOBJ(nm1, half1Faces, pp.points());
1389 /
"match2_"+
name()+
"_faceCentres.obj"
1391 Pout<<
"cyclicPolyPatch::order : "
1392 <<
"Dumping currently found cyclic match as lines between"
1393 <<
" corresponding face centres to file " << ccStr.
name()
1402 if (from1To0[i] != -1)
1406 const point& c0 = half0Ctrs[from1To0[i]];
1407 const point& c1 = half1Ctrs[i];
1427 label matchedFaceI = -1;
1431 label otherFaceI = pFaces[i];
1433 if (otherFaceI > faceI)
1441 matchedFaceI = otherFaceI;
1447 if (matchedFaceI != -1)
1449 half0ToPatch[baffleI] = faceI;
1450 half1ToPatch[baffleI] = matchedFaceI;
1455 if (baffleI == halfSize)
1461 getCentresAndAnchors
1486 Pout<<
"cyclicPolyPatch::order : test if baffles:"
1487 << matchedAll <<
endl;
1490 Pout<<
"cyclicPolyPatch::order : Writing half0"
1491 <<
" faces to OBJ file " << nm0 <<
endl;
1492 writeOBJ(nm0, half0Faces, ppPoints);
1495 Pout<<
"cyclicPolyPatch::order : Writing half1"
1496 <<
" faces to OBJ file " << nm1 <<
endl;
1497 writeOBJ(nm1, half1Faces, pp.points());
1502 /
"match3_"+
name() +
"_faceCentres.obj"
1504 Pout<<
"cyclicPolyPatch::order : "
1505 <<
"Dumping currently found cyclic match as lines between"
1506 <<
" corresponding face centres to file " << ccStr.
name()
1515 if (from1To0[i] != -1)
1519 const point& c0 = half0Ctrs[from1To0[i]];
1520 const point& c1 = half1Ctrs[i];
1535 bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1547 getCentresAndAnchors
1572 Pout<<
"cyclicPolyPatch::order : automatic ordering result:"
1573 << matchedAll <<
endl;
1576 Pout<<
"cyclicPolyPatch::order : Writing half0"
1577 <<
" faces to OBJ file " << nm0 <<
endl;
1578 writeOBJ(nm0, half0Faces, ppPoints);
1581 Pout<<
"cyclicPolyPatch::order : Writing half1"
1582 <<
" faces to OBJ file " << nm1 <<
endl;
1583 writeOBJ(nm1, half1Faces, pp.points());
1588 /
"match4_"+
name() +
"_faceCentres.obj"
1590 Pout<<
"cyclicPolyPatch::order : "
1591 <<
"Dumping currently found cyclic match as lines between"
1592 <<
" corresponding face centres to file " << ccStr.
name()
1601 if (from1To0[i] != -1)
1605 const point& c0 = half0Ctrs[from1To0[i]];
1606 const point& c1 = half1Ctrs[i];
1614 if (!matchedAll || debug)
1618 Pout<<
"cyclicPolyPatch::order : Writing half0"
1619 <<
" faces to OBJ file " << nm0 <<
endl;
1623 Pout<<
"cyclicPolyPatch::order : Writing half1"
1624 <<
" faces to OBJ file " << nm1 <<
endl;
1630 /
name() +
"_faceCentres.obj"
1632 Pout<<
"cyclicPolyPatch::order : "
1633 <<
"Dumping currently found cyclic match as lines between"
1634 <<
" corresponding face centres to file " << ccStr.
name()
1643 if (from1To0[i] != -1)
1647 const point& c0 = half0Ctrs[from1To0[i]];
1648 const point& c1 = half1Ctrs[i];
1659 "cyclicPolyPatch::order"
1660 "(const primitivePatch&, labelList&, labelList&) const"
1661 ) <<
"Patch:" <<
name() <<
" : "
1662 <<
"Cannot match vectors to faces on both sides of patch" << endl
1663 <<
" Perhaps your faces do not match?"
1664 <<
" The obj files written contain the current match." << endl
1665 <<
" Continuing with incorrect face ordering from now on!"
1692 if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1710 os.
writeKeyword(
"transform") << transformTypeNames[ROTATIONAL]
1720 os.
writeKeyword(
"transform") << transformTypeNames[TRANSLATIONAL]
1722 os.
writeKeyword(
"separationVector") << separationVector_