From 1372cf4c70e5598a7c11158271b47df0c0b40c37 Mon Sep 17 00:00:00 2001 From: mattijs <mattijs> Date: Thu, 21 Jul 2011 14:01:58 +0100 Subject: [PATCH] ENH: polyMeshTetDecomposition: restore face centre tet decomposition --- .../polyMeshTetDecomposition.C | 98 ++++++------- .../polyMeshGeometry/polyMeshGeometry.C | 137 +++++++++--------- 2 files changed, 117 insertions(+), 118 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C index d20a540a57f..b881a5ae774 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C @@ -380,11 +380,11 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets ) { const labelList& own = mesh.faceOwner(); - // const labelList& nei = mesh.faceNeighbour(); + const labelList& nei = mesh.faceNeighbour(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); const vectorField& cc = mesh.cellCentres(); - // const vectorField& fc = mesh.faceCentres(); + const vectorField& fc = mesh.faceCentres(); // Calculate coupled cell centre pointField neiCc(mesh.nFaces() - mesh.nInternalFaces()); @@ -398,62 +398,62 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets const faceList& fcs = mesh.faces(); - // const pointField& p = mesh.points(); + const pointField& p = mesh.points(); label nErrorTets = 0; forAll(fcs, faceI) { - // const face& f = fcs[faceI]; - - // forAll(f, fPtI) - // { - // scalar tetQual = tetPointRef - // ( - // p[f[fPtI]], - // p[f.nextLabel(fPtI)], - // fc[faceI], - // cc[own[faceI]] - // ).quality(); - - // if (tetQual > -tol) - // { - // if (setPtr) - // { - // setPtr->insert(faceI); - // } - - // nErrorTets++; - // break; // no need to check other tets - // } - // } + const face& f = fcs[faceI]; + + forAll(f, fPtI) + { + scalar tetQual = tetPointRef + ( + p[f[fPtI]], + p[f.nextLabel(fPtI)], + fc[faceI], + cc[own[faceI]] + ).quality(); + + if (tetQual > -tol) + { + if (setPtr) + { + setPtr->insert(faceI); + } + + nErrorTets++; + break; // no need to check other tets + } + } if (mesh.isInternalFace(faceI)) { // Create the neighbour tet - it will have positive volume - // const face& f = fcs[faceI]; - - // forAll(f, fPtI) - // { - // scalar tetQual = tetPointRef - // ( - // p[f[fPtI]], - // p[f.nextLabel(fPtI)], - // fc[faceI], - // cc[nei[faceI]] - // ).quality(); - - // if (tetQual < tol) - // { - // if (setPtr) - // { - // setPtr->insert(faceI); - // } - - // nErrorTets++; - // break; - // } - // } + const face& f = fcs[faceI]; + + forAll(f, fPtI) + { + scalar tetQual = tetPointRef + ( + p[f[fPtI]], + p[f.nextLabel(fPtI)], + fc[faceI], + cc[nei[faceI]] + ).quality(); + + if (tetQual < tol) + { + if (setPtr) + { + setPtr->insert(faceI); + } + + nErrorTets++; + break; + } + } if (findSharedBasePoint(mesh, faceI, tol, report) == -1) { diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index 8c84ac72f55..fa6d11d5009 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -309,6 +309,7 @@ Foam::scalar Foam::polyMeshGeometry::calcSkewness } +// Create the neighbour pyramid - it will have positive volume bool Foam::polyMeshGeometry::checkFaceTet ( const polyMesh& mesh, @@ -787,7 +788,7 @@ bool Foam::polyMeshGeometry::checkFaceTets // check whether decomposing each cell into tets results in // positive volume, non-flat tets const labelList& own = mesh.faceOwner(); - // const labelList& nei = mesh.faceNeighbour(); + const labelList& nei = mesh.faceNeighbour(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); // Calculate coupled cell centre @@ -802,50 +803,48 @@ bool Foam::polyMeshGeometry::checkFaceTets label nErrorTets = 0; - // bool tetError = false - forAll(checkFaces, i) { label faceI = checkFaces[i]; // Create the owner pyramid - note: exchange cell and face centre // to get positive volume. - // tetError = checkFaceTet - // ( - // mesh, - // report, - // minTetQuality, - // p, - // faceI, - // cellCentres[own[faceI]], // face centre - // faceCentres[faceI], // cell centre - // setPtr - // ); - - // if (tetError) - // { - // nErrorTets++; - // } + bool tetError = checkFaceTet + ( + mesh, + report, + minTetQuality, + p, + faceI, + cellCentres[own[faceI]], // face centre + faceCentres[faceI], // cell centre + setPtr + ); + + if (tetError) + { + nErrorTets++; + } if (mesh.isInternalFace(faceI)) { // Create the neighbour tets - they will have positive volume - // tetError = checkFaceTet - // ( - // mesh, - // report, - // minTetQuality, - // p, - // faceI, - // faceCentres[faceI], // face centre - // cellCentres[nei[faceI]], // cell centre - // setPtr - // ); - - // if (tetError) - // { - // nErrorTets++; - // } + bool tetError = checkFaceTet + ( + mesh, + report, + minTetQuality, + p, + faceI, + faceCentres[faceI], // face centre + cellCentres[nei[faceI]], // cell centre + setPtr + ); + + if (tetError) + { + nErrorTets++; + } if ( @@ -921,40 +920,40 @@ bool Foam::polyMeshGeometry::checkFaceTets label face0 = baffles[i].first(); label face1 = baffles[i].second(); - // tetError = checkFaceTet - // ( - // mesh, - // report, - // minTetQuality, - // p, - // face0, - // cellCentres[own[face0]], // face centre - // faceCentres[face0], // cell centre - // setPtr - // ); - - // if (tetError) - // { - // nErrorTets++; - // } - - // // Create the neighbour tets - they will have positive volume - // tetError = checkFaceTet - // ( - // mesh, - // report, - // minTetQuality, - // p, - // face0, - // faceCentres[face0], // face centre - // cellCentres[own[face1]], // cell centre - // setPtr - // ); - - // if (tetError) - // { - // nErrorTets++; - // } + bool tetError = checkFaceTet + ( + mesh, + report, + minTetQuality, + p, + face0, + cellCentres[own[face0]], // face centre + faceCentres[face0], // cell centre + setPtr + ); + + if (tetError) + { + nErrorTets++; + } + + // Create the neighbour tets - they will have positive volume + tetError = checkFaceTet + ( + mesh, + report, + minTetQuality, + p, + face0, + faceCentres[face0], // face centre + cellCentres[own[face1]], // cell centre + setPtr + ); + + if (tetError) + { + nErrorTets++; + } if ( -- GitLab