59 octantNormal[mXmYmZ] =
vector(-1, -1, -1);
60 octantNormal[pXmYmZ] =
vector(1, -1, -1);
61 octantNormal[mXpYmZ] =
vector(-1, 1, -1);
62 octantNormal[pXpYmZ] =
vector(1, 1, -1);
64 octantNormal[mXmYpZ] =
vector(-1, -1, 1);
65 octantNormal[pXmYpZ] =
vector(1, -1, 1);
66 octantNormal[mXpYpZ] =
vector(-1, 1, 1);
67 octantNormal[pXpYpZ] =
vector(1, 1, 1);
69 octantNormal /=
mag(octantNormal);
84 pn[pointi] = pointNormals[pointi];
89 <<
"Average point normal not visible for point:"
109 visOctant &= ~mXmYmZMask;
110 visOctant &= ~mXmYpZMask;
111 visOctant &= ~mXpYmZMask;
112 visOctant &= ~mXpYpZMask;
114 else if (
n.x() < -SMALL)
117 visOctant &= ~pXmYmZMask;
118 visOctant &= ~pXmYpZMask;
119 visOctant &= ~pXpYmZMask;
120 visOctant &= ~pXpYpZMask;
125 visOctant &= ~mXmYmZMask;
126 visOctant &= ~mXmYpZMask;
127 visOctant &= ~pXmYmZMask;
128 visOctant &= ~pXmYpZMask;
130 else if (
n.y() < -SMALL)
132 visOctant &= ~mXpYmZMask;
133 visOctant &= ~mXpYpZMask;
134 visOctant &= ~pXpYmZMask;
135 visOctant &= ~pXpYpZMask;
140 visOctant &= ~mXmYmZMask;
141 visOctant &= ~mXpYmZMask;
142 visOctant &= ~pXmYmZMask;
143 visOctant &= ~pXpYmZMask;
145 else if (
n.z() < -SMALL)
147 visOctant &= ~mXmYpZMask;
148 visOctant &= ~mXpYpZMask;
149 visOctant &= ~pXmYpZMask;
150 visOctant &= ~pXpYpZMask;
158 for (label octant = 0; octant < 8; octant++)
160 if (visOctant & mask)
174 pn[pointi] = octantNormal[visI];
181 <<
"No visible octant for point:" << pp.
meshPoints()[pointi]
183 <<
"Normal set to " << pn[pointi] <<
endl;
198 return mesh.edges()[edgeI].unitVec(
mesh.points());
208 os <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
220 const point& pt = pts[i];
221 os <<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
nl;
249 os <<
"v " << p1.
x() <<
' ' << p1.
y() <<
' ' << p1.
z() <<
nl;
250 os <<
"v " << p2.
x() <<
' ' << p2.
y() <<
' ' << p2.
z() <<
nl;
251 os <<
"l " << (count + 1) <<
" " << (count + 2) <<
endl;
264 os <<
"v " << p1.
x() <<
' ' << p1.
y() <<
' ' << p1.
z() <<
nl;
266 << (p2.
x() - p1.
x()) <<
' '
267 << (p2.
y() - p1.
y()) <<
' '
268 << (p2.
z() - p1.
z()) <<
endl;
280 for (
const edge&
e : treeBoundBox::edges)
282 os <<
"l " << (
e[0] + 1) <<
' ' << (
e[1] + 1) <<
nl;
298 for (
const label celli : cellLabels)
314 return mesh.edgeCells(edgeI).found(celli);
325 return mesh.faceEdges(facei).found(edgeI);
336 if (
mesh.isInternalFace(facei))
340 (
mesh.faceOwner()[facei] == celli)
341 || (
mesh.faceNeighbour()[facei] == celli)
349 if (
mesh.faceOwner()[facei] == celli)
368 const label edgeI = candidates[i];
370 const edge&
e = edges[edgeI];
372 if ((
e[0] == v0 &&
e[1] == v1) || (
e[0] == v1 &&
e[1] == v0))
394 label edgeI = v0Edges[i];
396 const edge&
e = edges[edgeI];
398 if ((
e.start() == v1) || (
e.end() == v1))
419 label edge0 = f0Edges[f0EdgeI];
423 label edge1 = f1Edges[f1EdgeI];
432 <<
"Faces " << f0 <<
" and " << f1 <<
" do not share an edge"
447 const cell& cFaces =
mesh.cells()[cell0I];
451 label facei = cFaces[cFacei];
455 mesh.isInternalFace(facei)
457 mesh.faceOwner()[facei] == cell1I
458 ||
mesh.faceNeighbour()[facei] == cell1I
468 <<
"No common face for"
469 <<
" cell0I:" << cell0I <<
" faces:" << cFaces
470 <<
" cell1I:" << cell1I <<
" faces:"
471 <<
mesh.cells()[cell1I]
494 label facei = eFaces[eFacei];
511 if ((face0 == -1) || (face1 == -1))
514 <<
"Can not find faces using edge " <<
mesh.edges()[edgeI]
524 const label thisEdgeI,
525 const label thisVertI
528 forAll(edgeLabels, edgeLabelI)
530 label edgeI = edgeLabels[edgeLabelI];
532 if (edgeI != thisEdgeI)
536 if ((
e.start() == thisVertI) || (
e.end() == thisVertI))
544 <<
"Can not find edge in "
546 <<
" connected to edge "
547 << thisEdgeI <<
" with vertices " <<
mesh.edges()[thisEdgeI]
581 const label otherCelli,
585 if (!
mesh.isInternalFace(facei))
588 <<
"Face " << facei <<
" is not internal"
592 label newCelli =
mesh.faceOwner()[facei];
594 if (newCelli == otherCelli)
596 newCelli =
mesh.faceNeighbour()[facei];
606 const label startEdgeI,
607 const label startVertI,
613 label edgeI = startEdgeI;
615 label vertI = startVertI;
617 for (label iter = 0; iter <
nEdges; iter++)
621 vertI =
mesh.edges()[edgeI].otherVertex(vertI);
641 if (dirs[cmpt] == -1)
643 pt[cmpt] = 0.5*(
min[cmpt] +
max[cmpt]);
660 bool isConstrained =
false;
663 if (dirs[cmpt] == -1)
665 isConstrained =
true;
676 if (dirs[cmpt] == -1)
678 pts[i][cmpt] = 0.5*(
min[cmpt] +
max[cmpt]);
695 if (dirs[cmpt] == -1)
710 bool isConstrained =
false;
713 if (dirs[cmpt] == -1)
715 isConstrained =
true;
726 if (dirs[cmpt] == -1)
747 label facei = meshTools::otherFace(
mesh, celli, -1, e0);
750 e1 = meshTools::walkFace(
mesh, facei, e0,
mesh.edges()[e0].end(), 2);
752 facei = meshTools::otherFace(
mesh, celli, facei, e1);
754 e2 = meshTools::walkFace(
mesh, facei, e1,
mesh.edges()[e1].end(), 2);
756 facei = meshTools::otherFace(
mesh, celli, facei, e2);
758 e3 = meshTools::walkFace(
mesh, facei, e2,
mesh.edges()[e2].end(), 2);
766 const label startEdgeI
769 if (!hexMatcher::test(
mesh, celli))
778 label edgeI = startEdgeI;
782 for (label i = 0; i < 3; i++)
785 facei = meshTools::otherFace(
mesh, celli, facei, edgeI);
789 if ((eVec & avgVec) > 0)
798 label vertI =
mesh.edges()[edgeI].end();
800 edgeI = meshTools::walkFace(
mesh, facei, edgeI, vertI, 2);
816 if (!hexMatcher::test(
mesh, celli))
826 scalar maxCos = -GREAT;
829 for (label i = 0; i < 4; i++)
833 label e0 = cEdges[cEdgeI];
835 if (!doneEdges.
found(e0))
839 scalar cosAngle =
mag(avgDir & cutDir);
841 if (cosAngle > maxCos)
861 if (!doneEdges.
found(cEdges[cEdgeI]))
864 <<
"Cell:" << celli <<
" edges:" << cEdges <<
endl
865 <<
"Edge:" << cEdges[cEdgeI] <<
" not yet handled"
873 <<
"Problem : did not find edge aligned with " << cutDir
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
List< Key > toc() const
The table of contents (the keys) in unsorted order.
bool found(const Key &key) const
Return true if hashed entry is found in table.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
A list of faces which address into the list of points.
label nPoints() const
Number of points supporting patch faces.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const Field< point_type > & pointNormals() const
Return point normals for patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & pointFaces() const
Return point-face addressing.
A List with indirect addressing. Like IndirectList but does not store addressing.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
void size(const label n)
Older name for setAddressableSize.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
A cell is defined as a list of faces with extra functionality.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
static constexpr direction nComponents
Number of components in bool is 1.
Mesh consisting of general polyhedral cells.
Cell-face mesh analysis engine.
Standard boundBox with extra functionality for use in octree.
tmp< pointField > points() const
Vertex coordinates. In octant coding.
Representation of a 3D Cartesian coordinate system as a Vector of row vectors.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
#define WarningInFunction
Report a warning using Foam::Warning.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
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.