From 0a26787282d6c98f6bfc1745f5c2b044f7b4ef3a Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 6 Oct 2010 10:22:12 +0100 Subject: [PATCH 01/12] ENH: Moving pointIndexHit to OpenFOAM, along-side pointHit. --- src/OpenFOAM/Make/files | 1 + .../meshes/primitiveShapes/objectHit}/PointIndexHit.H | 0 .../meshes/primitiveShapes/objectHit}/pointHitSort.H | 0 .../meshes/primitiveShapes/objectHit}/pointIndexHit.H | 0 .../meshes/primitiveShapes/objectHit}/pointIndexHitIOList.C | 0 .../meshes/primitiveShapes/objectHit}/pointIndexHitIOList.H | 0 src/meshTools/Make/files | 1 - 7 files changed, 1 insertion(+), 1 deletion(-) rename src/{meshTools/octree => OpenFOAM/meshes/primitiveShapes/objectHit}/PointIndexHit.H (100%) rename src/{meshTools/octree => OpenFOAM/meshes/primitiveShapes/objectHit}/pointHitSort.H (100%) rename src/{meshTools/octree => OpenFOAM/meshes/primitiveShapes/objectHit}/pointIndexHit.H (100%) rename src/{meshTools/octree => OpenFOAM/meshes/primitiveShapes/objectHit}/pointIndexHitIOList.C (100%) rename src/{meshTools/octree => OpenFOAM/meshes/primitiveShapes/objectHit}/pointIndexHitIOList.H (100%) diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index d9548b67b1..94bc86f2ef 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -296,6 +296,7 @@ primitiveShapes = meshes/primitiveShapes $(primitiveShapes)/line/line.C $(primitiveShapes)/plane/plane.C $(primitiveShapes)/triangle/intersection.C +$(primitiveShapes)/objectHit/pointIndexHitIOList.C meshShapes = meshes/meshShapes $(meshShapes)/edge/edge.C diff --git a/src/meshTools/octree/PointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H similarity index 100% rename from src/meshTools/octree/PointIndexHit.H rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H diff --git a/src/meshTools/octree/pointHitSort.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H similarity index 100% rename from src/meshTools/octree/pointHitSort.H rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H diff --git a/src/meshTools/octree/pointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H similarity index 100% rename from src/meshTools/octree/pointIndexHit.H rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H diff --git a/src/meshTools/octree/pointIndexHitIOList.C b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C similarity index 100% rename from src/meshTools/octree/pointIndexHitIOList.C rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C diff --git a/src/meshTools/octree/pointIndexHitIOList.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H similarity index 100% rename from src/meshTools/octree/pointIndexHitIOList.H rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 14e561fcba..2c9f81195d 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -44,7 +44,6 @@ octree/octreeDataFace.C octree/treeBoundBox.C octree/treeNodeName.C octree/treeLeafName.C -octree/pointIndexHitIOList.C indexedOctree/indexedOctreeName.C indexedOctree/treeDataCell.C From 6cdbd0ada7527bd3aec97bd657deec926eec85cd Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 6 Oct 2010 10:25:26 +0100 Subject: [PATCH 02/12] ENH: Making nearestPointClassify query for triangle. This is to access the face/edge/point status of the nearest at the same time to ensure a consistent result. Using getVolumeType query in distanceSurface, not simple normal dot-product comparison, fails on edges. --- .../primitiveShapes/triangle/triangle.H | 44 +- .../primitiveShapes/triangle/triangleI.H | 426 +++++------------- .../indexedOctree/treeDataTriSurface.C | 320 ++++++------- .../indexedOctree/treeDataTriSurface.H | 28 +- .../surfaceIntersection/edgeIntersections.C | 5 +- .../surfaceIntersection/surfaceIntersection.C | 2 +- .../octreeData/octreeDataTriSurface.C | 2 +- .../orientedSurface/orientedSurface.C | 4 +- .../triSurfaceTools/triSurfaceTools.C | 113 +++-- .../triSurfaceTools/triSurfaceTools.H | 4 +- .../distanceSurface/distanceSurface.C | 72 ++- 11 files changed, 439 insertions(+), 581 deletions(-) diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 64d9318306..abb5d7d8e4 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -79,21 +79,6 @@ class triangle PointRef a_, b_, c_; - // Private Member Functions - - //- Fast distance to triangle calculation. From - // "Distance Between Point and Trangle in 3D" - // David Eberly, Magic Software Inc. Aug. 2002. - // Works on function Q giving distance to point and tries to - // minimize this. - static pointHit nearestPoint - ( - const Point& baseVertex, - const vector& E0, - const vector& E1, - const point& P - ); - public: @@ -202,24 +187,27 @@ public: const scalar tol = 0.0 ) const; - //- Return nearest point to p on triangle - inline pointHit nearestPoint - ( - const point& p - ) const; - - //- Classify point in triangle plane w.r.t. triangle edges. - // - inside (true returned)/outside (false returned) - // - near point (nearType=POINT, nearLabel=0, 1, 2) - // - near edge (nearType=EDGE, nearLabel=0, 1, 2) + //- Find the nearest point to p on the triangle and classify it: + // + near point (nearType=POINT, nearLabel=0, 1, 2) + // + near edge (nearType=EDGE, nearLabel=0, 1, 2) // Note: edges are counted from starting // vertex so e.g. edge 2 is from f[2] to f[0] - // tol is fraction to account for truncation error. Is only used - // when comparing normalized (0..1) numbers. + pointHit nearestPointClassify + ( + const point& p, + label& nearType, + label& nearLabel + ) const; + + //- Return nearest point to p on triangle + inline pointHit nearestPoint(const point& p) const; + + //- Classify nearest point to p in triangle plane + // w.r.t. triangle edges and points. Returns inside + // (true)/outside (false). bool classify ( const point& p, - const scalar tol, label& nearType, label& nearLabel ) const; diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 10267a923f..7840d8ee7b 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -32,158 +32,6 @@ License namespace Foam { -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template -pointHit triangle::nearestPoint -( - const Point& baseVertex, - const vector& E0, - const vector& E1, - const point& P -) -{ - // Distance vector - const vector D(baseVertex - P); - - // Some geometrical factors - const scalar a = E0 & E0; - const scalar b = E0 & E1; - const scalar c = E1 & E1; - - // Precalculate distance factors - const scalar d = E0 & D; - const scalar e = E1 & D; - const scalar f = D & D; - - // Do classification - const scalar det = a*c - b*b; - scalar s = b*e - c*d; - scalar t = b*d - a*e; - - bool inside = false; - - if (s+t < det) - { - if (s < 0) - { - if (t < 0) - { - // Region 4 - if (e > 0) - { - // min on edge t = 0 - t = 0; - s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); - } - else - { - // min on edge s=0 - s = 0; - t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); - } - } - else - { - // Region 3. Min on edge s = 0 - s = 0; - t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); - } - } - else if (t < 0) - { - // Region 5 - t = 0; - s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); - } - else - { - // Region 0 - const scalar invDet = 1/det; - s *= invDet; - t *= invDet; - - inside = true; - } - } - else - { - if (s < 0) - { - // Region 2 - const scalar tmp0 = b + d; - const scalar tmp1 = c + e; - if (tmp1 > tmp0) - { - // min on edge s+t=1 - const scalar numer = tmp1 - tmp0; - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - t = 1 - s; - } - else - { - // min on edge s=0 - s = 0; - t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c)); - } - } - else if (t < 0) - { - // Region 6 - const scalar tmp0 = b + d; - const scalar tmp1 = c + e; - if (tmp1 > tmp0) - { - // min on edge s+t=1 - const scalar numer = tmp1 - tmp0; - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - t = 1 - s; - } - else - { - // min on edge t=0 - t = 0; - s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a)); - } - } - else - { - // Region 1 - const scalar numer = c+e-(b+d); - if (numer <= 0) - { - s = 0; - } - else - { - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - } - } - - t = 1 - s; - } - - // Calculate distance. - // Note: Foam::mag used since truncation error causes negative distances - // with points very close to one of the triangle vertices. - // (Up to -2.77556e-14 seen). Could use +SMALL but that not large enough. - - return pointHit - ( - inside, - baseVertex + s*E0 + t*E1, - Foam::sqrt - ( - Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f) - ), - !inside - ); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template @@ -247,7 +95,7 @@ inline Point triangle::centre() const template inline scalar triangle::mag() const { - return ::Foam::mag(normal()); + return Foam::mag(normal()); } @@ -536,7 +384,7 @@ inline pointHit triangle::ray inter.setMiss(eligible); // The miss point is the nearest point on the triangle - inter.setPoint(nearestPoint(a_, E0, E1, p).rawPoint()); + inter.setPoint(nearestPoint(p).rawPoint()); // The distance to the miss is the distance between the // original point and plane of intersection @@ -633,18 +481,130 @@ inline pointHit triangle::intersection } +template +pointHit triangle::nearestPointClassify +( + const point& p, + label& nearType, + label& nearLabel +) const +{ + // Adapted from: + // Real-time collision detection, Christer Ericson, 2005, 136-142 + + // Check if P in vertex region outside A + vector ab = b_ - a_; + vector ac = c_ - a_; + vector ap = p - a_; + + scalar d1 = ab & ap; + scalar d2 = ac & ap; + + if (d1 <= 0.0 && d2 <= 0.0) + { + // barycentric coordinates (1, 0, 0) + + nearType = POINT; + nearLabel = 0; + return pointHit(false, a_, Foam::mag(a_ - p), true); + } + + // Check if P in vertex region outside B + vector bp = p - b_; + scalar d3 = ab & bp; + scalar d4 = ac & bp; + + if (d3 >= 0.0 && d4 <= d3) + { + // barycentric coordinates (0, 1, 0) + + nearType = POINT; + nearLabel = 1; + return pointHit(false, b_, Foam::mag(b_ - p), true); + } + + // Check if P in edge region of AB, if so return projection of P onto AB + scalar vc = d1*d4 - d3*d2; + + if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) + { + // barycentric coordinates (1-v, v, 0) + scalar v = d1/(d1 - d3); + + point nearPt = a_ + v*ab; + nearType = EDGE; + nearLabel = 0; + return pointHit(false, nearPt, Foam::mag(nearPt - p), true); + } + + // Check if P in vertex region outside C + vector cp = p - c_; + scalar d5 = ab & cp; + scalar d6 = ac & cp; + + if (d6 >= 0.0 && d5 <= d6) + { + // barycentric coordinates (0, 0, 1) + + nearType = POINT; + nearLabel = 2; + return pointHit(false, c_, Foam::mag(c_ - p), true); + } + + // Check if P in edge region of AC, if so return projection of P onto AC + scalar vb = d5*d2 - d1*d6; + + if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) + { + // barycentric coordinates (1-w, 0, w) + scalar w = d2/(d2 - d6); + + point nearPt = a_ + w*ac; + nearType = EDGE; + nearLabel = 2; + return pointHit(false, nearPt, Foam::mag(nearPt - p), true); + } + + // Check if P in edge region of BC, if so return projection of P onto BC + scalar va = d3*d6 - d5*d4; + + if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) + { + // barycentric coordinates (0, 1-w, w) + scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6)); + + point nearPt = b_ + w*(c_ - b_); + nearType = EDGE; + nearLabel = 1; + return pointHit(false, nearPt, Foam::mag(nearPt - p), true); + } + + // P inside face region. Compute Q through its barycentric + // coordinates (u, v, w) + scalar denom = 1.0/(va + vb + vc); + scalar v = vb * denom; + scalar w = vc * denom; + + // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w + + point nearPt = a_ + ab*v + ac*w; + nearType = NONE, + nearLabel = -1; + return pointHit(true, nearPt, Foam::mag(nearPt - p), false); +} + + template inline pointHit triangle::nearestPoint ( const point& p ) const { - // Express triangle in terms of baseVertex (point a_) and - // two edge vectors - vector E0 = b_ - a_; - vector E1 = c_ - a_; + // Dummy labels + label nearType = -1; + label nearLabel = -1; - return nearestPoint(a_, E0, E1, p); + return nearestPointClassify(p, nearType, nearLabel); } @@ -652,160 +612,14 @@ template inline bool triangle::classify ( const point& p, - const scalar tol, label& nearType, label& nearLabel ) const { - const vector E0 = b_ - a_; - const vector E1 = c_ - a_; - const vector n = 0.5*(E0 ^ E1); - - // Get largest component of normal - scalar magX = Foam::mag(n.x()); - scalar magY = Foam::mag(n.y()); - scalar magZ = Foam::mag(n.z()); - - label i0 = -1; - if ((magX >= magY) && (magX >= magZ)) - { - i0 = 0; - } - else if ((magY >= magX) && (magY >= magZ)) - { - i0 = 1; - } - else - { - i0 = 2; - } - - // Get other components - label i1 = (i0 + 1) % 3; - label i2 = (i1 + 1) % 3; - - - scalar u1 = E0[i1]; - scalar v1 = E0[i2]; - - scalar u2 = E1[i1]; - scalar v2 = E1[i2]; - - scalar det = v2*u1 - u2*v1; - - scalar u0 = p[i1] - a_[i1]; - scalar v0 = p[i2] - a_[i2]; - - scalar alpha = 0; - scalar beta = 0; - - bool hit = false; - - if (Foam::mag(u1) < ROOTVSMALL) - { - beta = u0/u2; - - alpha = (v0 - beta*v2)/v1; - - hit = ((alpha >= 0) && ((alpha + beta) <= 1)); - } - else - { - beta = (v0*u1 - u0*v1)/det; - - alpha = (u0 - beta*u2)/u1; - - hit = ((alpha >= 0) && ((alpha + beta) <= 1)); - } - - // - // Now alpha, beta are the coordinates in the triangle local coordinate - // system E0, E1 - // - - //Pout<< "alpha:" << alpha << endl; - //Pout<< "beta:" << beta << endl; - //Pout<< "hit:" << hit << endl; - //Pout<< "tol:" << tol << endl; - - if (hit) - { - // alpha,beta might get negative due to precision errors - alpha = max(0.0, min(1.0, alpha)); - beta = max(0.0, min(1.0, beta)); - } - - nearType = NONE; - nearLabel = -1; - - if (Foam::mag(alpha+beta-1) <= tol) - { - // On edge between vert 1-2 (=edge 1) - - if (Foam::mag(alpha) <= tol) - { - nearType = POINT; - nearLabel = 2; - } - else if (Foam::mag(beta) <= tol) - { - nearType = POINT; - nearLabel = 1; - } - else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1)) - { - nearType = EDGE; - nearLabel = 1; - } - } - else if (Foam::mag(alpha) <= tol) - { - // On edge between vert 2-0 (=edge 2) - - if (Foam::mag(beta) <= tol) - { - nearType = POINT; - nearLabel = 0; - } - else if (Foam::mag(beta-1) <= tol) - { - nearType = POINT; - nearLabel = 2; - } - else if ((beta >= 0) && (beta <= 1)) - { - nearType = EDGE; - nearLabel = 2; - } - } - else if (Foam::mag(beta) <= tol) - { - // On edge between vert 0-1 (= edge 0) - - if (Foam::mag(alpha) <= tol) - { - nearType = POINT; - nearLabel = 0; - } - else if (Foam::mag(alpha-1) <= tol) - { - nearType = POINT; - nearLabel = 1; - } - else if ((alpha >= 0) && (alpha <= 1)) - { - nearType = EDGE; - nearLabel = 0; - } - } - - return hit; + return nearestPointClassify(p, nearType, nearLabel).hit(); } - - - // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // template diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C index e0f6d1a9df..b019a85304 100644 --- a/src/meshTools/indexedOctree/treeDataTriSurface.C +++ b/src/meshTools/indexedOctree/treeDataTriSurface.C @@ -35,145 +35,145 @@ defineTypeNameAndDebug(Foam::treeDataTriSurface, 0); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -// Fast distance to triangle calculation. From -// "Distance Between Point and Trangle in 3D" -// David Eberly, Magic Software Inc. Aug. 2003. -// Works on function Q giving distance to point and tries to minimize this. -Foam::scalar Foam::treeDataTriSurface::nearestCoords -( - const point& base, - const point& E0, - const point& E1, - const scalar a, - const scalar b, - const scalar c, - const point& P, - scalar& s, - scalar& t -) -{ - // distance vector - const vector D(base - P); +// // Fast distance to triangle calculation. From +// // "Distance Between Point and Triangle in 3D" +// // David Eberly, Magic Software Inc. Aug. 2003. +// // Works on function Q giving distance to point and tries to minimize this. +// Foam::scalar Foam::treeDataTriSurface::nearestCoords +// ( +// const point& base, +// const point& E0, +// const point& E1, +// const scalar a, +// const scalar b, +// const scalar c, +// const point& P, +// scalar& s, +// scalar& t +// ) +// { +// // distance vector +// const vector D(base - P); - // Precalculate distance factors. - const scalar d = E0 & D; - const scalar e = E1 & D; +// // Precalculate distance factors. +// const scalar d = E0 & D; +// const scalar e = E1 & D; - // Do classification - const scalar det = a*c - b*b; +// // Do classification +// const scalar det = a*c - b*b; - s = b*e - c*d; - t = b*d - a*e; +// s = b*e - c*d; +// t = b*d - a*e; - if (s+t < det) - { - if (s < 0) - { - if (t < 0) - { - //region 4 - if (e > 0) - { - //min on edge t = 0 - t = 0; - s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); - } - else - { - //min on edge s=0 - s = 0; - t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); - } - } - else - { - //region 3. Min on edge s = 0 - s = 0; - t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); - } - } - else if (t < 0) - { - //region 5 - t = 0; - s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); - } - else - { - //region 0 - const scalar invDet = 1/det; - s *= invDet; - t *= invDet; - } - } - else - { - if (s < 0) - { - //region 2 - const scalar tmp0 = b + d; - const scalar tmp1 = c + e; - if (tmp1 > tmp0) - { - //min on edge s+t=1 - const scalar numer = tmp1 - tmp0; - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - t = 1 - s; - } - else - { - //min on edge s=0 - s = 0; - t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c)); - } - } - else if (t < 0) - { - //region 6 - const scalar tmp0 = b + d; - const scalar tmp1 = c + e; - if (tmp1 > tmp0) - { - //min on edge s+t=1 - const scalar numer = tmp1 - tmp0; - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - t = 1 - s; - } - else - { - //min on edge t=0 - t = 0; - s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a)); - } - } - else - { - //region 1 - const scalar numer = c+e-(b+d); - if (numer <= 0) - { - s = 0; - } - else - { - const scalar denom = a-2*b+c; - s = (numer >= denom ? 1 : numer/denom); - } - } - t = 1 - s; - } +// if (s+t < det) +// { +// if (s < 0) +// { +// if (t < 0) +// { +// //region 4 +// if (e > 0) +// { +// //min on edge t = 0 +// t = 0; +// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); +// } +// else +// { +// //min on edge s=0 +// s = 0; +// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); +// } +// } +// else +// { +// //region 3. Min on edge s = 0 +// s = 0; +// t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c)); +// } +// } +// else if (t < 0) +// { +// //region 5 +// t = 0; +// s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a)); +// } +// else +// { +// //region 0 +// const scalar invDet = 1/det; +// s *= invDet; +// t *= invDet; +// } +// } +// else +// { +// if (s < 0) +// { +// //region 2 +// const scalar tmp0 = b + d; +// const scalar tmp1 = c + e; +// if (tmp1 > tmp0) +// { +// //min on edge s+t=1 +// const scalar numer = tmp1 - tmp0; +// const scalar denom = a-2*b+c; +// s = (numer >= denom ? 1 : numer/denom); +// t = 1 - s; +// } +// else +// { +// //min on edge s=0 +// s = 0; +// t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c)); +// } +// } +// else if (t < 0) +// { +// //region 6 +// const scalar tmp0 = b + d; +// const scalar tmp1 = c + e; +// if (tmp1 > tmp0) +// { +// //min on edge s+t=1 +// const scalar numer = tmp1 - tmp0; +// const scalar denom = a-2*b+c; +// s = (numer >= denom ? 1 : numer/denom); +// t = 1 - s; +// } +// else +// { +// //min on edge t=0 +// t = 0; +// s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a)); +// } +// } +// else +// { +// //region 1 +// const scalar numer = c+e-(b+d); +// if (numer <= 0) +// { +// s = 0; +// } +// else +// { +// const scalar denom = a-2*b+c; +// s = (numer >= denom ? 1 : numer/denom); +// } +// } +// t = 1 - s; +// } - // Calculate distance. - // Note: abs should not be needed but truncation error causes problems - // with points very close to one of the triangle vertices. - // (seen up to -9e-15). Alternatively add some small value. +// // Calculate distance. +// // Note: abs should not be needed but truncation error causes problems +// // with points very close to one of the triangle vertices. +// // (seen up to -9e-15). Alternatively add some small value. - const scalar f = D & D; - return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f); -} +// const scalar f = D & D; +// return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f); +// } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -234,9 +234,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType ( surface_, sample, - pHit.index(), - pHit.hitPoint(), - indexedOctree::perturbTol() + pHit.index() ); if (t == triSurfaceTools::UNKNOWN) @@ -353,39 +351,43 @@ void Foam::treeDataTriSurface::findNearest // ) //) { - // Get spanning vectors of triangle - vector base(p1); - vector E0(p0 - p1); - vector E1(p2 - p1); + // // Get spanning vectors of triangle + // vector base(p1); + // vector E0(p0 - p1); + // vector E1(p2 - p1); - scalar a(E0& E0); - scalar b(E0& E1); - scalar c(E1& E1); + // scalar a(E0& E0); + // scalar b(E0& E1); + // scalar c(E1& E1); - // Get nearest point in s,t coordinates (s is along E0, t is along - // E1) - scalar s; - scalar t; + // // Get nearest point in s,t coordinates (s is along E0, t is along + // // E1) + // scalar s; + // scalar t; - scalar distSqr = nearestCoords - ( - base, - E0, - E1, - a, - b, - c, - sample, + // scalar distSqr = nearestCoords + // ( + // base, + // E0, + // E1, + // a, + // b, + // c, + // sample, - s, - t - ); + // s, + // t + // ); + + pointHit pHit = triPointRef(p0, p1, p2).nearestPoint(sample); + + scalar distSqr = sqr(pHit.distance()); if (distSqr < nearestDistSqr) { nearestDistSqr = distSqr; minIndex = index; - nearestPoint = base + s*E0 + t*E1; + nearestPoint = pHit.rawPoint(); } } } diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.H b/src/meshTools/indexedOctree/treeDataTriSurface.H index f602783afe..dc959f78ff 100644 --- a/src/meshTools/indexedOctree/treeDataTriSurface.H +++ b/src/meshTools/indexedOctree/treeDataTriSurface.H @@ -60,20 +60,20 @@ class treeDataTriSurface // Private Member Functions - //- fast triangle nearest point calculation. Returns point in E0, E1 - // coordinate system: base + s*E0 + t*E1 - static scalar nearestCoords - ( - const point& base, - const point& E0, - const point& E1, - const scalar a, - const scalar b, - const scalar c, - const point& P, - scalar& s, - scalar& t - ); + // //- fast triangle nearest point calculation. Returns point in E0, E1 + // // coordinate system: base + s*E0 + t*E1 + // static scalar nearestCoords + // ( + // const point& base, + // const point& E0, + // const point& E1, + // const scalar a, + // const scalar b, + // const scalar c, + // const point& P, + // scalar& s, + // scalar& t + // ); public: diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C index 110f9dcd66..c8955f7470 100644 --- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C +++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C @@ -430,10 +430,7 @@ bool Foam::edgeIntersections::offsetPerturb point ctr = tri.centre(); - // Get measure for tolerance. - scalar tolDim = 0.001*mag(tri.a() - ctr); - - tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel); + tri.classify(pHit.hitPoint(), nearType, nearLabel); if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE) { diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C index 6602fbcf43..58c206aa75 100644 --- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C +++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C @@ -315,7 +315,7 @@ void Foam::surfaceIntersection::classifyHit surf2Pts[f2[0]], surf2Pts[f2[1]], surf2Pts[f2[2]] - ).classify(pHit.hitPoint(), tolDim, nearType, nearLabel); + ).classify(pHit.hitPoint(), nearType, nearLabel); // Classify points on edge of surface1 label edgeEnd = diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C index c80991819e..55e60a3111 100644 --- a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C +++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C @@ -43,7 +43,7 @@ Foam::scalar Foam::octreeDataTriSurface::tol(1E-6); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Fast distance to triangle calculation. From -// "Distance Between Point and Trangle in 3D" +// "Distance Between Point and Triangle in 3D" // David Eberly, Magic Software Inc. Aug. 2003. // Works on function Q giving distance to point and tries to minimize this. void Foam::octreeDataTriSurface::nearestCoords diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C index 6462df4cb3..38a7afd4b8 100644 --- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C +++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C @@ -234,9 +234,7 @@ void Foam::orientedSurface::propagateOrientation ( s, samplePoint, - nearestFaceI, - nearestPt, - 10*SMALL + nearestFaceI ); if (side == triSurfaceTools::UNKNOWN) diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 5a0cc03623..c2a88bf7e2 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -2121,12 +2121,13 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal label nearType; label nearLabel; + triPointRef ( points[f[0]], points[f[1]], points[f[2]] - ).classify(nearestPt, 1E-6, nearType, nearLabel); + ).classify(nearestPt, nearType, nearLabel); if (nearType == triPointRef::NONE) { @@ -2199,28 +2200,60 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide ( const triSurface& surf, const point& sample, - const label nearestFaceI, // nearest face - const point& nearestPoint, // nearest point on nearest face - const scalar tol + const label nearestFaceI ) { const labelledTri& f = surf[nearestFaceI]; const pointField& points = surf.points(); - // Find where point is on triangle. Note tolerance needed. Is relative - // one (used in comparing normalized [0..1] triangle coordinates). + // Find where point is on triangle. label nearType, nearLabel; - triPointRef + + pointHit pHit = triPointRef ( points[f[0]], points[f[1]], points[f[2]] - ).classify(nearestPoint, tol, nearType, nearLabel); + ).nearestPointClassify(sample, nearType, nearLabel); + + const point& nearestPoint(pHit.rawPoint()); if (nearType == triPointRef::NONE) { // Nearest to face interior. Use faceNormal to determine side +<<<<<<< HEAD scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI]; +======= + scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI]; + + // If the sample is essentially on the face, do not check for + // it being perpendicular. + + if (magSampleNearestVec > SMALL) + { + c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]); + + if (mag(c) < 0.99) + { + FatalErrorIn("triSurfaceTools::surfaceSide") + << "nearestPoint identified as being on triangle face " + << "but vector from nearestPoint to sample is not " + << "perpendicular to the normal." << nl + << "sample: " << sample << nl + << "nearestPoint: " << nearestPoint << nl + << "sample - nearestPoint: " << sample - nearestPoint << nl + << "normal: " << surf.faceNormals()[nearestFaceI] << nl + << "mag(sample - nearestPoint): " + << mag(sample - nearestPoint) << nl + << "normalised dot product: " << c << nl + << "triangle vertices: " << nl + << " " << points[f[0]] << nl + << " " << points[f[1]] << nl + << " " << points[f[2]] << nl + << abort(FatalError); + } + } +>>>>>>> 0bb6ebd... ENH: Making nearestPointClassify query for triangle. if (c > 0) { @@ -2239,27 +2272,27 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide // Get the edge. Assume order of faceEdges same as face vertices. label edgeI = surf.faceEdges()[nearestFaceI][nearLabel]; - //if (debug) - //{ - // // Check order of faceEdges same as face vertices. - // const edge& e = surf.edges()[edgeI]; - // const labelList& meshPoints = surf.meshPoints(); - // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]); - // - // if - // ( - // meshEdge - // != edge(f[nearLabel], f[f.fcIndex(nearLabel)]) - // ) - // { - // FatalErrorIn("triSurfaceTools::surfaceSide") - // << "Edge:" << edgeI << " local vertices:" << e - // << " mesh vertices:" << meshEdge - // << " not at position " << nearLabel - // << " in face " << f - // << abort(FatalError); - // } - //} + // if (debug) + { + // Check order of faceEdges same as face vertices. + const edge& e = surf.edges()[edgeI]; + const labelList& meshPoints = surf.meshPoints(); + const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]); + + if + ( + meshEdge + != edge(f[nearLabel], f[f.fcIndex(nearLabel)]) + ) + { + FatalErrorIn("triSurfaceTools::surfaceSide") + << "Edge:" << edgeI << " local vertices:" << e + << " mesh vertices:" << meshEdge + << " not at position " << nearLabel + << " in face " << f + << abort(FatalError); + } + } return edgeSide(surf, sample, nearestPoint, edgeI); } @@ -2717,7 +2750,14 @@ void Foam::triSurfaceTools::calcInterpolationWeights triPointRef tri(f.tri(points)); - pointHit nearest = tri.nearestPoint(samplePt); + label nearType, nearLabel; + + pointHit nearest = tri.nearestPointClassify + ( + samplePt, + nearType, + nearLabel + ); if (nearest.hit()) { @@ -2741,14 +2781,6 @@ void Foam::triSurfaceTools::calcInterpolationWeights minDistance = nearest.distance(); // Outside triangle. Store nearest. - label nearType, nearLabel; - tri.classify - ( - nearest.rawPoint(), - 1E-6, // relative tolerance - nearType, - nearLabel - ); if (nearType == triPointRef::POINT) { @@ -2779,12 +2811,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights max ( 0, - mag(nearest.rawPoint()-p0)/mag(p1-p0) + mag(nearest.rawPoint() - p0)/mag(p1 - p0) ) ); // Interpolate - weights[0] = 1-s; + weights[0] = 1 - s; weights[1] = s; weights[2] = -GREAT; @@ -2830,7 +2862,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify triPointRef(s[triI].tri(s.points())).classify ( trianglePoint, - 1E-6, elemType, index ); diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H index 62d66d031f..744a12217d 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H @@ -458,9 +458,7 @@ public: ( const triSurface& surf, const point& sample, - const label nearestFaceI, // nearest face - const point& nearestPt, // nearest point on nearest face - const scalar tol // tolerance for nearness test. + const label nearestFaceI ); // Triangulation of faces diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C index b777edbdb3..4a4c9219d6 100644 --- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C @@ -89,20 +89,29 @@ void Foam::distanceSurface::createGeometry() if (signed_) { - vectorField normal; - surfPtr_().getNormal(nearest, normal); + List volType; - forAll(nearest, i) + surfPtr_().getVolumeType(cc, volType); + + forAll(volType, i) { - vector d(cc[i]-nearest[i].hitPoint()); + searchableSurface::volumeType vT = volType[i]; - if ((d&normal[i]) > 0) + if (vT == searchableSurface::OUTSIDE) { - fld[i] = Foam::mag(d); + fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint()); + } + else if (vT == searchableSurface::INSIDE) + { + fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint()); } else { - fld[i] = -Foam::mag(d); + FatalErrorIn + ( + "void Foam::distanceSurface::createGeometry()" + ) << "getVolumeType failure, neither INSIDE or OUTSIDE" + << exit(FatalError); } } } @@ -132,20 +141,30 @@ void Foam::distanceSurface::createGeometry() if (signed_) { - vectorField normal; - surfPtr_().getNormal(nearest, normal); + List volType; - forAll(nearest, i) + surfPtr_().getVolumeType(cc, volType); + + forAll(volType, i) { - vector d(cc[i]-nearest[i].hitPoint()); + searchableSurface::volumeType vT = volType[i]; - if ((d&normal[i]) > 0) + if (vT == searchableSurface::OUTSIDE) { - fld[i] = Foam::mag(d); + fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint()); + } + else if (vT == searchableSurface::INSIDE) + { + fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint()); } else { - fld[i] = -Foam::mag(d); + FatalErrorIn + ( + "void Foam::distanceSurface::createGeometry()" + ) << "getVolumeType failure, " + << "neither INSIDE or OUTSIDE" + << exit(FatalError); } } } @@ -179,20 +198,31 @@ void Foam::distanceSurface::createGeometry() if (signed_) { - vectorField normal; - surfPtr_().getNormal(nearest, normal); + List volType; - forAll(nearest, i) + surfPtr_().getVolumeType(pts, volType); + + forAll(volType, i) { - vector d(pts[i]-nearest[i].hitPoint()); + searchableSurface::volumeType vT = volType[i]; - if ((d&normal[i]) > 0) + if (vT == searchableSurface::OUTSIDE) { - pointDistance_[i] = Foam::mag(d); + pointDistance_[i] = + Foam::mag(pts[i] - nearest[i].hitPoint()); + } + else if (vT == searchableSurface::INSIDE) + { + pointDistance_[i] = + -Foam::mag(pts[i] - nearest[i].hitPoint()); } else { - pointDistance_[i] = -Foam::mag(d); + FatalErrorIn + ( + "void Foam::distanceSurface::createGeometry()" + ) << "getVolumeType failure, neither INSIDE or OUTSIDE" + << exit(FatalError); } } } From e91c2d7d627fe4798371dd8273b65fe80b251098 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 6 Oct 2010 12:06:15 +0100 Subject: [PATCH 03/12] BUG: streamLine : store paths for stagnant particles. --- .../functionObjects/field/streamLine/streamLineParticle.C | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C index ce90b7a55d..56b05c88ca 100644 --- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C +++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C @@ -168,7 +168,6 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td) td.keepParticle && !td.switchProcessor && lifeTime_ > 0 - && tEnd > ROOTVSMALL ) { // TBD: implement subcycling so step through cells in more than @@ -191,6 +190,12 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td) tEnd -= dt; stepFraction() = 1.0 - tEnd/deltaT; + + if (tEnd <= ROOTVSMALL) + { + // Force removal + lifeTime_ = 0; + } } if (!td.keepParticle || lifeTime_ == 0) From 2b7f905fbd1389ae15331cc42c89550550915a34 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 6 Oct 2010 12:06:52 +0100 Subject: [PATCH 04/12] ENH: splitCyclics.txt : updated --- doc/changes/splitCyclic.txt | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/doc/changes/splitCyclic.txt b/doc/changes/splitCyclic.txt index 814dfe0850..3d0ffb71cf 100644 --- a/doc/changes/splitCyclic.txt +++ b/doc/changes/splitCyclic.txt @@ -20,7 +20,8 @@ The disadvantages: - a patch-wise loop now might need to store data to go to the neighbour half since it is no longer handled in a single patch. - decomposed cyclics now require overlapping communications so will -only work in non-blocking mode. Hence the underlying message passing library +only work in 'nonBlocking' mode or 'blocking' (=buffered) mode but not +in 'scheduled' mode. The underlying message passing library will require overlapping communications with message tags. - it is quite a code-change and there might be some oversights. - once converted (see foamUpgradeCyclics below) cases are not backwards @@ -103,19 +104,14 @@ type 'processorCyclic'. - processor patches use overlapping communication using a different message -tag. This maps straight through into the MPI message tag. -See processorCyclicPolyPatch::tag(). This needs to be calculated the -same on both sides so is calculated as - Pstream::nProcs()*max(myProcNo, neighbProcNo) - + min(myProcNo, neighbProcNo) -which is -- unique -- commutative -- does not interfere with the default tag (= 1) +tag. This maps straight through into the MPI message tag. Each processor +'interface' (processorPolyPatch, processorFvPatch, etc.) has a 'tag()' +to use for communication. - when constructing a GeometricField from a dictionary it will explicitly check for non-existing entries for cyclic patches and exit with an error message -warning to run foamUpgradeCyclics. +warning to run foamUpgradeCyclics. (1.7.x will check if you are trying +to run a case which has split cyclics) From 08a20c5bd19bdebb4c3c37f2e90cd5ebd12b2bbe Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 6 Oct 2010 12:45:42 +0100 Subject: [PATCH 05/12] BUG: indexedOctree : missing access function --- src/meshTools/indexedOctree/indexedOctree.C | 24 +++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C index 3f85828c09..9971bd97b4 100644 --- a/src/meshTools/indexedOctree/indexedOctree.C +++ b/src/meshTools/indexedOctree/indexedOctree.C @@ -2794,6 +2794,30 @@ Foam::indexedOctree::getVolumeType } +template +template +void Foam::indexedOctree::findNear +( + const scalar nearDist, + const indexedOctree& tree2, + CompareOp& cop +) const +{ + findNear + ( + nearDist, + true, + *this, + nodePlusOctant(0, 0), + bb(), + tree2, + nodePlusOctant(0, 0), + tree2.bb(), + cop + ); +} + + // Print contents of nodeI template void Foam::indexedOctree::print From 04c5454b211e333b835d544eb49316d100ed1617 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 6 Oct 2010 14:34:22 +0100 Subject: [PATCH 06/12] ENH: directMappedVelocityFlux : do not evaluate upon read --- .../directMappedVelocityFluxFixedValueFvPatchField.C | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C index 716c5b1235..17fa6f71fb 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C @@ -112,13 +112,6 @@ directMappedVelocityFluxFixedValueFvPatchField << " in file " << dimensionedInternalField().objectPath() << exit(FatalError); } - - // Force calculation of schedule (uses parallel comms) - const directMappedPolyPatch& mpp = refCast - ( - this->patch().patch() - ); - (void)mpp.map().schedule(); } From 557d91634253aa54b72b4d823201630b4b4c5e5a Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 6 Oct 2010 15:34:07 +0100 Subject: [PATCH 07/12] ENH: face::nearestPointClassify, similar to triangle. --- src/OpenFOAM/meshes/meshShapes/face/face.H | 22 +++++ .../meshes/meshShapes/face/faceIntersection.C | 55 ++++++++++- .../indexedOctree/treeDataTriSurface.C | 4 +- .../triSurfaceTools/triSurfaceTools.C | 96 +++++++++---------- 4 files changed, 125 insertions(+), 52 deletions(-) diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index 81ce90d8a9..dd5fc4a832 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/src/OpenFOAM/meshes/meshShapes/face/face.H @@ -126,6 +126,14 @@ class face public: + //- Return types for classify + enum proxType + { + NONE, + POINT, // Close to point + EDGE // Close to edge + }; + // Static data members static const char* const typeName; @@ -249,6 +257,20 @@ public: const pointField& meshPoints ) const; + //- Return nearest point to face and classify it: + // + near point (nearType=POINT, nearLabel=0, 1, 2) + // + near edge (nearType=EDGE, nearLabel=0, 1, 2) + // Note: edges are counted from starting vertex so + // e.g. edge n is from f[n] to f[0], where the face has n + 1 + // points + pointHit nearestPointClassify + ( + const point& p, + const pointField& meshPoints, + label& nearType, + label& nearLabel + ) const; + //- Return contact sphere diameter scalar contactSphereDiameter ( diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C index d4184e766a..f9f02eb612 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C @@ -181,6 +181,22 @@ Foam::pointHit Foam::face::nearestPoint const point& p, const pointField& meshPoints ) const +{ + // Dummy labels + label nearType = -1; + label nearLabel = -1; + + return nearestPointClassify(p, meshPoints, nearType, nearLabel); +} + + +Foam::pointHit Foam::face::nearestPointClassify +( + const point& p, + const pointField& meshPoints, + label& nearType, + label& nearLabel +) const { const face& f = *this; point ctr = centre(meshPoints); @@ -188,6 +204,9 @@ Foam::pointHit Foam::face::nearestPoint // Initialize to miss, distance=GREAT pointHit nearest(p); + nearType = -1; + nearLabel = -1; + label nPoints = f.size(); point nextPoint = ctr; @@ -196,8 +215,10 @@ Foam::pointHit Foam::face::nearestPoint { nextPoint = meshPoints[f[fcIndex(pI)]]; + label tmpNearType = -1; + label tmpNearLabel = -1; + // Note: for best accuracy, centre point always comes last - // triPointRef tri ( meshPoints[f[pI]], @@ -205,12 +226,42 @@ Foam::pointHit Foam::face::nearestPoint ctr ); - pointHit curHit = tri.nearestPoint(p); + pointHit curHit = tri.nearestPointClassify + ( + p, + tmpNearType, + tmpNearLabel + ); if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance())) { nearest.setDistance(curHit.distance()); + // Assume at first that the near type is NONE on the + // triangle (i.e. on the face of the triangle) then it is + // therefore also for the face. + + nearType = NONE; + + if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0) + { + // If the triangle edge label is 0, then this is also + // an edge of the face, if not, it is on the face + + nearType = EDGE; + + nearLabel = pI; + } + else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2) + { + // If the triangle point label is 0 or 1, then this is + // also a point of the face, if not, it is on the face + + nearType = POINT; + + nearLabel = pI + tmpNearLabel; + } + if (curHit.hit()) { nearest.setHit(); diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C index b019a85304..f380167174 100644 --- a/src/meshTools/indexedOctree/treeDataTriSurface.C +++ b/src/meshTools/indexedOctree/treeDataTriSurface.C @@ -360,8 +360,8 @@ void Foam::treeDataTriSurface::findNearest // scalar b(E0& E1); // scalar c(E1& E1); - // // Get nearest point in s,t coordinates (s is along E0, t is along - // // E1) + // // Get nearest point in s,t coordinates (s is along E0, t + // // is along E1) // scalar s; // scalar t; diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index c2a88bf7e2..8ade80cbf8 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -2220,40 +2220,40 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide if (nearType == triPointRef::NONE) { + vector sampleNearestVec = (sample - nearestPoint); + // Nearest to face interior. Use faceNormal to determine side -<<<<<<< HEAD - scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI]; -======= scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI]; - // If the sample is essentially on the face, do not check for - // it being perpendicular. + // // If the sample is essentially on the face, do not check for + // // it being perpendicular. - if (magSampleNearestVec > SMALL) - { - c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]); + // scalar magSampleNearestVec = mag(sampleNearestVec); - if (mag(c) < 0.99) - { - FatalErrorIn("triSurfaceTools::surfaceSide") - << "nearestPoint identified as being on triangle face " - << "but vector from nearestPoint to sample is not " - << "perpendicular to the normal." << nl - << "sample: " << sample << nl - << "nearestPoint: " << nearestPoint << nl - << "sample - nearestPoint: " << sample - nearestPoint << nl - << "normal: " << surf.faceNormals()[nearestFaceI] << nl - << "mag(sample - nearestPoint): " - << mag(sample - nearestPoint) << nl - << "normalised dot product: " << c << nl - << "triangle vertices: " << nl - << " " << points[f[0]] << nl - << " " << points[f[1]] << nl - << " " << points[f[2]] << nl - << abort(FatalError); - } - } ->>>>>>> 0bb6ebd... ENH: Making nearestPointClassify query for triangle. + // if (magSampleNearestVec > SMALL) + // { + // c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]); + + // if (mag(c) < 0.99) + // { + // FatalErrorIn("triSurfaceTools::surfaceSide") + // << "nearestPoint identified as being on triangle face " + // << "but vector from nearestPoint to sample is not " + // << "perpendicular to the normal." << nl + // << "sample: " << sample << nl + // << "nearestPoint: " << nearestPoint << nl + // << "sample - nearestPoint: " << sample - nearestPoint << nl + // << "normal: " << surf.faceNormals()[nearestFaceI] << nl + // << "mag(sample - nearestPoint): " + // << mag(sample - nearestPoint) << nl + // << "normalised dot product: " << c << nl + // << "triangle vertices: " << nl + // << " " << points[f[0]] << nl + // << " " << points[f[1]] << nl + // << " " << points[f[2]] << nl + // << abort(FatalError); + // } + // } if (c > 0) { @@ -2273,26 +2273,26 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide label edgeI = surf.faceEdges()[nearestFaceI][nearLabel]; // if (debug) - { - // Check order of faceEdges same as face vertices. - const edge& e = surf.edges()[edgeI]; - const labelList& meshPoints = surf.meshPoints(); - const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]); + // { + // // Check order of faceEdges same as face vertices. + // const edge& e = surf.edges()[edgeI]; + // const labelList& meshPoints = surf.meshPoints(); + // const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]); - if - ( - meshEdge - != edge(f[nearLabel], f[f.fcIndex(nearLabel)]) - ) - { - FatalErrorIn("triSurfaceTools::surfaceSide") - << "Edge:" << edgeI << " local vertices:" << e - << " mesh vertices:" << meshEdge - << " not at position " << nearLabel - << " in face " << f - << abort(FatalError); - } - } + // if + // ( + // meshEdge + // != edge(f[nearLabel], f[f.fcIndex(nearLabel)]) + // ) + // { + // FatalErrorIn("triSurfaceTools::surfaceSide") + // << "Edge:" << edgeI << " local vertices:" << e + // << " mesh vertices:" << meshEdge + // << " not at position " << nearLabel + // << " in face " << f + // << abort(FatalError); + // } + // } return edgeSide(surf, sample, nearestPoint, edgeI); } From dfbcf6eb463e3f99b4534069b046a137679be7e0 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 6 Oct 2010 16:20:55 +0100 Subject: [PATCH 08/12] BUG: singleCellFvMesh : locally (but not globally) zero sized patches --- .../singleCellFvMesh/singleCellFvMesh.C | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C index 4955f49007..8eaf3559f5 100644 --- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C +++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C @@ -42,24 +42,26 @@ void Foam::singleCellFvMesh::agglomerateMesh const polyBoundaryMesh& oldPatches = mesh.boundaryMesh(); // Check agglomeration within patch face range and continuous - labelList nAgglom(oldPatches.size()); + labelList nAgglom(oldPatches.size(), 0); forAll(oldPatches, patchI) { const polyPatch& pp = oldPatches[patchI]; - - nAgglom[patchI] = max(agglom[patchI])+1; - - forAll(pp, i) + if (pp.size() > 0) { - if (agglom[patchI][i] < 0 || agglom[patchI][i] >= pp.size()) + nAgglom[patchI] = max(agglom[patchI])+1; + + forAll(pp, i) { - FatalErrorIn - ( - "singleCellFvMesh::agglomerateMesh(..)" - ) << "agglomeration on patch " << patchI - << " is out of range 0.." << pp.size()-1 - << exit(FatalError); + if (agglom[patchI][i] < 0 || agglom[patchI][i] >= pp.size()) + { + FatalErrorIn + ( + "singleCellFvMesh::agglomerateMesh(..)" + ) << "agglomeration on patch " << patchI + << " is out of range 0.." << pp.size()-1 + << exit(FatalError); + } } } } @@ -155,6 +157,8 @@ void Foam::singleCellFvMesh::agglomerateMesh forAll(oldPatches, patchI) { + patchStarts[patchI] = coarseI; + const polyPatch& pp = oldPatches[patchI]; if (pp.size() > 0) @@ -170,8 +174,6 @@ void Foam::singleCellFvMesh::agglomerateMesh // From agglomeration to compact patch face labelList agglomToFace(nAgglom[patchI], -1); - patchStarts[patchI] = coarseI; - forAll(pp, i) { label myAgglom = agglom[patchI][i]; @@ -223,9 +225,9 @@ void Foam::singleCellFvMesh::agglomerateMesh ); } } - - patchSizes[patchI] = coarseI-patchStarts[patchI]; } + + patchSizes[patchI] = coarseI-patchStarts[patchI]; } //Pout<< "patchStarts:" << patchStarts << endl; From 95a2b3effda64a5cfcfe957911e74a9752bf3051 Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 6 Oct 2010 17:28:13 +0100 Subject: [PATCH 09/12] STYLE: 80 char lines. --- src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 8ade80cbf8..af30b3bcc8 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -2242,7 +2242,8 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide // << "perpendicular to the normal." << nl // << "sample: " << sample << nl // << "nearestPoint: " << nearestPoint << nl - // << "sample - nearestPoint: " << sample - nearestPoint << nl + // << "sample - nearestPoint: " + // << sample - nearestPoint << nl // << "normal: " << surf.faceNormals()[nearestFaceI] << nl // << "mag(sample - nearestPoint): " // << mag(sample - nearestPoint) << nl From aadacc5812c7c39be87d98246ac546bb414ba848 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 6 Oct 2010 17:41:40 +0100 Subject: [PATCH 10/12] BUG: corrected laminar clause in kappt wall function --- ...ppatJayatillekeWallFunctionFvPatchScalarField.C | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C index 7aa9fa2fe9..6efef7364b 100644 --- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C @@ -231,19 +231,17 @@ void kappatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs() scalar P = Psmooth(Prat); scalar yPlusTherm = this->yPlusTherm(P, Prat); - // Evaluate new effective thermal diffusivity - scalar kappaEff = 0.0; - if (yPlus < yPlusTherm) + // Update turbulent thermal conductivity + if (yPlus > yPlusTherm) { - kappaEff = Pr*yPlus; + scalar nu = nuw[faceI]; + scalar kt = nu*(yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P) - 1/Pr); + kappatw[faceI] = max(0.0, kt); } else { - kappaEff = nuw[faceI]*yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P); + kappatw[faceI] = 0.0; } - - // Update turbulent thermal diffusivity - kappatw[faceI] = max(0.0, kappaEff - nuw[faceI]/Pr); } fixedValueFvPatchField::updateCoeffs(); From 590e2c3a5253a9b5b69dca192a408cc78e388f73 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 6 Oct 2010 17:47:53 +0100 Subject: [PATCH 11/12] ENH: Updated heat transfer solvers and tutorials to use p_rgh --- .../buoyantBoussinesqPimpleFoam/TEqn.H | 2 +- .../buoyantBoussinesqPimpleFoam/UEqn.H | 9 +- .../buoyantBoussinesqPimpleFoam.C | 2 +- .../createFields.H | 56 ++- .../buoyantBoussinesqPimpleFoam/pEqn.H | 34 +- .../createFields.H | 18 +- .../buoyantBoussinesqSimpleFoam/pEqn.H | 15 +- .../heatTransfer/buoyantPimpleFoam/UEqn.H | 9 +- .../buoyantPimpleFoam/buoyantPimpleFoam.C | 2 +- .../buoyantPimpleFoam/createFields.H | 27 +- .../heatTransfer/buoyantPimpleFoam/hEqn.H | 2 +- .../heatTransfer/buoyantPimpleFoam/pEqn.H | 38 +- .../buoyantSimpleFoam/buoyantSimpleFoam.C | 3 - .../buoyantSimpleFoam/createFields.H | 68 ++- .../heatTransfer/buoyantSimpleFoam/pEqn.H | 9 +- .../buoyantSimpleRadiationFoam.C | 2 +- .../chtMultiRegionFoam/chtMultiRegionFoam.C | 2 + .../chtMultiRegionSimpleFoam/fluid/UEqn.H | 6 +- .../fluid/createFluidFields.H | 59 ++- .../chtMultiRegionSimpleFoam/fluid/hEqn.H | 4 +- .../chtMultiRegionSimpleFoam/fluid/pEqn.H | 53 ++- .../fluid/setRegionFluidFields.H | 5 +- .../fluid/solveFluid.H | 2 +- .../chtMultiRegionFoam/fluid/UEqn.H | 9 +- .../fluid/createFluidFields.H | 41 +- .../chtMultiRegionFoam/fluid/hEqn.H | 13 +- .../chtMultiRegionFoam/fluid/pEqn.H | 42 +- .../fluid/readFluidMultiRegionPISOControls.H | 17 - .../fluid/setRegionFluidFields.H | 14 +- .../chtMultiRegionFoam/fluid/solveFluid.H | 10 + .../fluid/storeOldFluidFields.H | 2 +- .../solid/createSolidFields.H | 6 +- .../chtMultiRegionFoam/solid/solveSolid.H | 12 +- .../lagrangian/coalChemistryFoam/UEqn.H | 2 +- .../coalChemistryFoam/coalChemistryFoam.C | 15 +- .../lagrangian/coalChemistryFoam/hsEqn.H | 2 +- .../lagrangian/coalChemistryFoam/pEqn.H | 30 +- .../hotRoom/0/T.org | 406 +----------------- .../buoyantBoussinesqPimpleFoam/hotRoom/0/p | 15 +- .../hotRoom/0/p_rgh | 45 ++ .../hotRoom/Allclean | 2 +- .../hotRoom/Allrun | 1 + .../hotRoom/system/fvSchemes | 16 +- .../hotRoom/system/fvSolution | 8 +- .../buoyantBoussinesqSimpleFoam/hotRoom/0/p | 15 +- .../hotRoom/Allclean | 2 +- .../hotRoom/Allrun | 1 + .../hotRoom/system/controlDict | 20 + .../hotRoom/system/fvSolution | 6 +- .../iglooWithFridges/0/p | 20 +- .../constant/polyMesh/boundary | 18 +- .../iglooWithFridges/system/controlDict | 20 + .../iglooWithFridges/system/fvSolution | 7 +- .../buoyantPimpleFoam/hotRoom/0/p | 12 +- .../buoyantPimpleFoam/hotRoom/0/p_rgh | 42 ++ .../hotRoom/system/fvSchemes | 4 +- .../hotRoom/system/fvSolution | 12 +- .../buoyantSimpleFoam/buoyantCavity/0/p | 16 +- .../buoyantCavity/system/controlDict | 20 + .../buoyantCavity/system/fvSolution | 11 +- .../buoyantSimpleFoam/hotRoom/0/T | 406 +----------------- .../buoyantSimpleFoam/hotRoom/0/T.org | 406 +----------------- .../buoyantSimpleFoam/hotRoom/0/p | 12 +- .../buoyantSimpleFoam/hotRoom/0/p_rgh | 42 ++ .../hotRoom/system/controlDict | 20 + .../hotRoom/system/fvSchemes | 14 +- .../hotRoom/system/fvSolution | 16 +- .../hotRadiationRoom/0/G | 8 +- .../hotRadiationRoom/0/p | 16 +- .../hotRadiationRoom/0/p_rgh | 48 +++ .../hotRadiationRoom/system/controlDict | 20 + .../hotRadiationRoom/system/fvSchemes | 4 +- .../hotRadiationRoom/system/fvSolution | 14 +- .../hotRadiationRoomFvDOM/0/IDefault | 2 +- .../hotRadiationRoomFvDOM/0/p | 8 +- .../hotRadiationRoomFvDOM/0/p_rgh | 48 +++ .../hotRadiationRoomFvDOM/system/controlDict | 20 + .../hotRadiationRoomFvDOM/system/fvSchemes | 4 +- .../hotRadiationRoomFvDOM/system/fvSolution | 9 +- .../multiRegionHeater/0/epsilon | 1 - .../chtMultiRegionFoam/multiRegionHeater/0/k | 1 - .../multiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 19 +- .../multiRegionHeater/system/controlDict | 9 +- .../system/heater/decomposeParDict | 2 +- .../system/heater/fvSolution | 6 + .../system/leftSolid/fvSolution | 6 + .../system/rightSolid/fvSolution | 6 + .../system/topAir/changeDictionaryDict | 39 +- .../multiRegionHeater/system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 43 +- .../snappyMultiRegionHeater/0/alphat | 1 - .../snappyMultiRegionHeater/0/epsilon | 1 - .../snappyMultiRegionHeater/0/k | 1 - .../snappyMultiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/decomposeParDict | 2 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 42 +- .../system/heater/decomposeParDict | 2 +- .../system/heater/fvSolution | 6 + .../system/leftSolid/decomposeParDict | 2 +- .../system/leftSolid/fvSolution | 7 + .../system/rightSolid/decomposeParDict | 2 +- .../system/rightSolid/fvSolution | 7 + .../system/topAir/changeDictionaryDict | 39 +- .../system/topAir/decomposeParDict | 2 +- .../system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 43 +- .../multiRegionHeater/0/U | 2 +- .../multiRegionHeater/0/epsilon | 1 - .../multiRegionHeater/0/k | 1 - .../multiRegionHeater/0/p_rgh | 29 ++ .../system/bottomAir/changeDictionaryDict | 20 +- .../system/bottomAir/decomposeParDict | 2 +- .../system/bottomAir/fvSchemes | 6 +- .../system/bottomAir/fvSolution | 13 +- .../multiRegionHeater/system/controlDict | 7 +- .../multiRegionHeater/system/decomposeParDict | 2 +- .../system/heater/decomposeParDict | 2 +- .../system/leftSolid/decomposeParDict | 2 +- .../system/rightSolid/decomposeParDict | 2 +- .../system/topAir/changeDictionaryDict | 34 +- .../system/topAir/decomposeParDict | 2 +- .../multiRegionHeater/system/topAir/fvSchemes | 6 +- .../system/topAir/fvSolution | 15 +- 128 files changed, 1255 insertions(+), 1780 deletions(-) delete mode 100644 applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H create mode 100644 tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p_rgh create mode 100644 tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/p_rgh create mode 100644 tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/p_rgh diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H index dbfc61739f..9a835792a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H @@ -12,7 +12,7 @@ ); TEqn.relax(); - TEqn.solve(); + TEqn.solve(mesh.solver(T.select(finalIter))); rhok = 1.0 - beta*(T - TRef); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H index 35387f4ff4..20a05e5cd4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H @@ -18,9 +18,10 @@ fvc::reconstruct ( ( - fvc::interpolate(rhok)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) - ) + - ghf*fvc::snGrad(rhok) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C index ebf68d409c..54519906a4 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C @@ -87,7 +87,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H index 23d20cfa96..342af079d8 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H @@ -14,12 +14,12 @@ mesh ); - Info<< "Reading field p\n" << endl; - volScalarField p + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( IOobject ( - "p", + "p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -52,6 +52,18 @@ incompressible::RASModel::New(U, phi, laminarTransport) ); + // Kinematic density for buoyancy force + volScalarField rhok + ( + IOobject + ( + "rhok", + runTime.timeName(), + mesh + ), + 1.0 - beta*(T - TRef) + ); + // kinematic turbulent thermal thermal conductivity m2/s Info<< "Reading field kappat\n" << endl; volScalarField kappat @@ -67,25 +79,41 @@ mesh ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rhok*gh + ); + label pRefCell = 0; scalar pRefValue = 0.0; setRefCell ( p, + p_rgh, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue ); - - // Kinematic density for buoyancy force - volScalarField rhok - ( - IOobject + if (p_rgh.needReference()) + { + p += dimensionedScalar ( - "rhok", - runTime.timeName(), - mesh - ), - 1.0 - beta*(T - TRef) - ); + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + } + diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H index 21be033f9b..60828e18dc 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H @@ -7,22 +7,23 @@ phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); - surfaceScalarField buoyancyPhi = - rUAf*fvc::interpolate(rhok)*(g & mesh.Sf()); - phi += buoyancyPhi; + surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf(); + phi -= buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvm::laplacian(rUAf, p) == fvc::div(phi) + fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - pEqn.solve + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -36,17 +37,30 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi -= pEqn.flux(); + phi -= p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf); + U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf); U.correctBoundaryConditions(); } } #include "continuityErrs.H" + + p = p_rgh + rhok*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - rhok*gh; + } } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H index 76cc4da8ba..477a322833 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H @@ -96,29 +96,23 @@ p_rgh + rhok*gh ); - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; + label pRefCell = 0; + scalar pRefValue = 0.0; setRefCell ( + p, p_rgh, mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + pRefCell, + pRefValue ); - scalar pRefValue = 0.0; - if (p_rgh.needReference()) { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - p += dimensionedScalar ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); } diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H index 5664bb2154..50149ed360 100644 --- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H @@ -18,17 +18,9 @@ fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); - // retain the residual from the first iteration - if (nonOrth == 0) - { - p_rghEqn.solve(); - } - else - { - p_rghEqn.solve(); - } + p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) { @@ -55,7 +47,8 @@ ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) + pRefValue - getRefCellValue(p, pRefCell) ); + p_rgh = p - rhok*gh; } } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H index d4878d063d..8c6a3f7671 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H @@ -17,8 +17,11 @@ == fvc::reconstruct ( - fvc::interpolate(rho)*(g & mesh.Sf()) - - fvc::snGrad(p)*mesh.magSf() - ) + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ), + mesh.solver(U.select(finalIter)) ); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C index 0075c1e3ed..cb4c4d34f6 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C @@ -80,7 +80,7 @@ int main(int argc, char *argv[]) if (nOuterCorr != 1) { - p.storePrevIter(); + p_rgh.storePrevIter(); } #include "UEqn.H" diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H index b8ac5595e4..e39d0bb803 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H @@ -53,15 +53,30 @@ ) ); + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + Info<< "Creating field DpDt\n" << endl; volScalarField DpDt ( "DpDt", fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p) ); - - thermo.correct(); - - dimensionedScalar initialMass = fvc::domainIntegrate(rho); - - dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H index 3125cc3ffa..94537508b3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H @@ -9,7 +9,7 @@ ); hEqn.relax(); - hEqn.solve(); + hEqn.solve(mesh.solver(h.select(finalIter))); thermo.correct(); } diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H index a1c3dbfed5..e625f122a3 100644 --- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H @@ -1,11 +1,9 @@ { - bool closedVolume = p.needReference(); - rho = thermo.rho(); // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution - done in 2 parts. Part 1: - thermo.rho() -= psi*p; + thermo.rho() -= psi*p_rgh; volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); @@ -18,24 +16,23 @@ + fvc::ddtPhiCorr(rUA, rho, U, phi) ); - surfaceScalarField buoyancyPhi = - rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf()); + surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi += buoyancyPhi; for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn + fvScalarMatrix p_rghEqn ( - fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phi) - - fvm::laplacian(rhorUAf, p) + - fvm::laplacian(rhorUAf, p_rgh) ); - pEqn.solve + p_rghEqn.solve ( mesh.solver ( - p.select + p_rgh.select ( ( finalIter @@ -49,34 +46,25 @@ if (nonOrth == nNonOrthCorr) { // Calculate the conservative fluxes - phi += pEqn.flux(); + phi += p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector - p.relax(); + p_rgh.relax(); // Correct the momentum source with the pressure gradient flux // calculated from the relaxed pressure - U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf); + U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf); U.correctBoundaryConditions(); } } + p = p_rgh + rho*gh; + // Second part of thermodynamic density update - thermo.rho() += psi*p; + thermo.rho() += psi*p_rgh; DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); #include "rhoEqn.H" #include "compressibleContinuityErrs.H" - - // For closed-volume cases adjust the pressure and density levels - // to obey overall mass continuity - if (closedVolume) - { - p += - (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); - thermo.rho() = psi*p; - rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume; - } } diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C index 06fa5d12a0..64cc110cec 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C @@ -62,10 +62,7 @@ int main(int argc, char *argv[]) { #include "UEqn.H" #include "hEqn.H" - for (int i=0; i<3; i++) - { #include "pEqn.H" - } } turbulence->correct(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H index cf23150332..d6fa9acee9 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H @@ -23,20 +23,6 @@ volScalarField& h = thermo.h(); const volScalarField& psi = thermo.psi(); - Info<< "Reading field p_rgh\n" << endl; - volScalarField p_rgh - ( - IOobject - ( - "p_rgh", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - Info<< "Reading field U\n" << endl; volVectorField U ( @@ -53,7 +39,6 @@ #include "compressibleCreatePhi.H" - Info<< "Creating turbulence model\n" << endl; autoPtr turbulence ( @@ -66,40 +51,39 @@ ) ); + Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("ghf", g & mesh.Cf()); - p = p_rgh + rho*gh; - thermo.correct(); - rho = thermo.rho(); - p_rgh = p - rho*gh; - - label p_rghRefCell = 0; - scalar p_rghRefValue = 0.0; - setRefCell + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( - p_rgh, - mesh.solutionDict().subDict("SIMPLE"), - p_rghRefCell, - p_rghRefValue + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh ); + // Force p_rgh to be consistent with p + p_rgh = p - rho*gh; + + + label pRefCell = 0; scalar pRefValue = 0.0; - - if (p_rgh.needReference()) - { - pRefValue = readScalar - ( - mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue") - ); - - p += dimensionedScalar - ( - "p", - p.dimensions(), - pRefValue - getRefCellValue(p, p_rghRefCell) - ); - } + setRefCell + ( + p, + p_rgh, + mesh.solutionDict().subDict("SIMPLE"), + pRefCell, + pRefValue + ); dimensionedScalar initialMass = fvc::domainIntegrate(rho); + dimensionedScalar totalVolume = sum(mesh.V()); diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H index 3088b42c69..f1a62d4770 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H +++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H @@ -1,11 +1,12 @@ { rho = thermo.rho(); + rho.relax(); volScalarField rUA = 1.0/UEqn().A(); surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA)); U = rUA*UEqn().H(); - //UEqn.clear(); + UEqn.clear(); phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); bool closedVolume = adjustPhi(phi, U, p_rgh); @@ -20,7 +21,7 @@ fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(); if (nonOrth == nNonOrthCorr) @@ -42,13 +43,13 @@ p = p_rgh + rho*gh; - // For closed-volume cases adjust the pressure and density levels + // For closed-volume cases adjust the pressure level // to obey overall mass continuity if (closedVolume) { p += (initialMass - fvc::domainIntegrate(psi*p)) /fvc::domainIntegrate(psi); - p_rgh == p - rho*gh; + p_rgh = p - rho*gh; } rho = thermo.rho(); diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C index 0dbe80c67c..6c35536333 100644 --- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C +++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) #include "readSIMPLEControls.H" - p.storePrevIter(); + p_rgh.storePrevIter(); rho.storePrevIter(); // Pressure-velocity SIMPLE corrector diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index a81b0faaf3..36b3f1c3b0 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -93,6 +93,8 @@ int main(int argc, char *argv[]) // --- PIMPLE loop for (int oCorr=0; oCorr phiFluid(fluidRegions.size()); PtrList gFluid(fluidRegions.size()); PtrList turbulence(fluidRegions.size()); - PtrList DpDtf(fluidRegions.size()); + PtrList p_rghFluid(fluidRegions.size()); + PtrList ghFluid(fluidRegions.size()); + PtrList ghfFluid(fluidRegions.size()); List initialMassFluid(fluidRegions.size()); List