54 const UList<face>& faces,
62 ctrs[facei] = faces[facei].centre(
points);
71 const UList<face>& faces,
79 anchors[facei] =
points[faces[facei][0]];
86Foam::label Foam::oldCyclicPolyPatch::findMaxArea
93 scalar maxAreaSqr = -GREAT;
97 scalar areaSqr =
magSqr(faces[facei].areaNormal(
points));
99 if (maxAreaSqr < areaSqr)
101 maxAreaSqr = areaSqr;
109bool Foam::oldCyclicPolyPatch::getGeometricHalves
126 boolList regionEdge(pp.nEdges(),
false);
130 label nRegionEdges = 0;
134 const labelList& eFaces = edgeFaces[edgeI];
138 if (eFaces.size() == 2)
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());
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;
286void 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);
436bool 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!"
518Foam::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)
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"
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;
832 writeOBJ(nm1, half1Faces, pp.
points());
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];
857 writeOBJ(ccStr, c0, c1, vertI);
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];
937 writeOBJ(ccStr, c0, c1, vertI);
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];
1048 writeOBJ(ccStr, c0, c1, vertI);
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];
1132 writeOBJ(ccStr, c0, c1, vertI);
1139 if (!matchedAll || debug)
1143 Pout<<
"oldCyclicPolyPatch::order : Writing half0"
1144 <<
" faces to OBJ file " << nm0 <<
endl;
1145 writeOBJ(nm0, half0Faces, pp.
points());
1148 Pout<<
"oldCyclicPolyPatch::order : Writing half1"
1149 <<
" faces to OBJ file " << nm1 <<
endl;
1150 writeOBJ(nm1, half1Faces, pp.
points());
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];
1174 writeOBJ(ccStr, c0, c1, vertI);
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."
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
void setSize(const label n)
Alias for resize()
Output to file stream, using an OSstream.
virtual const fileName & name() const
Read/write access to the name of the stream.
virtual const fileName & name() const
Get the name of the stream.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
A list of faces which address into the list of points.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & pointFaces() const
Return point-face addressing.
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
A List with indirect addressing. Like IndirectList but does not store addressing.
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
virtual transformType transform() const
Type of transform.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
void calcGeometry()
Calculate the geometry for the patches.
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
virtual bool write()
Write the output fields.
virtual void initMovePoints()
Initialise the patches for moving points.
order
Enumeration specifying required accuracy.
void movePoints()
Update for new mesh geometry.
void updateMesh()
Update for new mesh topology.
'old' style cyclic polyPatch with all faces in single patch. Does ordering but cannot be used to run....
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
virtual void initOrder(PstreamBuffers &, const primitivePatch &) const
Initialize ordering for primitivePatch. Does not.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
virtual ~oldCyclicPolyPatch()
const word & name() const noexcept
The patch name.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
A patch is a list of labels that address the faces in the global face list.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Tensor of scalars, i.e. Tensor<scalar>.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
Determine correspondence between points. See below.
#define WarningInFunction
Report a warning using Foam::Warning.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
List< label > labelList
A List of labels.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
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.
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
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]
#define forAll(list, i)
Loop across all elements in list.
static const char *const typeName
The type name used in ensight case files.