diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2D.C b/applications/utilities/mesh/generation/CV2DMesher/CV2D.C index 64c484bf7e..e9d0c1b284 100644 --- a/applications/utilities/mesh/generation/CV2DMesher/CV2D.C +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2D.C @@ -344,9 +344,59 @@ void Foam::CV2D::newPoints(const scalar relaxation) const vectorField& faceNormals = qSurf_.faceNormals(); - // Initialise the total displacement and its distance for writing out - vector2D totalDisp = vector2D::zero; - scalar totalDist = 0; + Field dualVertices(number_of_faces()); + + label dualVerti = 0; + + // Find the dual point of each tetrahedron and assign it an index. + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + fit->faceIndex() = -1; + + if + ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + { + fit->faceIndex() = dualVerti; + + dualVertices[dualVerti] = toPoint2D(circumcenter(fit)); + + dualVerti++; + } + } + + dualVertices.setSize(dualVerti); + + Field displacementAccumulator + ( + startOfSurfacePointPairs_, + vector2D::zero + ); + + // Calculate target size and alignment for vertices + scalarField sizes + ( + number_of_vertices(), + controls_.minCellSize + ); + + vector2D align(3,-6); + + align /= mag(align); + + Field alignments + ( + number_of_vertices(), + align + ); for ( @@ -355,186 +405,213 @@ void Foam::CV2D::newPoints(const scalar relaxation) ++vit ) { - if (vit->internalPoint()) + if (vit->internalOrBoundaryPoint()) { - // Current dual-cell defining vertex ("centre") - point2DFromPoint defVert0 = toPoint2D(vit->point()); + point2D vert = toPoint2D(vit->point()); - Triangulation::Edge_circulator ec = incident_edges(vit); - Triangulation::Edge_circulator ecStart = ec; + // alignment determination + pointIndexHit pHit = qSurf_.tree().findNearest + ( + toPoint3D(vert), + tols_.span2 + ); - // Circulate around the edges to find the first which is not - // infinite - do + if (pHit.hit()) { - if (!is_infinite(ec)) break; - } while (++ec != ecStart); - - // Store the start-end of the first non-infinte edge - point2D de0 = toPoint2D(circumcenter(ec->first)); - - // Keep track of the maximum edge length^2 - scalar maxEdgeLen2 = 0.0; - - // Keep track of the index of the longest edge - label edgecd0i = -1; - - // Edge counter - label edgei = 0; - - do - { - if (!is_infinite(ec)) - { - // Get the end of the current edge - point2D de1 = toPoint2D - ( - circumcenter(ec->first->neighbor(ec->second)) - ); - - // Store the current edge vector - edges[edgei] = de1 - de0; - - // Store the edge mid-point in the vertices array - vertices[edgei] = 0.5*(de1 + de0); - - // Move the current edge end into the edge start for the - // next iteration - de0 = de1; - - // Keep track of the longest edge - - scalar edgeLen2 = magSqr(edges[edgei]); - - if (edgeLen2 > maxEdgeLen2) - { - maxEdgeLen2 = edgeLen2; - edgecd0i = edgei; - } - - edgei++; - } - } while (++ec != ecStart); - - // Initialise cd0 such that the mesh will align - // in in the x-y directions - vector2D cd0(1, 0); - - if (controls_.relaxOrientation) - { - // Get the longest edge from the array and use as the primary - // direction of the coordinate system of the "square" cell - cd0 = edges[edgecd0i]; + alignments[vit->index()] = toPoint2D(faceNormals[pHit.index()]); } - if (controls_.nearWallAlignedDist > 0) + // size determination + + if (vert.x() > 0) { - pointIndexHit pHit = qSurf_.tree().findNearest + sizes[vit->index()] *= 0.5; + } + } + } + + // Info<< "Calculated alignments" << endl; + + scalar cosAlignmentAcceptanceAngle = 0.68; + + // Upper and lower edge length ratios for weight + scalar u = 1.0; + scalar l = 0.7; + + PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true); + + DynamicList pointsToInsert; + + for + ( + Triangulation::Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + eit++ + ) + { + Vertex_handle vA = eit->first->vertex(cw(eit->second)); + Vertex_handle vB = eit->first->vertex(ccw(eit->second)); + + if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint()) + { + continue; + } + + const point2D& dualV1 = dualVertices[eit->first->faceIndex()]; + const point2D& dualV2 = + dualVertices[eit->first->neighbor(eit->second)->faceIndex()]; + + scalar dualEdgeLength = mag(dualV1 - dualV2); + + point2D dVA = toPoint2D(vA->point()); + point2D dVB = toPoint2D(vB->point()); + + Field alignmentDirsA(2); + + alignmentDirsA[0] = alignments[vA->index()]; + alignmentDirsA[1] = vector2D + ( + -alignmentDirsA[0].y(), + alignmentDirsA[0].x() + ); + + Field alignmentDirsB(2); + + alignmentDirsB[0] = alignments[vB->index()]; + alignmentDirsB[1] = vector2D + ( + -alignmentDirsB[0].y(), + alignmentDirsB[0].x() + ); + + Field alignmentDirs(2); + + forAll(alignmentDirsA, aA) + { + const vector2D& a(alignmentDirsA[aA]); + + scalar maxDotProduct = 0.0; + + forAll(alignmentDirsB, aB) + { + const vector2D& b(alignmentDirsB[aB]); + + scalar dotProduct = a & b; + + if (mag(dotProduct) > maxDotProduct) + { + maxDotProduct = mag(dotProduct); + + alignmentDirs[aA] = a + sign(dotProduct)*b; + + alignmentDirs[aA] /= mag(alignmentDirs[aA]); + } + } + } + + vector2D rAB = dVA - dVB; + + scalar rABMag = mag(rAB); + + forAll(alignmentDirs, aD) + { + vector2D& alignmentDir = alignmentDirs[aD]; + + if ((rAB & alignmentDir) < 0) + { + // swap the direction of the alignment so that has the + // same sense as rAB + alignmentDir *= -1; + } + + scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir); + + if (alignmentDotProd > cosAlignmentAcceptanceAngle) + { + scalar targetFaceSize = + 0.5*(sizes[vA->index()] + sizes[vB->index()]); + + alignmentDir *= 0.5*targetFaceSize; + + vector2D delta = alignmentDir - 0.5*rAB; + + if (dualEdgeLength < 0.7*targetFaceSize) + { + delta *= 0; + } + else if (dualEdgeLength < targetFaceSize) + { + delta *= + ( + dualEdgeLength + /(targetFaceSize*(u - l)) + - 1/((u/l) - 1) + ); + } + + if ( - toPoint3D(defVert0), - controls_.nearWallAlignedDist2 - ); - - if (pHit.hit()) + vA->internalPoint() + && vB->internalPoint() + && rABMag > 1.75*targetFaceSize + && dualEdgeLength > 0.05*targetFaceSize + && alignmentDotProd > 0.93 + ) { - cd0 = toPoint2D(faceNormals[pHit.index()]); + // Point insertion + pointsToInsert.append(0.5*(dVA + dVB)); } - } - - // Rotate by 45deg needed to create an averaging procedure which - // encourages the cells to be square - cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); - - // Normalise the primary coordinate direction - cd0 /= mag(cd0); - - // Calculate the orthogonal coordinate direction - vector2D cd1(-cd0.y(), cd0.x()); - - - // Restart the circulator - ec = ecStart; - - // ... and the counter - edgei = 0; - - // Initialise the displacement for the centre and sum-weights - vector2D disp = vector2D::zero; - scalar sumw = 0; - - do - { - if (!is_infinite(ec)) + else if + ( + (vA->internalPoint() || vB->internalPoint()) + && rABMag < 0.65*targetFaceSize + ) { - // Pick up the current edge - const vector2D& ei = edges[edgei]; + // Point removal - // Calculate the centre to edge-centre vector - vector2D deltai = vertices[edgei] - defVert0; - - // Set the weight for this edge contribution - scalar w = 1; - - if (controls_.squares) + // Only insert a point at the midpoint of the short edge + // if neither attached point has already been identified + // to be removed. + if + ( + pointToBeRetained[vA->index()] == true + && pointToBeRetained[vB->index()] == true + ) { - w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x()); - // alternative weights - //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x()); - //w = magSqr(ei)*mag(deltai); - - // Use the following for an ~square mesh - // Find the coordinate contributions for this edge delta - scalar cd0deltai = cd0 & deltai; - scalar cd1deltai = cd1 & deltai; - - // Create a "square" displacement - if (mag(cd0deltai) > mag(cd1deltai)) - { - disp += (w*cd0deltai)*cd0; - } - else - { - disp += (w*cd1deltai)*cd1; - } - } - else - { - // Use this for a hexagon/pentagon mesh - disp += w*deltai; + pointsToInsert.append(0.5*(dVA + dVB)); } - // Sum the weights - sumw += w; + if (vA->internalPoint()) + { + pointToBeRetained[vA->index()] = false; + } + + if (vB->internalPoint()) + { + pointToBeRetained[vB->index()] = false; + } } else { - FatalErrorIn("CV2D::newPoints() const") - << "Infinite triangle found in internal mesh" - << exit(FatalError); + if (vA->internalPoint()) + { + displacementAccumulator[vA->index()] += delta; + } + + if (vB->internalPoint()) + { + displacementAccumulator[vB->index()] += -delta; + } } - - edgei++; - - } while (++ec != ecStart); - - // Calculate the average displacement - disp /= sumw; - totalDisp += disp; - totalDist += mag(disp); - - // Move the point by a fraction of the average displacement - movePoint(vit, defVert0 + relaxation*disp); + } } } - Info << "\nTotal displacement = " << totalDisp - << " total distance = " << totalDist << endl; -} + vector2D totalDisp = sum(displacementAccumulator); + scalar totalDist = sum(mag(displacementAccumulator)); - -void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) -{ - label pointI = 0; + // Relax the calculated displacement + displacementAccumulator *= relaxation; for ( @@ -545,12 +622,269 @@ void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) { if (vit->internalPoint()) { - movePoint(vit, newPoints[pointI++]); + if (!pointToBeRetained[vit->index()]) + { + remove(vit); + } + else + { + movePoint + ( + vit, + vit->point() + + K::Vector_2 + ( + displacementAccumulator[vit->index()].x(), + displacementAccumulator[vit->index()].y() + ) + ); + } } } + + removeSurfacePointPairs(); + + // Re-index internal points + + label nVert = startOfInternalPoints_; + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalPoint()) + { + vit->index() = nVert++; + } + } + + // Insert new points + + Info<< "Inserting " << pointsToInsert.size() << " new points" << endl; + + forAll(pointsToInsert, i) + { + insertPoint(pointsToInsert[i], Vb::INTERNAL_POINT); + } + + Info<< " Total displacement = " << totalDisp << nl + << " Total distance = " << totalDist << nl + << " Points added = " << pointsToInsert.size() + << endl; + + insertSurfacePointPairs(); + boundaryConform(); + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Old Method + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // for + // ( + // Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + // vit != finite_vertices_end(); + // ++vit + // ) + // { + // if (vit->internalPoint()) + // { + // // Current dual-cell defining vertex ("centre") + // point2DFromPoint defVert0 = toPoint2D(vit->point()); + + // Triangulation::Edge_circulator ec = incident_edges(vit); + // Triangulation::Edge_circulator ecStart = ec; + + // // Circulate around the edges to find the first which is not + // // infinite + // do + // { + // if (!is_infinite(ec)) break; + // } while (++ec != ecStart); + + // // Store the start-end of the first non-infinte edge + // point2D de0 = toPoint2D(circumcenter(ec->first)); + + // // Keep track of the maximum edge length^2 + // scalar maxEdgeLen2 = 0.0; + + // // Keep track of the index of the longest edge + // label edgecd0i = -1; + + // // Edge counter + // label edgei = 0; + + // do + // { + // if (!is_infinite(ec)) + // { + // // Get the end of the current edge + // point2D de1 = toPoint2D + // ( + // circumcenter(ec->first->neighbor(ec->second)) + // ); + + // // Store the current edge vector + // edges[edgei] = de1 - de0; + + // // Store the edge mid-point in the vertices array + // vertices[edgei] = 0.5*(de1 + de0); + + // // Move the current edge end into the edge start for the + // // next iteration + // de0 = de1; + + // // Keep track of the longest edge + + // scalar edgeLen2 = magSqr(edges[edgei]); + + // if (edgeLen2 > maxEdgeLen2) + // { + // maxEdgeLen2 = edgeLen2; + // edgecd0i = edgei; + // } + + // edgei++; + // } + // } while (++ec != ecStart); + + // // Initialise cd0 such that the mesh will align + // // in in the x-y directions + // vector2D cd0(1, 0); + + // if (controls_.relaxOrientation) + // { + // // Get the longest edge from the array and use as the primary + // // direction of the coordinate system of the "square" cell + // cd0 = edges[edgecd0i]; + // } + + // if (controls_.nearWallAlignedDist > 0) + // { + // pointIndexHit pHit = qSurf_.tree().findNearest + // ( + // toPoint3D(defVert0), + // controls_.nearWallAlignedDist2 + // ); + + // if (pHit.hit()) + // { + // cd0 = toPoint2D(faceNormals[pHit.index()]); + // } + // } + + // // Rotate by 45deg needed to create an averaging procedure which + // // encourages the cells to be square + // cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); + + // // Normalise the primary coordinate direction + // cd0 /= mag(cd0); + + // // Calculate the orthogonal coordinate direction + // vector2D cd1(-cd0.y(), cd0.x()); + + + // // Restart the circulator + // ec = ecStart; + + // // ... and the counter + // edgei = 0; + + // // Initialise the displacement for the centre and sum-weights + // vector2D disp = vector2D::zero; + // scalar sumw = 0; + + // do + // { + // if (!is_infinite(ec)) + // { + // // Pick up the current edge + // const vector2D& ei = edges[edgei]; + + // // Calculate the centre to edge-centre vector + // vector2D deltai = vertices[edgei] - defVert0; + + // // Set the weight for this edge contribution + // scalar w = 1; + + // if (controls_.squares) + // { + // w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x()); + // // alternative weights + // //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x()); + // //w = magSqr(ei)*mag(deltai); + + // // Use the following for an ~square mesh + // // Find the coordinate contributions for this edge delta + // scalar cd0deltai = cd0 & deltai; + // scalar cd1deltai = cd1 & deltai; + + // // Create a "square" displacement + // if (mag(cd0deltai) > mag(cd1deltai)) + // { + // disp += (w*cd0deltai)*cd0; + // } + // else + // { + // disp += (w*cd1deltai)*cd1; + // } + // } + // else + // { + // // Use this for a hexagon/pentagon mesh + // disp += w*deltai; + // } + + // // Sum the weights + // sumw += w; + // } + // else + // { + // FatalErrorIn("CV2D::newPoints() const") + // << "Infinite triangle found in internal mesh" + // << exit(FatalError); + // } + + // edgei++; + + // } while (++ec != ecStart); + + // // Calculate the average displacement + // disp /= sumw; + // totalDisp += disp; + // totalDist += mag(disp); + + // // Move the point by a fraction of the average displacement + // movePoint(vit, defVert0 + relaxation*disp); + // } + // } + + // Info << "\nTotal displacement = " << totalDisp + // << " total distance = " << totalDist << endl; } +// void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) +// { +// label pointI = 0; + +// for +// ( +// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); +// vit != finite_vertices_end(); +// ++vit +// ) +// { +// if (vit->internalPoint()) +// { +// movePoint(vit, newPoints[pointI++]); +// } +// } +// } + + void Foam::CV2D::write() const { if (controls_.writeFinalTriangulation) diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2D.H b/applications/utilities/mesh/generation/CV2DMesher/CV2D.H index d959816034..60900028ce 100644 --- a/applications/utilities/mesh/generation/CV2DMesher/CV2D.H +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2D.H @@ -29,7 +29,7 @@ Description Conformal-Voronoi 2D automatic mesher with grid or read initial points and point position relaxation with optional "squarification". - There are a substantial number of options to this mesher read from + There are a substantial number of options to this mesher read from CV2DMesherDict file e.g.: // Min cell size used in tolerances when inserting points for @@ -48,7 +48,7 @@ Description // Should the mesh be square-dominated or of unbiased hexagons squares yes; - // Near-wall region where cells are aligned with the wall specified as a + // Near-wall region where cells are aligned with the wall specified as a // number of cell layers nearWallAlignedDist 3; @@ -124,6 +124,7 @@ SourceFiles #include "point2DFieldFwd.H" #include "dictionary.H" #include "Switch.H" +#include "PackedBoolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -338,7 +339,7 @@ private: const Triangulation::Finite_vertices_iterator& vit ) const; - //- Insert point-pairs at the nearest points on the surface to the + //- Insert point-pairs at the nearest points on the surface to the // control vertex of dual-cells which intersect the boundary in order // to provide a boundary-layer mesh. // NB: This is not guaranteed to close the boundary @@ -374,7 +375,7 @@ private: void markNearBoundaryPoints(); - //- Restore the Delaunay contraint + //- Restore the Delaunay contraint void fast_restore_Delaunay(Vertex_handle vh); // Flip operations used by fast_restore_Delaunay @@ -434,7 +435,7 @@ public: const scalar nearness ); - //- Create the initial mesh from the internal points in the given + //- Create the initial mesh from the internal points in the given // file. Points outside the geometry are ignored. void insertPoints(const fileName& pointFileName); @@ -461,11 +462,11 @@ public: // Point motion - inline void movePoint(const Vertex_handle& vh, const point2D& p); + inline void movePoint(const Vertex_handle& vh, const Point& P); //- Move the internal points to the given new locations and update // the triangulation to ensure it is Delaunay - void moveInternalPoints(const point2DField& newPoints); + // void moveInternalPoints(const point2DField& newPoints); //- Calculate the displacements to create the new points void newPoints(const scalar relaxation); diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2DI.H b/applications/utilities/mesh/generation/CV2DMesher/CV2DI.H index 02cded78e8..eaa01fc26f 100644 --- a/applications/utilities/mesh/generation/CV2DMesher/CV2DI.H +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2DI.H @@ -85,7 +85,7 @@ inline void Foam::CV2D::insertPointPair surfPt - ppDistn, number_of_vertices() + 1 ); - + insertPoint(surfPt + ppDistn, master); } @@ -136,10 +136,16 @@ inline Foam::point Foam::CV2D::toPoint3D(const Point& P) const } -inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const point2D& p) +inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const Point& P) { - vh->set_point(toPoint(p)); - fast_restore_Delaunay(vh); + label i = vh->index(); + + move(vh, P); + + vh->index() = i; + + // vh->set_point(toPoint(p)); + // fast_restore_Delaunay(vh); } diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C b/applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C index c6562adb69..c1c14c90a0 100644 --- a/applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C @@ -99,10 +99,9 @@ int main(int argc, char *argv[]) ) *scalar(iter)/scalar(nIterations); + Info<< "Relaxation = " << relax << endl; + mesh.newPoints(relax); - mesh.removeSurfacePointPairs(); - mesh.insertSurfacePointPairs(); - mesh.boundaryConform(); } mesh.write(); diff --git a/applications/utilities/mesh/generation/CV2DMesher/Make/options b/applications/utilities/mesh/generation/CV2DMesher/Make/options index ff54c4ad91..ae6f0b006e 100755 --- a/applications/utilities/mesh/generation/CV2DMesher/Make/options +++ b/applications/utilities/mesh/generation/CV2DMesher/Make/options @@ -1,4 +1,5 @@ -//EXE_DEBUG = -DFULLDEBUG -g -O0 + +// EXE_DEBUG = -DFULLDEBUG -g -O0 EXE_FROUNDING_MATH = -frounding-math EXE_NDEBUG = -DNDEBUG