54 const UList<face>& faces,
62 ctrs[facei] = faces[facei].centre(
points);
71 const UList<face>& faces,
79 anchors[facei] =
points[faces[facei][0]];
86 Foam::label Foam::oldCyclicPolyPatch::findMaxArea
93 scalar maxAreaSqr = -GREAT;
97 scalar areaSqr =
magSqr(faces[facei].areaNormal(
points));
99 if (maxAreaSqr < areaSqr)
101 maxAreaSqr = areaSqr;
109 bool Foam::oldCyclicPolyPatch::getGeometricHalves
126 boolList regionEdge(pp.nEdges(),
false);
130 label nRegionEdges = 0;
134 const labelList& eFaces = edgeFaces[edgeI];
138 if (eFaces.size() == 2)
140 if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
142 regionEdge[edgeI] =
true;
152 patchZones ppZones(pp, regionEdge);
156 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : " 157 <<
"Found " << nRegionEdges <<
" edges on patch " <<
name()
158 <<
" where the cos of the angle between two connected faces" 159 <<
" was less than " << featureCos_ <<
nl 160 <<
"Patch divided by these and by single sides edges into " 161 << ppZones.nZones() <<
" parts." <<
endl;
168 for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
175 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing zone " 176 << zoneI <<
" face centres to OBJ file " << stream.
name()
183 writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
186 nZoneFaces[zoneI] = zoneFaces.size();
191 if (ppZones.nZones() == 2)
200 Pout<<
"oldCyclicPolyPatch::getGeometricHalves :" 201 <<
" falling back to face-normal comparison" <<
endl;
204 half0ToPatch.setSize(pp.size());
207 half1ToPatch.setSize(pp.size());
210 forAll(faceNormals, facei)
212 if ((faceNormals[facei] & faceNormals[0]) > 0)
214 half0ToPatch[n0Faces++] = facei;
218 half1ToPatch[n1Faces++] = facei;
221 half0ToPatch.setSize(n0Faces);
222 half1ToPatch.setSize(n1Faces);
226 Pout<<
"oldCyclicPolyPatch::getGeometricHalves :" 227 <<
" Number of faces per zone:(" 228 << n0Faces <<
' ' << n1Faces <<
')' <<
endl;
232 if (half0ToPatch.size() != half1ToPatch.size())
234 fileName casePath(boundaryMesh().
mesh().time().
path());
238 fileName nm0(casePath/
name()+
"_half0_faces.obj");
239 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half0" 240 <<
" faces to OBJ file " << nm0 <<
endl;
241 writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
243 fileName nm1(casePath/
name()+
"_half1_faces.obj");
244 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half1" 245 <<
" faces to OBJ file " << nm1 <<
endl;
246 writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
251 OFstream str0(casePath/
name()+
"_half0.obj");
252 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half0" 253 <<
" face centres to OBJ file " << str0.
name() <<
endl;
257 writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
260 OFstream str1(casePath/
name()+
"_half1.obj");
261 Pout<<
"oldCyclicPolyPatch::getGeometricHalves : Writing half1" 262 <<
" face centres to OBJ file " << str1.
name() <<
endl;
265 writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
270 <<
"Patch " <<
name() <<
" gets decomposed in two zones of" 271 <<
"inequal size: " << half0ToPatch.size()
272 <<
" and " << half1ToPatch.size() <<
endl 273 <<
"This means that the patch is either not two separate regions" 274 <<
" or one region where the angle between the different regions" 275 <<
" is not sufficiently sharp." <<
endl 276 <<
"Please adapt the featureCos setting." <<
endl 277 <<
"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);
314 (half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_
320 (half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_
325 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 326 <<
" Specified rotation :" 327 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
368 label max0I = findMaxArea(pp.points(), half0Faces);
369 const vector n0 = half0Faces[max0I].unitNormal(pp.points());
371 label max1I = findMaxArea(pp.points(), half1Faces);
372 const vector n1 = half1Faces[max1I].unitNormal(pp.points());
374 if (
mag(n0 & n1) < 1-matchTolerance())
378 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 379 <<
" Detected rotation :" 380 <<
" n0:" << n0 <<
" n1:" << n1 <<
endl;
407 const pointField& half0Pts = half0.localPoints();
408 const point ctr0(
sum(half0Pts)/half0Pts.size());
411 const pointField& half1Pts = half1.localPoints();
412 const point ctr1(
sum(half1Pts)/half1Pts.size());
416 Pout<<
"oldCyclicPolyPatch::getCentresAndAnchors :" 417 <<
" Detected translation :" 418 <<
" n0:" << n0 <<
" n1:" << n1
419 <<
" ctr0:" << ctr0 <<
" ctr1:" << ctr1 <<
endl;
422 half0Ctrs += ctr1 - ctr0;
423 anchors0 += ctr1 - ctr0;
424 ppPoints = pp.points() + ctr1 - ctr0;
432 tols = matchTolerance()*calcFaceTol(half1Faces, pp.points(), half1Ctrs);
436 bool Foam::oldCyclicPolyPatch::matchAnchors
456 forAll(half0ToPatch, half0Facei)
459 label patchFacei = half0ToPatch[half0Facei];
461 faceMap[patchFacei] = half0Facei;
464 rotation[patchFacei] = 0;
467 bool fullMatch =
true;
469 forAll(from1To0, half1Facei)
471 label patchFacei = half1ToPatch[half1Facei];
474 label half0Facei = from1To0[half1Facei];
476 label newFacei = half0Facei + pp.size()/2;
478 faceMap[patchFacei] = newFacei;
484 const point& wantedAnchor = anchors0[half0Facei];
486 rotation[newFacei] = getRotation
489 half1Faces[half1Facei],
494 if (rotation[newFacei] == -1)
500 const face&
f = half1Faces[half1Facei];
502 <<
"Patch:" <<
name() <<
" : " 503 <<
"Cannot find point on face " <<
f 505 << UIndirectList<point>(pp.points(),
f)
506 <<
" that matches point " << wantedAnchor
507 <<
" when matching the halves of cyclic patch " <<
name()
509 <<
"Continuing with incorrect face ordering from now on!" 518 Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
525 magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
529 (faceCentres - rotationCentre_) & rotationAxis_
531 axisLen = axisLen -
min(axisLen);
534 magRadSqr + axisLen*axisLen
538 scalar maxMagLenSqr = -GREAT;
539 scalar maxMagRadSqr = -GREAT;
542 if (magLenSqr[i] >= maxMagLenSqr)
544 if (magRadSqr[i] > maxMagRadSqr)
547 maxMagLenSqr = magLenSqr[i];
548 maxMagRadSqr = magRadSqr[i];
555 Info<<
"getConsistentRotationFace(const pointField&)" <<
nl 556 <<
" rotFace = " << rotFace <<
nl 557 <<
" point = " << faceCentres[rotFace] <<
endl;
573 const word& patchType,
580 rotationCentre_(
Zero),
581 separationVector_(
Zero)
591 const word& patchType
597 rotationCentre_(
Zero),
598 separationVector_(
Zero)
600 if (
dict.found(
"neighbourPatch"))
603 <<
"Found \"neighbourPatch\" entry when reading cyclic patch " 605 <<
"Is this mesh already with split cyclics?" <<
endl 606 <<
"If so run a newer version that supports it" 607 <<
", if not comment out the \"neighbourPatch\" entry and re-run" 611 dict.readIfPresent(
"featureCos", featureCos_);
617 dict.readEntry(
"rotationAxis", rotationAxis_);
618 dict.readEntry(
"rotationCentre", rotationCentre_);
623 dict.readEntry(
"separationVector", separationVector_);
641 featureCos_(pp.featureCos_),
642 rotationAxis_(pp.rotationAxis_),
643 rotationCentre_(pp.rotationCentre_),
644 separationVector_(pp.separationVector_)
658 featureCos_(pp.featureCos_),
659 rotationAxis_(pp.rotationAxis_),
660 rotationCentre_(pp.rotationCentre_),
661 separationVector_(pp.separationVector_)
763 <<
"Size of cyclic " <<
name() <<
" should be a multiple of 2" 767 label halfSize = pp.size()/2;
785 half1ToPatch = half0ToPatch + halfSize;
792 pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
820 Pout<<
"oldCyclicPolyPatch::order : test if already ordered:" 821 << matchedAll <<
endl;
825 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 826 <<
" faces to OBJ file " << nm0 <<
endl;
827 writeOBJ(nm0, half0Faces, ppPoints);
830 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 831 <<
" faces to OBJ file " << nm1 <<
endl;
837 /
"match1_"+
name() +
"_faceCentres.obj" 839 Pout<<
"oldCyclicPolyPatch::order : " 840 <<
"Dumping currently found cyclic match as lines between" 841 <<
" corresponding face centres to file " << ccStr.
name()
855 const point& c0 = half0Ctrs[i];
856 const point&
c1 = half1Ctrs[i];
869 for (label i = 0; i < halfSize; i++)
871 half0ToPatch[i] = facei++;
872 half1ToPatch[i] = facei++;
904 Pout<<
"oldCyclicPolyPatch::order : test if pairwise ordered:" 905 << matchedAll <<
endl;
908 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 909 <<
" faces to OBJ file " << nm0 <<
endl;
910 writeOBJ(nm0, half0Faces, ppPoints);
913 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 914 <<
" faces to OBJ file " << nm1 <<
endl;
915 writeOBJ(nm1, half1Faces, pp.points());
920 /
"match2_"+
name()+
"_faceCentres.obj" 922 Pout<<
"oldCyclicPolyPatch::order : " 923 <<
"Dumping currently found cyclic match as lines between" 924 <<
" corresponding face centres to file " << ccStr.
name()
932 if (from1To0[i] != -1)
935 const point& c0 = half0Ctrs[from1To0[i]];
936 const point&
c1 = half1Ctrs[i];
956 label matchedFacei = -1;
960 label otherFacei =
pFaces[i];
962 if (otherFacei > facei)
970 matchedFacei = otherFacei;
976 if (matchedFacei != -1)
978 half0ToPatch[baffleI] = facei;
979 half1ToPatch[baffleI] = matchedFacei;
984 if (baffleI == halfSize)
1015 Pout<<
"oldCyclicPolyPatch::order : test if baffles:" 1016 << matchedAll <<
endl;
1019 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1020 <<
" faces to OBJ file " << nm0 <<
endl;
1021 writeOBJ(nm0, half0Faces, ppPoints);
1024 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1025 <<
" faces to OBJ file " << nm1 <<
endl;
1026 writeOBJ(nm1, half1Faces, pp.points());
1031 /
"match3_"+
name() +
"_faceCentres.obj" 1033 Pout<<
"oldCyclicPolyPatch::order : " 1034 <<
"Dumping currently found cyclic match as lines between" 1035 <<
" corresponding face centres to file " << ccStr.
name()
1043 if (from1To0[i] != -1)
1046 const point& c0 = half0Ctrs[from1To0[i]];
1047 const point&
c1 = half1Ctrs[i];
1062 bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1074 getCentresAndAnchors
1099 Pout<<
"oldCyclicPolyPatch::order : automatic ordering result:" 1100 << matchedAll <<
endl;
1103 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1104 <<
" faces to OBJ file " << nm0 <<
endl;
1105 writeOBJ(nm0, half0Faces, ppPoints);
1108 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1109 <<
" faces to OBJ file " << nm1 <<
endl;
1110 writeOBJ(nm1, half1Faces, pp.points());
1115 /
"match4_"+
name() +
"_faceCentres.obj" 1117 Pout<<
"oldCyclicPolyPatch::order : " 1118 <<
"Dumping currently found cyclic match as lines between" 1119 <<
" corresponding face centres to file " << ccStr.
name()
1127 if (from1To0[i] != -1)
1130 const point& c0 = half0Ctrs[from1To0[i]];
1131 const point&
c1 = half1Ctrs[i];
1139 if (!matchedAll ||
debug)
1143 Pout<<
"oldCyclicPolyPatch::order : Writing half0" 1144 <<
" faces to OBJ file " << nm0 <<
endl;
1148 Pout<<
"oldCyclicPolyPatch::order : Writing half1" 1149 <<
" faces to OBJ file " << nm1 <<
endl;
1155 /
name() +
"_faceCentres.obj" 1157 Pout<<
"oldCyclicPolyPatch::order : " 1158 <<
"Dumping currently found cyclic match as lines between" 1159 <<
" corresponding face centres to file " << ccStr.
name()
1168 if (from1To0[i] != -1)
1172 const point& c0 = half0Ctrs[from1To0[i]];
1173 const point&
c1 = half1Ctrs[i];
1183 <<
"Patch:" <<
name() <<
" : " 1184 <<
"Cannot match vectors to faces on both sides of patch" <<
endl 1185 <<
" Perhaps your faces do not match?" 1186 <<
" The obj files written contain the current match." <<
endl 1187 <<
" Continuing with incorrect face ordering from now on!" 1214 if (
faceMap[facei] != facei || rotation[facei] != 0)
1254 <<
"Please run foamUpgradeCyclics to convert these old-style" 1255 <<
" cyclics into two separate cyclics patches."
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
A List of labelList.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
A class for handling file names.
virtual bool order(PstreamBuffers &, const primitivePatch &, labelList &faceMap, labelList &rotation) const
Return new ordering for primitivePatch.
void write(Ostream &os) const
Write (physicalType, inGroups) dictionary entries (without surrounding braces)
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Output to file stream, using an OSstream.
constexpr char nl
The newline '\n' character (0x0a)
List< face > faceList
A List of faces.
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...
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
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.
virtual void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
Determine correspondence between points. See below.
Macros for easy insertion into run-time selection tables.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
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.
#define forAll(list, i)
Loop across all elements in list.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
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.
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
const Field< point_type > & points() const noexcept
Return reference to global points.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
virtual void initOrder(PstreamBuffers &, const primitivePatch &) const
Initialize ordering for primitivePatch. Does not.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
const labelListList & pointFaces() const
Return point-face addressing.
const word & name() const noexcept
The patch name.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
defineTypeNameAndDebug(combustionModel, 0)
Buffers for inter-processor communications streams (UOPstream, UIPstream).
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::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].reset(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]
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
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.
A List with indirect addressing. Like IndirectList but does not store addressing. ...
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
messageStream Info
Information stream (stdout output on master, null elsewhere)
virtual transformType transform() const
Type of transform.
Field< vector > vectorField
Specialisation of Field<T> for vector.
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
List< label > labelList
A List of labels.
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
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.
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)