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()