50 createEdgePointGroupByCirculating(feMesh, edHit, pts);
54 label edgeI = edHit.
index();
63 createExternalEdgePointGroup(feMesh, edHit, pts);
68 createInternalEdgePointGroup(feMesh, edHit, pts);
73 createFlatEdgePointGroup(feMesh, edHit, pts);
78 createOpenEdgePointGroup(feMesh, edHit, pts);
83 createMultipleEdgePointGroup(feMesh, edHit, pts);
86 case extendedFeatureEdgeMesh::NONE:
95bool Foam::conformalVoronoiMesh::meshableRegion
127bool Foam::conformalVoronoiMesh::regionIsInside
146 const bool meshableRegionA = meshableRegion(sideA, volTypeA);
147 const bool meshableRegionB = meshableRegion(sideB, volTypeB);
149 if (meshableRegionA == meshableRegionB)
151 return meshableRegionA;
160void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
171 const label edgeI = edHit.
index();
173 scalar ppDist = pointPairDistance(edgePt);
193 bool addedMasterPreviously =
false;
194 label initialRegion = -1;
198 const sideVolumeType volType = normalVolumeTypes[circ.curr()];
199 const sideVolumeType nextVolType = normalVolumeTypes[circ.next()];
201 const vector& normal = feNormals[circ.curr()];
202 const vector& nextNormal = feNormals[circ.next()];
204 vector normalDir = (normal ^ edDir);
205 normalDir *= circNormalDirs.curr()/
mag(normalDir);
207 vector nextNormalDir = (nextNormal ^ edDir);
208 nextNormalDir *= circNormalDirs.next()/
mag(nextNormalDir);
226 ((normalDir ^ nextNormalDir) & edDir) < SMALL
227 ||
mag(masterPtVec) < SMALL
231 addedMasterPreviously =
false;
236 &&
mag((normal & nextNormal) - 1) < SMALL
239 const vector n = 0.5*(normal + nextNormal);
241 const vector s = ppDist*(edDir ^
n);
245 createBafflePointPair(ppDist, edgePt +
s,
n,
true, pts);
246 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
251 <<
"Failed to insert flat/open edge as volType is "
265 const Foam::point masterPt = edgePt + ppDist*masterPtVec;
279 if (
mag((normalDir & nextNormalDir) - 1) < SMALL)
285 vector s = ppDist*(edDir ^ normal);
287 plane facePlane(edgePt, normal);
305 <<
"Faces are parallel but master point is not inside"
310 if (!addedMasterPreviously)
312 if (initialRegion == -1)
314 initialRegion = circ.nRotations();
317 addedMasterPreviously =
true;
319 masterPoints.
insert(circ.nRotations(), masterPt);
328 masterPointReflectionsPrev.
insert
331 plane(edgePt, normal)
334 masterPointReflectionsNext.
insert
337 plane(edgePt, nextNormal)
340 else if (addedMasterPreviously)
342 addedMasterPreviously =
true;
344 masterPointReflectionsNext.
erase(circ.nRotations() - 1);
351 plane p(masterPoints[circ.nRotations() - 1], normalDir);
354 scalar cutPoint =
p.normalIntersect(r);
359 edgePt + cutPoint*(masterPt - edgePt)
370 masterPointReflectionsNext.
insert
373 plane(edgePt, nextNormal)
379 masterPoints.
size() > 1
381 && circ.nRotations() == circ.size() - 1
384 if (initialRegion == 0)
386 plane p(masterPoints[initialRegion], nextNormalDir);
389 scalar cutPoint =
p.normalIntersect(r);
391 masterPoints[circ.nRotations()] =
392 edgePt + cutPoint*(masterPt - edgePt);
397 masterPointReflectionsPrev.
erase(initialRegion);
398 masterPointReflectionsNext.
erase(circ.nRotations());
416 const vertexType ptType = masterPointsTypes[iter.key()];
423 const vertexType reflectedPtType =
430 if (masterPointReflectionsPrev.
found(iter.key()))
433 masterPointReflectionsPrev[iter.key()].mirror(pt);
439 pts.
append(
Vb(reflectedPt, reflectedPtType));
442 if (masterPointReflectionsNext.
found(iter.key()))
445 masterPointReflectionsNext[iter.key()].mirror(pt);
451 pts.
append(
Vb(reflectedPt, reflectedPtType));
459void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
468 scalar ppDist = pointPairDistance(edgePt);
476 const vector& nA = feNormals[edNormalIs[0]];
477 const vector& nB = feNormals[edNormalIs[1]];
480 normalVolumeTypes[edNormalIs[0]];
483 normalVolumeTypes[edNormalIs[1]];
493 vector refVec((nA + nB)/(1 + (nA & nB)));
498 ppDist *= 5.0/
mag(refVec);
516 if (!geometryToConformTo_.inside(refPt))
526 vertexCount() + pts.
size(),
541 vertexCount() + pts.
size(),
557 vertexCount() + pts.
size(),
567 ptPairs_.addPointPair
569 pts[pts.
size() - 3].index(),
570 pts[pts.
size() - 1].index()
573 ptPairs_.addPointPair
575 pts[pts.
size() - 3].index(),
576 pts[pts.
size() - 2].index()
581void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
590 scalar ppDist = pointPairDistance(edgePt);
598 const vector& nA = feNormals[edNormalIs[0]];
599 const vector& nB = feNormals[edNormalIs[1]];
602 normalVolumeTypes[edNormalIs[0]];
615 vector refVec((nA + nB)/(1 + (nA & nB)));
620 ppDist *= 5.0/
mag(refVec);
637 Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
640 Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
642 Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
648 int nQuads = int(totalAngle/foamyHexMeshControls().maxQuadAngle()) + 1;
652 int nAddPoints =
min(
max(nQuads - 2, 0), 2);
678 !geometryToConformTo_.inside(reflectedA)
679 || !geometryToConformTo_.inside(reflectedB)
691 vertexCount() + pts.
size(),
703 vertexCount() + pts.
size(),
715 vertexCount() + pts.
size(),
725 ptPairs_.addPointPair
727 pts[pts.
size() - 2].index(),
728 pts[pts.
size() - 1].index()
731 ptPairs_.addPointPair
733 pts[pts.
size() - 3].index(),
734 pts[pts.
size() - 1].index()
746 vertexCount() + pts.
size(),
752 else if (nAddPoints == 2)
760 vertexCount() + pts.
size(),
772 vertexCount() + pts.
size(),
781void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
790 const scalar ppDist = pointPairDistance(edgePt);
798 const vector& nA = feNormals[edNormalIs[0]];
799 const vector& nB = feNormals[edNormalIs[1]];
803 const vector n = 0.5*(nA + nB);
812 createPointPair(ppDist, edgePt +
s, -
n,
true, pts);
813 createPointPair(ppDist, edgePt -
s, -
n,
true, pts);
817 createBafflePointPair(ppDist, edgePt +
s,
n,
true, pts);
818 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
822 createPointPair(ppDist, edgePt +
s,
n,
true, pts);
823 createPointPair(ppDist, edgePt -
s,
n,
true, pts);
828void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
838 const scalar ppDist = pointPairDistance(edgePt);
843 if (edNormalIs.
size() == 1)
850 const vector&
n = feNormals[edNormalIs[0]];
854 plane facePlane(edgePt,
n);
856 const label initialPtsSize = pts.
size();
860 !geometryToConformTo_.inside(edgePt)
866 createBafflePointPair(ppDist, edgePt -
s,
n,
true, pts);
867 createBafflePointPair(ppDist, edgePt +
s,
n,
false, pts);
869 for (label ptI = initialPtsSize; ptI < pts.
size(); ++ptI)
876 Info<<
"NOT INSERTING OPEN EDGE POINT GROUP WITH MORE THAN 1 "
877 <<
"EDGE NORMAL, NOT IMPLEMENTED" <<
endl;
882void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
893 const scalar ppDist = pointPairDistance(edgePt);
906 forAll(edNormalIs, edgeNormalI)
909 normalVolumeTypes[edNormalIs[edgeNormalI]];
911 nNormalTypes[sType]++;
916 label masterEdgeNormalIndex = -1;
918 forAll(edNormalIs, edgeNormalI)
921 normalVolumeTypes[edNormalIs[edgeNormalI]];
925 masterEdgeNormalIndex = edgeNormalI;
930 const vector&
n = feNormals[edNormalIs[masterEdgeNormalIndex]];
932 label nDir = normalDirs[masterEdgeNormalIndex];
935 (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
936 normalDir *= nDir/
mag(normalDir);
941 plane plane3(edgePt, normalDir);
951 vertexCount() + pts.
size(),
961 vertexCount() + pts.
size(),
967 ptPairs_.addPointPair
969 pts[pts.
size() - 2].index(),
970 pts[pts.
size() - 1].index()
978 vertexCount() + pts.
size(),
984 ptPairs_.addPointPair
986 pts[pts.
size() - 3].index(),
987 pts[pts.
size() - 1].index()
995 vertexCount() + pts.
size(),
1001 ptPairs_.addPointPair
1003 pts[pts.
size() - 3].index(),
1004 pts[pts.
size() - 1].index()
1007 ptPairs_.addPointPair
1009 pts[pts.
size() - 2].index(),
1010 pts[pts.
size() - 1].index()
1019 label masterEdgeNormalIndex = -1;
1021 forAll(edNormalIs, edgeNormalI)
1024 normalVolumeTypes[edNormalIs[edgeNormalI]];
1028 masterEdgeNormalIndex = edgeNormalI;
1033 const vector&
n = feNormals[edNormalIs[masterEdgeNormalIndex]];
1035 label nDir = normalDirs[masterEdgeNormalIndex];
1038 (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
1039 normalDir *= nDir/
mag(normalDir);
1041 const label nextNormalI =
1042 (masterEdgeNormalIndex + 1) % edNormalIs.size();
1043 if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
1048 Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*
n;
1049 Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*
n;
1051 plane plane3(edgePt, normalDir);
1061 vertexCount() + pts.
size(),
1071 vertexCount() + pts.
size(),
1077 ptPairs_.addPointPair
1079 pts[pts.
size() - 2].index(),
1080 pts[pts.
size() - 1].index()
1088 vertexCount() + pts.
size(),
1094 ptPairs_.addPointPair
1096 pts[pts.
size() - 3].index(),
1097 pts[pts.
size() - 1].index()
1105 vertexCount() + pts.
size(),
1111 ptPairs_.addPointPair
1113 pts[pts.
size() - 3].index(),
1114 pts[pts.
size() - 1].index()
1137void Foam::conformalVoronoiMesh::insertFeaturePoints(
bool distribute)
1140 <<
"Inserting feature points" <<
endl;
1142 const label preFeaturePointSize(number_of_vertices());
1146 ftPtConformer_.distribute(decomposition());
1149 const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
1155 ftPtConformer_.reIndexPointPairs(oldToNewIndices);
1157 label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
1160 Info<<
" Inserted " << nFeatureVertices <<
" feature vertices" <<
endl;
CGAL::indexedVertex< K > Vb
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Like Foam::Circulator, but with const-access iterators.
Map< label > insertPoints(const List< Vb > &vertices, const bool reIndex)
Insert the list of vertices (calls rangeInsertWithInfo)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
bool found(const Key &key) const
Return true if hashed entry is found in table.
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
label size() const noexcept
The number of elements in table.
bool erase(const iterator &iter)
Erase an entry specified by given iterator.
A HashTable to objects of type <T> with a label key.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
label index() const noexcept
Return the hit index.
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
void size(const label n)
Older name for setAddressableSize.
static bool & parRun() noexcept
Test if this a parallel run.
const labelListList & normalDirections() const
const vectorField & edgeDirections() const
Return the edgeDirection vectors.
const List< sideVolumeType > & normalVolumeTypes() const
Return.
const labelListList & edgeNormals() const
Return the indices of the normals that are adjacent to the.
@ FLAT
Neither concave or convex, on a flat surface.
@ OPEN
Only connected to a single face.
@ MULTIPLE
Multiply connected (connected to more than two faces)
static const Enum< sideVolumeType > sideVolumeTypeNames_
const vectorField & normals() const
Return the normals of the surfaces adjacent to the feature edges.
sideVolumeType
Normals point to the outside.
@ NEITHER
not sure when this may be used
edgeStatus getEdgeStatus(label edgeI) const
Return the edgeStatus of a specified edge.
@ BOTH
limit by both minimum and maximum values
A reference point and direction.
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
@ FRONT
The front (positive normal) side of the plane.
@ BACK
The back (negative normal) side of the plane.
int myProcNo() const noexcept
Return processor number.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar pi(M_PI)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
static constexpr const zero Zero
Global zero (0)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.