From 70116c004928d5dddb00f93f3ec2bb598ff0968e Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 20 Jan 2010 18:54:56 +0000 Subject: [PATCH] Handling cases where there are repeated eigenvalues for a face, i.e. a square. The subdeterminants in the eigenvector calculation in tensor.C ran off the bottom of SMALL. Handling by identifying the longest edge as the collapseAxis. Adding function to return the index of the longest edge for a face. --- .../conformalVoronoiMesh.H | 3 + .../conformalVoronoiMeshCalcDualMesh.C | 119 ++++++++++++------ 2 files changed, 84 insertions(+), 38 deletions(-) diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index ee6fcf10a0..619f9b4fb2 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -549,6 +549,9 @@ private: label maxFC ) const; + // Identify the index of the longest edge on the face + label longestEdge(const face& f, const pointField& pts) const; + //- Identify the face labels of the deferred collapse faces void deferredCollapseFaceSet ( diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C index f46b03d4c6..ea10b6ced5 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C @@ -865,6 +865,8 @@ Foam::conformalVoronoiMesh::collapseFace maxFC ); + labelList facePts(f); + const vector fC = f.centre(pts); vector fN = f.normal(pts); @@ -898,59 +900,50 @@ Foam::conformalVoronoiMesh::collapseFace vector collapseAxis = vector::zero; - scalar aspectRatio = 1; + scalar aspectRatio = 1.0; if (detJ < 1e-5) { - const edgeList& eds = f.edges(); - - label longestEdgeI = -1; - - scalar longestEdgeLength = -SMALL; - - forAll(eds, edI) - { - scalar edgeLength = eds[edI].mag(pts); - - if (edgeLength > longestEdgeLength) - { - longestEdgeI = edI; - - longestEdgeLength = edgeLength; - } - } - - collapseAxis = eds[longestEdgeI].vec(pts); - - if (mag(collapseAxis) < VSMALL) - { - Info<< "if (mag(collapseAxis) < VSMALL) " << collapseAxis << endl; - } + collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts); collapseAxis /= mag(collapseAxis); // Empirical correlation for high aspect ratio faces aspectRatio = sqrt(0.35/detJ); - - // Info<< "# Longest edge determined collapseAxis" << endl; } else { vector eVals = eigenValues(J); - // The maximum eigenvalue (z()) must be the direction of the - // normal, as it has the greatest value. The minimum eigenvalue - // is the dominant collapse axis for high aspect ratio faces. + if (mag(eVals.y() - eVals.x()) < 100*SMALL) + { + // First two eigenvalues are the same: i.e. a square face - collapseAxis = eigenVector(J, eVals.x()); + // Cannot necessarily determine linearly independent + // eigenvectors, or any at all, use longest edge direction. - // The inertia calculation describes the mass distribution as a - // function of distance squared to the axis, so the square root of - // the ratio of face-plane moments gives a good indication of the - // aspect ratio. + collapseAxis = f.edges()[longestEdge(f, pts)].vec(pts); - aspectRatio = sqrt(eVals.y()/max(eVals.x(), SMALL)); + collapseAxis /= mag(collapseAxis); + + aspectRatio = 1.0; + } + else + { + // The maximum eigenvalue (z()) must be the direction of the + // normal, as it has the greatest value. The minimum eigenvalue + // is the dominant collapse axis for high aspect ratio faces. + + collapseAxis = eigenVector(J, eVals.x()); + + // The inertia calculation describes the mass distribution as a + // function of distance squared to the axis, so the square root of + // the ratio of face-plane moments gives a good indication of the + // aspect ratio. + + aspectRatio = sqrt(eVals.y()/max(eVals.x(), SMALL)); + } } if (magSqr(collapseAxis) < VSMALL) @@ -962,6 +955,30 @@ Foam::conformalVoronoiMesh::collapseFace << "No collapse axis found for face, not collapsing." << endl; + // Output face and collapse axis for visualisation + + Info<< "# Aspect ratio = " << aspectRatio << nl + << "# inertia = " << J << nl + << "# determinant = " << detJ << nl + << "# eigenvalues = " << eigenValues(J) << nl + << "# collapseAxis = " << collapseAxis << nl + << "# facePts = " << facePts << nl + << endl; + + forAll(f, fPtI) + { + meshTools::writeOBJ(Info, pts[f[fPtI]]); + } + + Info<< "f"; + + forAll(f, fPtI) + { + Info << " " << fPtI + 1; + } + + Info<< endl; + return fcmNone; } @@ -980,8 +997,6 @@ Foam::conformalVoronoiMesh::collapseFace // Sort the projected distances and the corresponding vertex // indices along the collapse axis - labelList facePts(f); - labelList oldToNew; sortedOrder(d, oldToNew); @@ -1214,6 +1229,34 @@ Foam::conformalVoronoiMesh::collapseFace } +Foam::label Foam::conformalVoronoiMesh::longestEdge +( + const face& f, + const pointField& pts +) const +{ + const edgeList& eds = f.edges(); + + label longestEdgeI = -1; + + scalar longestEdgeLength = -SMALL; + + forAll(eds, edI) + { + scalar edgeLength = eds[edI].mag(pts); + + if (edgeLength > longestEdgeLength) + { + longestEdgeI = edI; + + longestEdgeLength = edgeLength; + } + } + + return longestEdgeI; +} + + void Foam::conformalVoronoiMesh::deferredCollapseFaceSet ( labelList& owner,