52 const UList<face>& faces,
60 ctrs[facei] = faces[facei].centre(points);
69 const UList<face>& faces,
77 anchors[facei] = points[faces[facei][0]];
91 scalar maxAreaSqr = -GREAT;
95 scalar areaSqr =
magSqr(faces[facei].normal(points));
97 if (areaSqr > maxAreaSqr)
107 bool Foam::oldCyclicPolyPatch::getGeometricHalves
124 boolList regionEdge(pp.nEdges(),
false);
128 label nRegionEdges = 0;
132 const labelList& eFaces = edgeFaces[edgeI];
136 if (eFaces.size() == 2)
138 if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
140 regionEdge[edgeI] =
true;
150 patchZones ppZones(pp, regionEdge);
154 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : " 155 <<
"Found " << nRegionEdges <<
" edges on patch " <<
name()
156 <<
" where the cos of the angle between two connected faces" 157 <<
" was less than " << featureCos_ <<
nl 158 <<
"Patch divided by these and by single sides edges into " 159 << ppZones.nZones() <<
" parts." <<
endl;
166 for (
label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
173 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing zone " 174 << zoneI <<
" face centres to OBJ file " << stream.
name()
181 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
184 nZoneFaces[zoneI] = zoneFaces.size();
189 if (ppZones.nZones() == 2)
198 Pout<<
"oldCyclicPolyPatch::getGeometricHalves :" 199 <<
" falling back to face-normal comparison" <<
endl;
202 half0ToPatch.setSize(pp.size());
205 half1ToPatch.setSize(pp.size());
208 forAll(faceNormals, facei)
210 if ((faceNormals[facei] & faceNormals[0]) > 0)
212 half0ToPatch[n0Faces++] = facei;
216 half1ToPatch[n1Faces++] = facei;
219 half0ToPatch.setSize(n0Faces);
220 half1ToPatch.setSize(n1Faces);
224 Pout<<
"oldCyclicPolyPatch::getGeometricHalves :" 225 <<
" Number of faces per zone:(" 226 << n0Faces <<
' ' << n1Faces <<
')' <<
endl;
230 if (half0ToPatch.size() != half1ToPatch.size())
232 fileName casePath(boundaryMesh().
mesh().time().
path());
236 fileName nm0(casePath/
name()+
"_half0_faces.obj");
237 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half0" 238 <<
" faces to OBJ file " << nm0 <<
endl;
239 writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
241 fileName nm1(casePath/
name()+
"_half1_faces.obj");
242 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half1" 243 <<
" faces to OBJ file " << nm1 <<
endl;
244 writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
249 OFstream str0(casePath/
name()+
"_half0.obj");
250 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half0" 251 <<
" face centres to OBJ file " << str0.
name() <<
endl;
255 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
258 OFstream str1(casePath/
name()+
"_half1.obj");
259 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half1" 260 <<
" face centres to OBJ file " << str1.
name() <<
endl;
263 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
268 <<
"Patch " <<
name() <<
" gets decomposed in two zones of" 269 <<
"inequal size: " << half0ToPatch.size()
270 <<
" and " << half1ToPatch.size() << endl
271 <<
"This means that the patch is either not two separate regions" 272 <<
" or one region where the angle between the different regions" 273 <<
" is not sufficiently sharp." << endl
274 <<
"Please adapt the featureCos setting." << endl
275 <<
"Continuing with incorrect face ordering from now on!" <<
endl;
286 void Foam::oldCyclicPolyPatch::getCentresAndAnchors
300 half0Ctrs = calcFaceCentres(half0Faces, pp.points());
301 anchors0 = getAnchorPoints(half0Faces, pp.points());
302 half1Ctrs = calcFaceCentres(half1Faces, pp.points());
308 label face0 = getConsistentRotationFace(half0Ctrs);
309 label face1 = getConsistentRotationFace(half1Ctrs);
311 vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
312 vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
313 n0 /=
mag(n0) + VSMALL;
314 n1 /=
mag(n1) + VSMALL;
318 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 319 <<
" Specified rotation :" 320 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
361 label max0I = findMaxArea(pp.points(), half0Faces);
362 vector n0 = half0Faces[max0I].normal(pp.points());
363 n0 /=
mag(n0) + VSMALL;
365 label max1I = findMaxArea(pp.points(), half1Faces);
366 vector n1 = half1Faces[max1I].normal(pp.points());
367 n1 /=
mag(n1) + VSMALL;
369 if (
mag(n0 & n1) < 1-matchTolerance())
373 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 374 <<
" Detected rotation :" 375 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
402 const pointField& half0Pts = half0.localPoints();
403 const point ctr0(
sum(half0Pts)/half0Pts.size());
406 const pointField& half1Pts = half1.localPoints();
407 const point ctr1(
sum(half1Pts)/half1Pts.size());
411 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 412 <<
" Detected translation :" 413 <<
" n0:" << n0 <<
" n1:" << n1
414 <<
" ctr0:" << ctr0 <<
" ctr1:" << ctr1 <<
endl;
417 half0Ctrs += ctr1 - ctr0;
418 anchors0 += ctr1 - ctr0;
419 ppPoints = pp.points() + ctr1 - ctr0;
427 tols = matchTolerance()*calcFaceTol(half1Faces, pp.points(), half1Ctrs);
431 bool Foam::oldCyclicPolyPatch::matchAnchors
451 forAll(half0ToPatch, half0Facei)
454 label patchFacei = half0ToPatch[half0Facei];
456 faceMap[patchFacei] = half0Facei;
459 rotation[patchFacei] = 0;
462 bool fullMatch =
true;
464 forAll(from1To0, half1Facei)
466 label patchFacei = half1ToPatch[half1Facei];
469 label half0Facei = from1To0[half1Facei];
471 label newFacei = half0Facei + pp.size()/2;
473 faceMap[patchFacei] = newFacei;
479 const point& wantedAnchor = anchors0[half0Facei];
481 rotation[newFacei] = getRotation
484 half1Faces[half1Facei],
489 if (rotation[newFacei] == -1)
495 const face&
f = half1Faces[half1Facei];
497 <<
"Patch:" <<
name() <<
" : " 498 <<
"Cannot find point on face " << f
500 << UIndirectList<point>(pp.points(),
f)()
501 <<
" that matches point " << wantedAnchor
502 <<
" when matching the halves of cyclic patch " <<
name()
504 <<
"Continuing with incorrect face ordering from now on!" 513 Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
520 magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
524 (faceCentres - rotationCentre_) & rotationAxis_
526 axisLen = axisLen -
min(axisLen);
529 magRadSqr + axisLen*axisLen
533 scalar maxMagLenSqr = -GREAT;
534 scalar maxMagRadSqr = -GREAT;
537 if (magLenSqr[i] >= maxMagLenSqr)
539 if (magRadSqr[i] > maxMagRadSqr)
542 maxMagLenSqr = magLenSqr[i];
543 maxMagRadSqr = magRadSqr[i];
550 Info<<
"getConsistentRotationFace(const pointField&)" <<
nl 551 <<
" rotFace = " << rotFace <<
nl 552 <<
" point = " << faceCentres[rotFace] <<
endl;
568 const word& patchType,
575 rotationCentre_(
Zero),
576 separationVector_(
Zero)
586 const word& patchType
592 rotationCentre_(
Zero),
593 separationVector_(
Zero)
595 if (dict.
found(
"neighbourPatch"))
600 ) <<
"Found \"neighbourPatch\" entry when reading cyclic patch " 602 <<
"Is this mesh already with split cyclics?" << endl
603 <<
"If so run a newer version that supports it" 604 <<
", if not comment out the \"neighbourPatch\" entry and re-run" 614 dict.
lookup(
"rotationAxis") >> rotationAxis_;
615 dict.
lookup(
"rotationCentre") >> rotationCentre_;
620 dict.
lookup(
"separationVector") >> separationVector_;
638 featureCos_(pp.featureCos_),
639 rotationAxis_(pp.rotationAxis_),
640 rotationCentre_(pp.rotationCentre_),
641 separationVector_(pp.separationVector_)
655 featureCos_(pp.featureCos_),
656 rotationAxis_(pp.rotationAxis_),
657 rotationCentre_(pp.rotationCentre_),
658 separationVector_(pp.separationVector_)
760 <<
"Size of cyclic " <<
name() <<
" should be a multiple of 2" 764 label halfSize = pp.size()/2;
782 half1ToPatch = half0ToPatch + halfSize;
789 pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
817 Pout<<
"oldCyclicPolyPatch::order : test if already ordered:" 818 << matchedAll <<
endl;
822 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 823 <<
" faces to OBJ file " << nm0 <<
endl;
824 writeOBJ(nm0, half0Faces, ppPoints);
827 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 828 <<
" faces to OBJ file " << nm1 <<
endl;
834 /
"match1_"+
name() +
"_faceCentres.obj" 836 Pout<<
"oldCyclicPolyPatch::order : " 837 <<
"Dumping currently found cyclic match as lines between" 838 <<
" corresponding face centres to file " << ccStr.
name()
852 const point& c0 = half0Ctrs[i];
853 const point&
c1 = half1Ctrs[i];
866 for (
label i = 0; i < halfSize; i++)
868 half0ToPatch[i] = facei++;
869 half1ToPatch[i] = facei++;
901 Pout<<
"oldCyclicPolyPatch::order : test if pairwise ordered:" 902 << matchedAll <<
endl;
905 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 906 <<
" faces to OBJ file " << nm0 <<
endl;
907 writeOBJ(nm0, half0Faces, ppPoints);
910 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 911 <<
" faces to OBJ file " << nm1 <<
endl;
912 writeOBJ(nm1, half1Faces, pp.points());
917 /
"match2_"+
name()+
"_faceCentres.obj" 919 Pout<<
"oldCyclicPolyPatch::order : " 920 <<
"Dumping currently found cyclic match as lines between" 921 <<
" corresponding face centres to file " << ccStr.
name()
929 if (from1To0[i] != -1)
932 const point& c0 = half0Ctrs[from1To0[i]];
933 const point&
c1 = half1Ctrs[i];
953 label matchedFacei = -1;
957 label otherFacei = pFaces[i];
959 if (otherFacei > facei)
967 matchedFacei = otherFacei;
973 if (matchedFacei != -1)
975 half0ToPatch[baffleI] = facei;
976 half1ToPatch[baffleI] = matchedFacei;
981 if (baffleI == halfSize)
1012 Pout<<
"oldCyclicPolyPatch::order : test if baffles:" 1013 << matchedAll <<
endl;
1016 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1017 <<
" faces to OBJ file " << nm0 <<
endl;
1018 writeOBJ(nm0, half0Faces, ppPoints);
1021 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1022 <<
" faces to OBJ file " << nm1 <<
endl;
1023 writeOBJ(nm1, half1Faces, pp.points());
1028 /
"match3_"+
name() +
"_faceCentres.obj" 1030 Pout<<
"oldCyclicPolyPatch::order : " 1031 <<
"Dumping currently found cyclic match as lines between" 1032 <<
" corresponding face centres to file " << ccStr.
name()
1040 if (from1To0[i] != -1)
1043 const point& c0 = half0Ctrs[from1To0[i]];
1044 const point&
c1 = half1Ctrs[i];
1059 bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1071 getCentresAndAnchors
1096 Pout<<
"oldCyclicPolyPatch::order : automatic ordering result:" 1097 << matchedAll <<
endl;
1100 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1101 <<
" faces to OBJ file " << nm0 <<
endl;
1102 writeOBJ(nm0, half0Faces, ppPoints);
1105 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1106 <<
" faces to OBJ file " << nm1 <<
endl;
1107 writeOBJ(nm1, half1Faces, pp.points());
1112 /
"match4_"+
name() +
"_faceCentres.obj" 1114 Pout<<
"oldCyclicPolyPatch::order : " 1115 <<
"Dumping currently found cyclic match as lines between" 1116 <<
" corresponding face centres to file " << ccStr.
name()
1124 if (from1To0[i] != -1)
1127 const point& c0 = half0Ctrs[from1To0[i]];
1128 const point&
c1 = half1Ctrs[i];
1136 if (!matchedAll || debug)
1140 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1141 <<
" faces to OBJ file " << nm0 <<
endl;
1145 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1146 <<
" faces to OBJ file " << nm1 <<
endl;
1152 /
name() +
"_faceCentres.obj" 1154 Pout<<
"oldCyclicPolyPatch::order : " 1155 <<
"Dumping currently found cyclic match as lines between" 1156 <<
" corresponding face centres to file " << ccStr.
name()
1165 if (from1To0[i] != -1)
1169 const point& c0 = half0Ctrs[from1To0[i]];
1170 const point&
c1 = half1Ctrs[i];
1180 <<
"Patch:" <<
name() <<
" : " 1181 <<
"Cannot match vectors to faces on both sides of patch" << endl
1182 <<
" Perhaps your faces do not match?" 1183 <<
" The obj files written contain the current match." << endl
1184 <<
" Continuing with incorrect face ordering from now on!" 1211 if (faceMap[facei] != facei || rotation[facei] != 0)
1244 os.
writeKeyword(
"separationVector") << separationVector_
1255 <<
"Please run foamUpgradeCyclics to convert these old-style" 1256 <<
" cyclics into two separate cyclics patches."
virtual const fileName & name() const
Return the name of the stream.
List< labelList > labelListList
A List of labelList.
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
virtual bool order(PstreamBuffers &, const primitivePatch &, labelList &faceMap, labelList &rotation) const
Return new ordering for primitivePatch.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
A list of keyword definitions, which are a keyword followed by any number of values (e...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Ostream & endl(Ostream &os)
Add newline and flush stream.
'old' style cyclic polyPatch with all faces in single patch. Does ordering but cannot be used to run...
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
Vector< scalar > vector
A scalar version of the templated Vector.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurences of given element. Linear search.
virtual void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
Determine correspondence between points. See below.
Macros for easy insertion into run-time selection tables.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual void movePoints(PstreamBuffers &, const pointField &p)
Correct patches after moving points.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
List< bool > boolList
Bool container classes.
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
virtual ~oldCyclicPolyPatch()
virtual void movePoints(PstreamBuffers &, const pointField &)
Correct patches after moving points.
A list of faces which address into the list of points.
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
vectorField pointField
pointField is a vectorField.
A class for handling words, derived from string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
const Field< PointType > & points() const
Return reference to global points.
List< label > labelList
A List of labels.
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
errorManip< error > abort(error &err)
virtual void initOrder(PstreamBuffers &, const primitivePatch &) const
Initialize ordering for primitivePatch. Does not.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
const labelListList & pointFaces() const
Return point-face addressing.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
defineTypeNameAndDebug(combustionModel, 0)
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
PrimitivePatch< face, List, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
word name(const complex &)
Return a string representation of a complex.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
void write(Ostream &) const
Write patchIdentifier as a dictionary.
void setSize(const label)
Reset size of List.
Template functions to aid in the implementation of demand driven data.
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
prefixOSstream Pout(cout, "Pout")
A List with indirect addressing.
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
Tensor< scalar > tensor
Tensor of scalars.
dimensionSet transform(const dimensionSet &)
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
oldCyclicPolyPatch(const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm, const word &patchType, const transformType transform=UNKNOWN)
Construct from components.