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].area(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 const label max0I = findMaxArea(pp.points(), half0Faces);
362 const vector n0 = half0Faces[max0I].normal(pp.points());
364 const label max1I = findMaxArea(pp.points(), half1Faces);
365 const vector n1 = half1Faces[max1I].normal(pp.points());
367 if (
mag(n0 & n1) < 1-matchTolerance())
371 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 372 <<
" Detected rotation :" 373 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
400 const pointField& half0Pts = half0.localPoints();
401 const point ctr0(
sum(half0Pts)/half0Pts.size());
404 const pointField& half1Pts = half1.localPoints();
405 const point ctr1(
sum(half1Pts)/half1Pts.size());
409 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 410 <<
" Detected translation :" 411 <<
" n0:" << n0 <<
" n1:" << n1
412 <<
" ctr0:" << ctr0 <<
" ctr1:" << ctr1 <<
endl;
415 half0Ctrs += ctr1 - ctr0;
416 anchors0 += ctr1 - ctr0;
417 ppPoints = pp.points() + ctr1 - ctr0;
425 tols = matchTolerance()*calcFaceTol(half1Faces, pp.points(), half1Ctrs);
429 bool Foam::oldCyclicPolyPatch::matchAnchors
449 forAll(half0ToPatch, half0Facei)
452 label patchFacei = half0ToPatch[half0Facei];
454 faceMap[patchFacei] = half0Facei;
457 rotation[patchFacei] = 0;
460 bool fullMatch =
true;
462 forAll(from1To0, half1Facei)
464 label patchFacei = half1ToPatch[half1Facei];
467 label half0Facei = from1To0[half1Facei];
469 label newFacei = half0Facei + pp.size()/2;
471 faceMap[patchFacei] = newFacei;
477 const point& wantedAnchor = anchors0[half0Facei];
479 rotation[newFacei] = getRotation
482 half1Faces[half1Facei],
487 if (rotation[newFacei] == -1)
493 const face&
f = half1Faces[half1Facei];
495 <<
"Patch:" <<
name() <<
" : " 496 <<
"Cannot find point on face " << f
498 << UIndirectList<point>(pp.points(),
f)()
499 <<
" that matches point " << wantedAnchor
500 <<
" when matching the halves of cyclic patch " <<
name()
502 <<
"Continuing with incorrect face ordering from now on!" 511 Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
518 magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
522 (faceCentres - rotationCentre_) & rotationAxis_
524 axisLen = axisLen -
min(axisLen);
527 magRadSqr + axisLen*axisLen
531 scalar maxMagLenSqr = -great;
532 scalar maxMagRadSqr = -great;
535 if (magLenSqr[i] >= maxMagLenSqr)
537 if (magRadSqr[i] > maxMagRadSqr)
540 maxMagLenSqr = magLenSqr[i];
541 maxMagRadSqr = magRadSqr[i];
548 Info<<
"getConsistentRotationFace(const pointField&)" <<
nl 549 <<
" rotFace = " << rotFace <<
nl 550 <<
" point = " << faceCentres[rotFace] <<
endl;
566 const word& patchType,
573 rotationCentre_(
Zero),
574 separationVector_(
Zero)
584 const word& patchType
590 rotationCentre_(
Zero),
591 separationVector_(
Zero)
593 if (dict.
found(
"neighbourPatch"))
598 ) <<
"Found \"neighbourPatch\" entry when reading cyclic patch " 600 <<
"Is this mesh already with split cyclics?" << endl
601 <<
"If so run a newer version that supports it" 602 <<
", if not comment out the \"neighbourPatch\" entry and re-run" 612 dict.
lookup(
"rotationAxis") >> rotationAxis_;
613 dict.
lookup(
"rotationCentre") >> rotationCentre_;
618 dict.
lookup(
"separationVector") >> separationVector_;
636 featureCos_(pp.featureCos_),
637 rotationAxis_(pp.rotationAxis_),
638 rotationCentre_(pp.rotationCentre_),
639 separationVector_(pp.separationVector_)
653 featureCos_(pp.featureCos_),
654 rotationAxis_(pp.rotationAxis_),
655 rotationCentre_(pp.rotationCentre_),
656 separationVector_(pp.separationVector_)
758 <<
"Size of cyclic " <<
name() <<
" should be a multiple of 2" 762 label halfSize = pp.size()/2;
780 half1ToPatch = half0ToPatch + halfSize;
787 pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
815 Pout<<
"oldCyclicPolyPatch::order : test if already ordered:" 816 << matchedAll <<
endl;
820 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 821 <<
" faces to OBJ file " << nm0 <<
endl;
822 writeOBJ(nm0, half0Faces, ppPoints);
825 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 826 <<
" faces to OBJ file " << nm1 <<
endl;
832 /
"match1_"+
name() +
"_faceCentres.obj" 834 Pout<<
"oldCyclicPolyPatch::order : " 835 <<
"Dumping currently found cyclic match as lines between" 836 <<
" corresponding face centres to file " << ccStr.
name()
850 const point& c0 = half0Ctrs[i];
851 const point&
c1 = half1Ctrs[i];
864 for (
label i = 0; i < halfSize; i++)
866 half0ToPatch[i] = facei++;
867 half1ToPatch[i] = facei++;
899 Pout<<
"oldCyclicPolyPatch::order : test if pairwise ordered:" 900 << matchedAll <<
endl;
903 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 904 <<
" faces to OBJ file " << nm0 <<
endl;
905 writeOBJ(nm0, half0Faces, ppPoints);
908 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 909 <<
" faces to OBJ file " << nm1 <<
endl;
910 writeOBJ(nm1, half1Faces, pp.points());
915 /
"match2_"+
name()+
"_faceCentres.obj" 917 Pout<<
"oldCyclicPolyPatch::order : " 918 <<
"Dumping currently found cyclic match as lines between" 919 <<
" corresponding face centres to file " << ccStr.
name()
927 if (from1To0[i] != -1)
930 const point& c0 = half0Ctrs[from1To0[i]];
931 const point&
c1 = half1Ctrs[i];
951 label matchedFacei = -1;
955 label otherFacei = pFaces[i];
957 if (otherFacei > facei)
965 matchedFacei = otherFacei;
971 if (matchedFacei != -1)
973 half0ToPatch[baffleI] = facei;
974 half1ToPatch[baffleI] = matchedFacei;
979 if (baffleI == halfSize)
1010 Pout<<
"oldCyclicPolyPatch::order : test if baffles:" 1011 << matchedAll <<
endl;
1014 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1015 <<
" faces to OBJ file " << nm0 <<
endl;
1016 writeOBJ(nm0, half0Faces, ppPoints);
1019 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1020 <<
" faces to OBJ file " << nm1 <<
endl;
1021 writeOBJ(nm1, half1Faces, pp.points());
1026 /
"match3_"+
name() +
"_faceCentres.obj" 1028 Pout<<
"oldCyclicPolyPatch::order : " 1029 <<
"Dumping currently found cyclic match as lines between" 1030 <<
" corresponding face centres to file " << ccStr.
name()
1038 if (from1To0[i] != -1)
1041 const point& c0 = half0Ctrs[from1To0[i]];
1042 const point&
c1 = half1Ctrs[i];
1057 bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1069 getCentresAndAnchors
1094 Pout<<
"oldCyclicPolyPatch::order : automatic ordering result:" 1095 << matchedAll <<
endl;
1098 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1099 <<
" faces to OBJ file " << nm0 <<
endl;
1100 writeOBJ(nm0, half0Faces, ppPoints);
1103 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1104 <<
" faces to OBJ file " << nm1 <<
endl;
1105 writeOBJ(nm1, half1Faces, pp.points());
1110 /
"match4_"+
name() +
"_faceCentres.obj" 1112 Pout<<
"oldCyclicPolyPatch::order : " 1113 <<
"Dumping currently found cyclic match as lines between" 1114 <<
" corresponding face centres to file " << ccStr.
name()
1122 if (from1To0[i] != -1)
1125 const point& c0 = half0Ctrs[from1To0[i]];
1126 const point&
c1 = half1Ctrs[i];
1134 if (!matchedAll || debug)
1138 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1139 <<
" faces to OBJ file " << nm0 <<
endl;
1143 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1144 <<
" faces to OBJ file " << nm1 <<
endl;
1150 /
name() +
"_faceCentres.obj" 1152 Pout<<
"oldCyclicPolyPatch::order : " 1153 <<
"Dumping currently found cyclic match as lines between" 1154 <<
" corresponding face centres to file " << ccStr.
name()
1163 if (from1To0[i] != -1)
1167 const point& c0 = half0Ctrs[from1To0[i]];
1168 const point&
c1 = half1Ctrs[i];
1178 <<
"Patch:" <<
name() <<
" : " 1179 <<
"Cannot match vectors to faces on both sides of patch" << endl
1180 <<
" Perhaps your faces do not match?" 1181 <<
" The obj files written contain the current match." << endl
1182 <<
" Continuing with incorrect face ordering from now on!" 1209 if (faceMap[facei] != facei || rotation[facei] != 0)
1242 os.
writeKeyword(
"separationVector") << separationVector_
1253 <<
"Please run foamUpgradeCyclics to convert these old-style" 1254 <<
" 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 occurrences 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.
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.
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]
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.