42 void Foam::CV2D::insertBoundingBox()
44 Info<<
"insertBoundingBox: creating bounding mesh" <<
endl;
53 void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh)
56 Face_handle
f = vh->face(), next, start(
f);
63 if (!internal_flip(
f, cw(i))) external_flip(
f, i);
64 if (
f->neighbor(i) == start) start =
f;
66 f =
f->neighbor(cw(i));
71 void Foam::CV2D::external_flip(Face_handle&
f,
int i)
73 Face_handle
n =
f->neighbor(i);
77 CGAL::ON_POSITIVE_SIDE
78 != side_of_oriented_circle(
n,
f->vertex(i)->point())
82 i =
n->index(
f->vertex(i));
87 bool Foam::CV2D::internal_flip(Face_handle&
f,
int i)
89 Face_handle
n =
f->neighbor(i);
93 CGAL::ON_POSITIVE_SIDE
94 != side_of_oriented_circle(
n,
f->vertex(i)->point())
111 const dictionary& cvMeshDict
116 rndGen_(64293*Pstream::myProcNo()),
121 "cvSearchableSurfaces",
128 cvMeshDict.subDict(
"geometry"),
129 cvMeshDict.getOrDefault(
"singleRegionName", true)
136 cvMeshDict.subDict(
"surfaceConformation")
138 controls_(cvMeshDict, qSurf_.globalBounds()),
142 cvMeshDict.subDict(
"motionControl").subDict(
"shapeControlFunctions"),
144 controls_.minCellSize()
150 cvMeshDict.subDict(
"motionControl"),
156 cvMeshDict.subDict(
"surfaceConformation")
159 startOfInternalPoints_(0),
160 startOfSurfacePointPairs_(0),
161 startOfBoundaryConformPointPairs_(0),
167 insertFeaturePoints();
182 const scalar nearness
185 Info<<
"insertInitialPoints(const point2DField& points): ";
187 startOfInternalPoints_ = number_of_vertices();
188 label nVert = startOfInternalPoints_;
193 if (qSurf_.wellInside(toPoint3D(
p), nearness))
200 <<
"Rejecting point " <<
p <<
" outside surface" <<
endl;
204 Info<< nVert <<
" vertices inserted" <<
endl;
206 if (meshControls().objOutput())
211 writeTriangles(
"initial_triangles.obj",
true);
212 writeFaces(
"initial_faces.obj",
true);
219 IFstream pointsFile(pointFileName);
221 if (pointsFile.good())
226 0.5*meshControls().minCellSize2()
232 <<
"Could not open pointsFile " << pointFileName
240 Info<<
"insertInitialGrid: ";
242 startOfInternalPoints_ = number_of_vertices();
243 label nVert = startOfInternalPoints_;
245 scalar x0 = qSurf_.globalBounds().min().x();
246 scalar xR = qSurf_.globalBounds().max().x() - x0;
247 int ni = int(xR/meshControls().minCellSize()) + 1;
248 scalar deltax = xR/ni;
250 scalar
y0 = qSurf_.globalBounds().min().y();
251 scalar yR = qSurf_.globalBounds().max().y() -
y0;
252 int nj = int(yR/meshControls().minCellSize()) + 1;
253 scalar deltay = yR/nj;
256 scalar pert = meshControls().randomPerturbation()*
min(deltax, deltay);
258 for (
int i=0; i<ni; i++)
260 for (
int j=0; j<nj; j++)
264 if (meshControls().randomiseInitialGrid())
270 if (qSurf_.wellInside(
p, 0.5*meshControls().minCellSize2()))
277 Info<< nVert <<
" vertices inserted" <<
endl;
279 if (meshControls().objOutput())
284 writeTriangles(
"initial_triangles.obj",
true);
285 writeFaces(
"initial_faces.obj",
true);
292 startOfSurfacePointPairs_ = number_of_vertices();
294 if (meshControls().insertSurfaceNearestPointPairs())
296 insertSurfaceNearestPointPairs();
303 if (meshControls().insertSurfaceNearPointPairs())
305 insertSurfaceNearPointPairs();
308 startOfBoundaryConformPointPairs_ = number_of_vertices();
314 if (!meshControls().insertSurfaceNearestPointPairs())
316 markNearBoundaryPoints();
322 Triangulation::Finite_faces_iterator fit = finite_faces_begin();
323 fit != finite_faces_end();
330 for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++)
332 label nIntersections = insertBoundaryConformPointPairs
334 "surfaceIntersections_" +
Foam::name(iter) +
".obj" 337 if (nIntersections == 0)
343 Info<<
"BC iteration " << iter <<
": " 344 << nIntersections <<
" point-pairs inserted" <<
endl;
352 Triangulation::Finite_faces_iterator fit = finite_faces_begin();
353 fit != finite_faces_end();
378 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
379 vit != finite_vertices_end();
383 if (vit->index() >= startOfSurfacePointPairs_)
393 const scalar relaxation = relaxationModel_->relaxation();
395 Info<<
"Relaxation = " << relaxation <<
endl;
397 Field<point2D> dualVertices(number_of_faces());
404 Triangulation::Finite_faces_iterator fit = finite_faces_begin();
405 fit != finite_faces_end();
409 fit->faceIndex() = -1;
413 fit->vertex(0)->internalOrBoundaryPoint()
414 || fit->vertex(1)->internalOrBoundaryPoint()
415 || fit->vertex(2)->internalOrBoundaryPoint()
418 fit->faceIndex() = dualVerti;
420 dualVertices[dualVerti] = toPoint2D(circumcenter(fit));
426 dualVertices.setSize(dualVerti);
428 Field<vector2D> displacementAccumulator
430 startOfSurfacePointPairs_,
437 number_of_vertices(),
438 meshControls().minCellSize()
441 Field<vector2D> alignments
443 number_of_vertices(),
449 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
450 vit != finite_vertices_end();
454 if (vit->internalOrBoundaryPoint())
456 point2D vert = toPoint2D(vit->point());
460 label hitSurface = -1;
462 qSurf_.findSurfaceNearest
465 meshControls().span2(),
473 allGeometry_[hitSurface].getNormal
475 List<pointIndexHit>(1, pHit),
479 alignments[vit->index()] = toPoint2D(norm[0]);
481 sizes[vit->index()] =
482 cellSizeControl_.cellSize
484 toPoint3D(vit->point())
492 scalar cosAlignmentAcceptanceAngle = 0.68;
498 bitSet pointToBeRetained(startOfSurfacePointPairs_,
true);
500 std::list<Point> pointsToInsert;
504 Triangulation::Finite_edges_iterator eit = finite_edges_begin();
505 eit != finite_edges_end();
509 Vertex_handle vA = eit->first->vertex(cw(eit->second));
510 Vertex_handle vB = eit->first->vertex(ccw(eit->second));
512 if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint())
517 const point2D& dualV1 = dualVertices[eit->first->faceIndex()];
519 dualVertices[eit->first->neighbor(eit->second)->faceIndex()];
521 scalar dualEdgeLength =
mag(dualV1 - dualV2);
523 point2D dVA = toPoint2D(vA->point());
524 point2D dVB = toPoint2D(vB->point());
526 Field<vector2D> alignmentDirsA(2);
528 alignmentDirsA[0] = alignments[vA->index()];
531 -alignmentDirsA[0].
y(),
532 alignmentDirsA[0].
x()
535 Field<vector2D> alignmentDirsB(2);
537 alignmentDirsB[0] = alignments[vB->index()];
540 -alignmentDirsB[0].
y(),
541 alignmentDirsB[0].
x()
544 Field<vector2D> alignmentDirs(alignmentDirsA);
546 forAll(alignmentDirsA, aA)
548 const vector2D& a(alignmentDirsA[aA]);
550 scalar maxDotProduct = 0.0;
552 forAll(alignmentDirsB, aB)
556 scalar dotProduct = a &
b;
558 if (
mag(dotProduct) > maxDotProduct)
560 maxDotProduct =
mag(dotProduct);
562 alignmentDirs[aA] = a +
sign(dotProduct)*
b;
564 alignmentDirs[aA].normalise();
571 scalar rABMag =
mag(rAB);
575 vector2D& alignmentDir = alignmentDirs[aD];
577 if ((rAB & alignmentDir) < 0)
584 scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
586 if (alignmentDotProd > cosAlignmentAcceptanceAngle)
588 scalar targetFaceSize =
589 0.5*(sizes[vA->index()] + sizes[vB->index()]);
598 alignmentDir *= 0.5*targetFaceSize;
602 if (dualEdgeLength < 0.7*targetFaceSize)
606 else if (dualEdgeLength < targetFaceSize)
611 /(targetFaceSize*(u - l))
619 && vB->internalPoint()
620 && rABMag > 1.75*targetFaceSize
621 && dualEdgeLength > 0.05*targetFaceSize
622 && alignmentDotProd > 0.93
626 pointsToInsert.push_back(
toPoint(0.5*(dVA + dVB)));
630 (vA->internalPoint() || vB->internalPoint())
631 && rABMag < 0.65*targetFaceSize
641 pointToBeRetained.test(vA->index())
642 && pointToBeRetained.test(vB->index())
645 pointsToInsert.push_back(
toPoint(0.5*(dVA + dVB)));
648 if (vA->internalPoint())
650 pointToBeRetained.unset(vA->index());
653 if (vB->internalPoint())
655 pointToBeRetained.unset(vB->index());
660 if (vA->internalPoint())
662 displacementAccumulator[vA->index()] +=
delta;
665 if (vB->internalPoint())
667 displacementAccumulator[vB->index()] += -
delta;
675 scalar totalDist =
sum(
mag(displacementAccumulator));
678 displacementAccumulator *= relaxation;
680 label numberOfNewPoints = pointsToInsert.size();
684 Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
685 vit != finite_vertices_end();
689 if (vit->internalPoint())
691 if (pointToBeRetained.test(vit->index()))
693 pointsToInsert.push_front
697 toPoint2D(vit->point())
698 + displacementAccumulator[vit->index()]
711 reinsertFeaturePoints();
713 startOfInternalPoints_ = number_of_vertices();
715 label nVert = startOfInternalPoints_;
717 Info<<
"Inserting " << numberOfNewPoints <<
" new points" <<
endl;
720 insert(pointsToInsert.begin(), pointsToInsert.end());
724 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
725 vit != finite_vertices_end();
735 vit->index() = nVert++;
739 Info<<
" Total displacement = " << totalDisp <<
nl 740 <<
" Total distance = " << totalDist <<
nl 741 <<
" Points added = " << pointsToInsert.size()
746 insertSurfacePointPairs();
959 if (meshControls().objOutput())
961 writeFaces(
"allFaces.obj",
false);
962 writeFaces(
"faces.obj",
true);
963 writeTriangles(
"allTriangles.obj",
false);
964 writeTriangles(
"triangles.obj",
true);
965 writePatch(
"patch.pch");
972 if (meshControls().objOutput())
981 + runTime_.timeName()
989 +
"Triangles/allTriangles_" 990 + runTime_.timeName()
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
dimensionedScalar sign(const dimensionedScalar &ds)
CGAL::Delaunay_triangulation_3< K, Tds, CompactLocator > Delaunay
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
vector2D point2D
Point2D is a vector.
errorManipArg< error, int > exit(error &err, const int errNo=1)
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...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void insertSurfacePointPairs()
Insert all surface point-pairs from.
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
constexpr char nl
The newline '\n' character (0x0a)
scalar span() const
Return the span.
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
const cv2DControls & meshControls() const
void newPoints()
Move the internal points to the given new locations and update.
dimensionedScalar y0(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
void insertGrid()
Create the initial mesh as a regular grid of points.
Type sample01()
Return a sample whose components lie in the range [0,1].
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
PointFrompoint toPoint(const Foam::point &p)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void insertPoints(const point2DField &points, const scalar nearness)
Create the initial mesh from the given internal points.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
defineTypeNameAndDebug(combustionModel, 0)
messageStream Warning
Warning stream (stdout output on master, null elsewhere), with additional 'FOAM Warning' header text...
vector2DField point2DField
point2DField is a vector2DField.
vector point
Point is a vector.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
void removeSurfacePointPairs()
Remove the point-pairs introduced by insertSurfacePointPairs.
static const Vector2D< scalar > zero
void boundaryConform()
Insert point-pairs where there are protrusions into.