diff --git a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options index 14c5e10a3e..abee4af7eb 100644 --- a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options +++ b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options @@ -19,7 +19,6 @@ EXE_LIBS = \ -ldynamicMesh \ -ledgeMesh \ -lengine \ - -lextrudeModel \ -lfieldFunctionObjects \ -lfileFormats \ -lfiniteVolume \ @@ -43,7 +42,6 @@ EXE_LIBS = \ -lliquidMixtureProperties \ -lliquidProperties \ -lmeshTools \ - -lMGridGenGAMGAgglomeration \ -lmolecularMeasurements \ -lmolecule \ -lmultiphaseInterFoam \ diff --git a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C index 4392395af0..c9ab1ffbbe 100644 --- a/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C +++ b/applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -100,7 +100,7 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer ( new passiveParticle ( - positions_, + procMesh, ppi.position(), procCelli, false @@ -112,7 +112,7 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer particleIndices_.setSize(pi); - IOPosition(positions_).write(); + IOPosition >(positions_).write(); } diff --git a/applications/utilities/preProcessing/mapFields/mapLagrangian.C b/applications/utilities/preProcessing/mapFields/mapLagrangian.C index 36dd981719..97a999b62c 100644 --- a/applications/utilities/preProcessing/mapFields/mapLagrangian.C +++ b/applications/utilities/preProcessing/mapFields/mapLagrangian.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,8 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "MapLagrangianFields.H" -#include "Cloud.H" -#include "passiveParticle.H" +#include "passiveParticleCloud.H" #include "meshSearch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,7 +43,9 @@ static label findCell(const Cloud& cloud, const point& pt) label tetFaceI = -1; label tetPtI = -1; - cloud.findCellFacePt(pt, cellI, tetFaceI, tetPtI); + const polyMesh& mesh = cloud.pMesh(); + + mesh.findCellFacePt(pt, cellI, tetFaceI, tetPtI); if (cellI >= 0) { @@ -55,8 +56,6 @@ static label findCell(const Cloud& cloud, const point& pt) // See if particle on face by finding nearest face and shifting // particle. - const polyMesh& mesh = cloud.pMesh(); - meshSearch meshSearcher(mesh, false); label faceI = meshSearcher.findNearestBoundaryFace(pt); @@ -67,7 +66,7 @@ static label findCell(const Cloud& cloud, const point& pt) const point perturbPt = (1-perturbFactor)*pt+perturbFactor*cc; - cloud.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI); + mesh.findCellFacePt(perturbPt, cellI, tetFaceI, tetPtI); return cellI; } @@ -124,7 +123,7 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp) Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl; // Read positions & cell - Cloud sourceParcels + passiveParticleCloud sourceParcels ( meshSource, cloudDirs[cloudI], @@ -134,13 +133,15 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp) << " parcels from source mesh." << endl; // Construct empty target cloud - Cloud targetParcels + passiveParticleCloud targetParcels ( meshTarget, cloudDirs[cloudI], IDLList() ); + particle::TrackingData td(targetParcels); + label sourceParticleI = 0; // Indices of source particles that get added to targetParcels @@ -176,15 +177,14 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp) ( new passiveParticle ( - targetParcels, + meshTarget, targetCc[targetCells[i]], targetCells[i] ) ); passiveParticle& newP = newPtr(); - scalar fraction = 0; - label faceI = newP.track(iter().position(), fraction); + label faceI = newP.track(iter().position(), td); if (faceI < 0 && newP.cell() >= 0) { @@ -246,7 +246,7 @@ void mapLagrangian(const meshToMesh& meshToMeshInterp) if (addParticles.size()) { - IOPosition(targetParcels).write(); + IOPosition(targetParcels).write(); // addParticles now contains the indices of the sourceMesh // particles that were appended to the target mesh. diff --git a/src/Allwmake b/src/Allwmake index e9a922ab36..37313b73df 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -38,9 +38,9 @@ parallel/decompose/AllwmakeLnInclude # dummyThirdParty (dummy metisDecomp, scotchDecomp etc) needed by e.g. meshTools dummyThirdParty/Allwmake $* +wmake $makeOption finiteVolume wmake $makeOption lagrangian/basic wmake $makeOption lagrangian/distributionModels -wmake $makeOption finiteVolume wmake $makeOption genericPatchFields # Build the proper scotchDecomp, metisDecomp etc. diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 18522689e8..432f6ccc0e 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -495,6 +495,8 @@ $(pointBoundaryMesh)/pointBoundaryMesh.C meshes/boundBox/boundBox.C +meshes/treeBoundBox/treeBoundBox.C + meshTools = meshes/meshTools $(meshTools)/matchPoints.C $(meshTools)/mergePoints.C @@ -572,6 +574,11 @@ $(interpolations)/patchToPatchInterpolation/PatchToPatchInterpolationName.C algorithms/MeshWave/MeshWaveName.C algorithms/MeshWave/FaceCellWaveName.C +algorithms/indexedOctree/indexedOctreeName.C +algorithms/indexedOctree/treeDataCell.C + + + graph/curve/curve.C graph/graph.C diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C similarity index 98% rename from src/meshTools/indexedOctree/indexedOctree.C rename to src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C index df489a6f95..6020d24d19 100644 --- a/src/meshTools/indexedOctree/indexedOctree.C +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C @@ -25,8 +25,8 @@ License #include "indexedOctree.H" #include "linePointRef.H" -#include "meshTools.H" #include "OFstream.H" +#include "ListOps.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -1730,26 +1730,6 @@ void Foam::indexedOctree::traverseNode point perturbedEnd(pushPoint(octantBb, end, false)); - - //if (debug) - { - // Dump octantBb to obj - writeOBJ(nodeI, octant); - // Dump ray to obj as well - { - OFstream str("ray.obj"); - meshTools::writeOBJ(str, start); - meshTools::writeOBJ(str, end); - str << "l 1 2" << nl; - } - WarningIn("indexedOctree::traverseNode(..)") - << "Did not intersect ray from endpoint:" << end - << " to startpoint:" << start - << " with bounding box:" << octantBb << nl - << "Re-intersecting with perturbed endpoint:" << perturbedEnd - << endl; - } - traverseNode ( findAny, @@ -2316,12 +2296,6 @@ void Foam::indexedOctree::writeOBJ pointField bbPoints(subBb.points()); label pointVertI = vertI; - forAll(bbPoints, i) - { - meshTools::writeOBJ(str, bbPoints[i]); - vertI++; - } - forAll(treeBoundBox::edges, i) { const edge& e = treeBoundBox::edges[i]; diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H similarity index 100% rename from src/meshTools/indexedOctree/indexedOctree.H rename to src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H diff --git a/src/meshTools/indexedOctree/indexedOctreeName.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctreeName.C similarity index 95% rename from src/meshTools/indexedOctree/indexedOctreeName.C rename to src/OpenFOAM/algorithms/indexedOctree/indexedOctreeName.C index 519d0d5992..4483d42068 100644 --- a/src/meshTools/indexedOctree/indexedOctreeName.C +++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctreeName.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/meshTools/indexedOctree/labelBits.H b/src/OpenFOAM/algorithms/indexedOctree/labelBits.H similarity index 98% rename from src/meshTools/indexedOctree/labelBits.H rename to src/OpenFOAM/algorithms/indexedOctree/labelBits.H index 221881f5bb..58893cbaeb 100644 --- a/src/meshTools/indexedOctree/labelBits.H +++ b/src/OpenFOAM/algorithms/indexedOctree/labelBits.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C similarity index 100% rename from src/meshTools/indexedOctree/treeDataCell.C rename to src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C diff --git a/src/meshTools/indexedOctree/treeDataCell.H b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H similarity index 98% rename from src/meshTools/indexedOctree/treeDataCell.H rename to src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H index 06342311e8..7cebddafa9 100644 --- a/src/meshTools/indexedOctree/treeDataCell.H +++ b/src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/OpenFOAM/fields/cloud/cloud.C b/src/OpenFOAM/fields/cloud/cloud.C index c4260f02ee..40eb070bcd 100644 --- a/src/OpenFOAM/fields/cloud/cloud.C +++ b/src/OpenFOAM/fields/cloud/cloud.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,7 +41,7 @@ Foam::cloud::cloud(const objectRegistry& obr, const word& cloudName) ( IOobject ( - ( cloudName.size() ? cloudName : defaultName ), + (cloudName.size() ? cloudName : defaultName), obr.time().timeName(), prefix, obr, @@ -58,4 +58,12 @@ Foam::cloud::~cloud() {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::cloud::autoMap(const mapPolyMesh&) +{ + notImplemented("cloud::autoMap(const mapPolyMesh&)"); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/fields/cloud/cloud.H b/src/OpenFOAM/fields/cloud/cloud.H index a4c4eae51e..0755c6dd61 100644 --- a/src/OpenFOAM/fields/cloud/cloud.H +++ b/src/OpenFOAM/fields/cloud/cloud.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -91,7 +91,7 @@ public: //- Remap the cells of particles corresponding to the // mesh topology change - virtual void autoMap(const mapPolyMesh&) = 0; + virtual void autoMap(const mapPolyMesh&); }; diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index f83861aae9..da7ffdc266 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -1183,10 +1183,123 @@ void Foam::polyMesh::removeFiles(const fileName& instanceDir) const } } + void Foam::polyMesh::removeFiles() const { removeFiles(instance()); } +void Foam::polyMesh::findCellFacePt +( + const point& pt, + label& cellI, + label& tetFaceI, + label& tetPtI +) const +{ + cellI = -1; + tetFaceI = -1; + tetPtI = -1; + + const indexedOctree& tree = cellTree(); + + // Find nearest cell to the point + + pointIndexHit info = tree.findNearest(pt, sqr(GREAT)); + + if (info.hit()) + { + label nearestCellI = tree.shapes().cellLabels()[info.index()]; + + // Check the nearest cell to see if the point is inside. + findTetFacePt(nearestCellI, pt, tetFaceI, tetPtI); + + if (tetFaceI != -1) + { + // Point was in the nearest cell + + cellI = nearestCellI; + + return; + } + else + { + // Check the other possible cells that the point may be in + + labelList testCells = tree.findIndices(pt); + + forAll(testCells, pCI) + { + label testCellI = tree.shapes().cellLabels()[testCells[pCI]]; + + if (testCellI == nearestCellI) + { + // Don't retest the nearest cell + + continue; + } + + // Check the test cell to see if the point is inside. + findTetFacePt(testCellI, pt, tetFaceI, tetPtI); + + if (tetFaceI != -1) + { + // Point was in the test cell + + cellI = testCellI; + + return; + } + } + } + } + else + { + FatalErrorIn + ( + "void Foam::polyMesh::findCellFacePt" + "(" + "const point&, " + "label&, " + "label&, " + "label&" + ") const" + ) << "Did not find nearest cell in search tree." + << abort(FatalError); + } +} + + +void Foam::polyMesh::findTetFacePt +( + const label cellI, + const point& pt, + label& tetFaceI, + label& tetPtI +) const +{ + const polyMesh& mesh = *this; + + tetFaceI = -1; + tetPtI = -1; + + List cellTets = + polyMeshTetDecomposition::cellTetIndices(mesh, cellI); + + forAll(cellTets, tetI) + { + const tetIndices& cellTetIs = cellTets[tetI]; + + if (cellTetIs.tet(mesh).inside(pt)) + { + tetFaceI = cellTetIs.face(); + tetPtI = cellTetIs.tetPt(); + + return; + } + } +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H index 65cf88dbe5..2be783a25e 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,6 +34,7 @@ SourceFiles polyMeshFromShapeMesh.C polyMeshIO.C polyMeshUpdate.C + polyMeshFindCell.C \*---------------------------------------------------------------------------*/ @@ -391,6 +392,7 @@ public: return *this; } + // Mesh motion //- Is mesh moving @@ -508,6 +510,28 @@ public: //- Remove all files from mesh instance() void removeFiles() const; + + + // Helper functions + + //- Find the cell, tetFaceI and tetPtI for the given position + void findCellFacePt + ( + const point& pt, + label& cellI, + label& tetFaceI, + label& tetPtI + ) const; + + //- Find the tetFaceI and tetPtI for the given position in + // the supplied cell, tetFaceI and tetPtI = -1 if not found + void findTetFacePt + ( + const label cellI, + const point& pt, + label& tetFaceI, + label& tetPtI + ) const; }; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C index a0fd190585..b0e83ce670 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,8 +26,6 @@ License #include "primitiveMesh.H" #include "demandDrivenData.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::primitiveMesh, 0); @@ -63,6 +61,8 @@ Foam::primitiveMesh::primitiveMesh() ppPtr_(NULL), cpPtr_(NULL), + cellTreePtr_(NULL), + labels_(0), cellCentresPtr_(NULL), @@ -105,6 +105,8 @@ Foam::primitiveMesh::primitiveMesh ppPtr_(NULL), cpPtr_(NULL), + cellTreePtr_(NULL), + labels_(0), cellCentresPtr_(NULL), @@ -347,4 +349,36 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const } +const Foam::indexedOctree& +Foam::primitiveMesh::cellTree() const +{ + if (!cellTreePtr_) + { + treeBoundBox overallBb(points()); + + Random rndGen(261782); + + overallBb = overallBb.extend(rndGen, 1E-4); + overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + + cellTreePtr_ = + new indexedOctree + ( + treeDataCell + ( + false, // not cache bb + *this + ), + overallBb, + 8, // maxLevel + 10, // leafsize + 3.0 // duplicity + ); + } + + return *cellTreePtr_; +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H index 23eefb0e0b..cc3bc2c80a 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,6 +67,8 @@ SourceFiles #include "Map.H" #include "EdgeMap.H" #include "boundBox.H" +#include "indexedOctree.H" +#include "treeDataCell.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -155,6 +157,9 @@ class primitiveMesh //- Cell-points mutable labelListList* cpPtr_; + //- Search tree to allow spatial tet searching + mutable indexedOctree* cellTreePtr_; + // On-the-fly edge addresing storage @@ -482,6 +487,10 @@ public: const labelListList& cellPoints() const; + //- Build (if necessary) and return the cell search tree + const indexedOctree& cellTree() const; + + // Geometric data (raw!) const vectorField& cellCentres() const; @@ -814,6 +823,9 @@ public: //- Clear topological data void clearAddressing(); + //- Clear cell tree data + void clearCellTree(); + //- Clear all geometry and addressing unnecessary for CFD void clearOut(); }; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C index f11a901c7d..cc232c162f 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -99,6 +99,11 @@ void Foam::primitiveMesh::printAllocated() const Pout<< " Cell-point" << endl; } + if (cellTreePtr_) + { + Pout<< " Cell-tree" << endl; + } + // Geometry if (cellCentresPtr_) { @@ -165,6 +170,14 @@ void Foam::primitiveMesh::clearAddressing() deleteDemandDrivenData(pePtr_); deleteDemandDrivenData(ppPtr_); deleteDemandDrivenData(cpPtr_); + + deleteDemandDrivenData(cellTreePtr_); +} + + +void Foam::primitiveMesh::clearCellTree() +{ + deleteDemandDrivenData(cellTreePtr_); } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C index a42befb41a..c7087aa70f 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,10 +26,8 @@ License #include "primitiveMesh.H" #include "cell.H" - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -// Is the point in the cell bounding box bool Foam::primitiveMesh::pointInCellBB ( const point& p, @@ -60,7 +58,6 @@ bool Foam::primitiveMesh::pointInCellBB } -// Is the point in the cell bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const { const labelList& f = cells()[celli]; @@ -86,7 +83,6 @@ bool Foam::primitiveMesh::pointInCell(const point& p, label celli) const } -// Find the cell with the nearest cell centre Foam::label Foam::primitiveMesh::findNearestCell(const point& location) const { const vectorField& centres = cellCentres(); @@ -109,7 +105,6 @@ Foam::label Foam::primitiveMesh::findNearestCell(const point& location) const } -// Find cell enclosing this location Foam::label Foam::primitiveMesh::findCell(const point& location) const { if (nCells() == 0) diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index 7e55cd872c..53e92c7f66 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -174,10 +174,10 @@ public: ) const; //- Return nearest point to p on tetrahedron - inline pointHit nearestPoint - ( - const point& p - ) const; + inline pointHit nearestPoint(const point& p) const; + + //- Return true if point is inside tetrahedron + inline bool inside(const point& pt) const; //- Return (min)containment sphere, i.e. the smallest sphere with // all points inside. Returns pointHit with: diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 38c9eae41c..0051951e6b 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -380,6 +380,82 @@ inline Foam::pointHit Foam::tetrahedron::nearestPoint } +template +bool Foam::tetrahedron::inside(const point& pt) const +{ + // For robustness, assuming that the point is in the tet unless + // "definitively" shown otherwise by obtaining a positive dot + // product greater than a tolerance of SMALL. + + // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal + // vectors and base points for the half-space planes are: + // area[0] = Sa(); + // area[1] = Sb(); + // area[2] = Sc(); + // area[3] = Sd(); + // planeBase[0] = tetBasePt = b_ + // planeBase[1] = ptA = c_ + // planeBase[2] = tetBasePt = b_ + // planeBase[3] = tetBasePt = b_ + + vector n = vector::zero; + + { + // 0, a + const point& basePt = b_; + + n = Sa(); + n /= (Foam::mag(n) + VSMALL); + + if (((pt - basePt) & n) > SMALL) + { + return false; + } + } + + { + // 1, b + const point& basePt = c_; + + n = Sb(); + n /= (Foam::mag(n) + VSMALL); + + if (((pt - basePt) & n) > SMALL) + { + return false; + } + } + + { + // 2, c + const point& basePt = b_; + + n = Sc(); + n /= (Foam::mag(n) + VSMALL); + + if (((pt - basePt) & n) > SMALL) + { + return false; + } + } + + { + // 3, d + const point& basePt = b_; + + n = Sd(); + n /= (Foam::mag(n) + VSMALL); + + if (((pt - basePt) & n) > SMALL) + { + return false; + } + } + + return true; +} + + // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // template diff --git a/src/meshTools/octree/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C similarity index 100% rename from src/meshTools/octree/treeBoundBox.C rename to src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C diff --git a/src/meshTools/octree/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H similarity index 100% rename from src/meshTools/octree/treeBoundBox.H rename to src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H diff --git a/src/meshTools/octree/treeBoundBoxI.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H similarity index 99% rename from src/meshTools/octree/treeBoundBoxI.H rename to src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H index 7426e04708..1f069feb48 100644 --- a/src/meshTools/octree/treeBoundBoxI.H +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/meshTools/octree/treeBoundBoxList.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxList.H similarity index 96% rename from src/meshTools/octree/treeBoundBoxList.H rename to src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxList.H index 610df2c1b1..b82eede2f0 100644 --- a/src/meshTools/octree/treeBoundBoxList.H +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxList.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/meshTools/octree/treeBoundBoxTemplates.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxTemplates.C similarity index 96% rename from src/meshTools/octree/treeBoundBoxTemplates.C rename to src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxTemplates.C index c8cb5fab5d..09e0a9be7a 100644 --- a/src/meshTools/octree/treeBoundBoxTemplates.C +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C index c140c449d5..425c8221cd 100644 --- a/src/lagrangian/basic/Cloud/Cloud.C +++ b/src/lagrangian/basic/Cloud/Cloud.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -32,12 +32,6 @@ License #include "OFstream.H" #include "wallPolyPatch.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -template -const Foam::scalar Foam::Cloud::trackingCorrectionTol = 1e-5; - - // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // template @@ -47,7 +41,7 @@ void Foam::Cloud::calcCellWallFaces() const PackedBoolList& cellWallFaces = cellWallFacesPtr_(); - const polyBoundaryMesh& patches = pMesh().boundaryMesh(); + const polyBoundaryMesh& patches = polyMesh_.boundaryMesh(); forAll(patches, patchI) { @@ -78,9 +72,7 @@ Foam::Cloud::Cloud cloud(pMesh), IDLList(), polyMesh_(pMesh), - particleCount_(0), labels_(), - cellTree_(), nTrackingRescues_(), cellWallFacesPtr_() { @@ -104,9 +96,7 @@ Foam::Cloud::Cloud cloud(pMesh, cloudName), IDLList(), polyMesh_(pMesh), - particleCount_(0), labels_(), - cellTree_(), nTrackingRescues_(), cellWallFacesPtr_() { @@ -121,236 +111,6 @@ Foam::Cloud::Cloud // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template -void Foam::Cloud::findCellFacePt -( - const point& pt, - label& cellI, - label& tetFaceI, - label& tetPtI -) const -{ - cellI = -1; - tetFaceI = -1; - tetPtI = -1; - - const indexedOctree& tree = cellTree(); - - // Find nearest cell to the point - - pointIndexHit info = tree.findNearest(pt, sqr(GREAT)); - - if (info.hit()) - { - label nearestCellI = tree.shapes().cellLabels()[info.index()]; - - // Check the nearest cell to see if the point is inside. - findFacePt(nearestCellI, pt, tetFaceI, tetPtI); - - if (tetFaceI != -1) - { - // Point was in the nearest cell - - cellI = nearestCellI; - - return; - } - else - { - // Check the other possible cells that the point may be in - - labelList testCells = tree.findIndices(pt); - - forAll(testCells, pCI) - { - label testCellI = tree.shapes().cellLabels()[testCells[pCI]]; - - if (testCellI == nearestCellI) - { - // Don't retest the nearest cell - - continue; - } - - // Check the test cell to see if the point is inside. - findFacePt(testCellI, pt, tetFaceI, tetPtI); - - if (tetFaceI != -1) - { - // Point was in the test cell - - cellI = testCellI; - - return; - } - } - } - } - else - { - FatalErrorIn - ( - "void Foam::Cloud::findCellFacePt" - "(" - "const point& pt, " - "label& cellI, " - "label& tetFaceI, " - "label& tetPtI" - ") const" - ) << "Did not find nearest cell in search tree." - << abort(FatalError); - } -} - - -template -void Foam::Cloud::findFacePt -( - label cellI, - const point& pt, - label& tetFaceI, - label& tetPtI -) const -{ - tetFaceI = -1; - tetPtI = -1; - - List cellTets = polyMeshTetDecomposition::cellTetIndices - ( - polyMesh_, - cellI - ); - - forAll(cellTets, tetI) - { - const tetIndices& cellTetIs = cellTets[tetI]; - - if (inTet(pt, cellTetIs.tet(polyMesh_))) - { - tetFaceI = cellTetIs.face(); - tetPtI = cellTetIs.tetPt(); - - return; - } - } -} - - -template -bool Foam::Cloud::inTet -( - const point& pt, - const tetPointRef& tet -) const -{ - // For robustness, assuming that the point is in the tet unless - // "definitively" shown otherwise by obtaining a positive dot - // product greater than a tolerance of SMALL. - - // The tet is defined: tet(Cc, tetBasePt, pA, pB) where the normal - // vectors and base points for the half-space planes are: - // area[0] = tet.Sa(); - // area[1] = tet.Sb(); - // area[2] = tet.Sc(); - // area[3] = tet.Sd(); - // planeBase[0] = tetBasePt = tet.b() - // planeBase[1] = ptA = tet.c() - // planeBase[2] = tetBasePt = tet.b() - // planeBase[3] = tetBasePt = tet.b() - - vector n = vector::zero; - - { - // 0, a - const point& basePt = tet.b(); - - n = tet.Sa(); - n /= (mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 1, b - const point& basePt = tet.c(); - - n = tet.Sb(); - n /= (mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 2, c - const point& basePt = tet.b(); - - n = tet.Sc(); - n /= (mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - { - // 3, d - const point& basePt = tet.b(); - - n = tet.Sd(); - n /= (mag(n) + VSMALL); - - if (((pt - basePt) & n) > SMALL) - { - return false; - } - } - - return true; -} - - -template -const Foam::indexedOctree& -Foam::Cloud::cellTree() const -{ - if (cellTree_.empty()) - { - treeBoundBox overallBb(polyMesh_.points()); - - Random rndGen(261782); - - overallBb = overallBb.extend(rndGen, 1E-4); - overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); - overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); - - cellTree_.reset - ( - new indexedOctree - ( - treeDataCell - ( - false, // not cache bb - polyMesh_ - ), - overallBb, - 8, // maxLevel - 10, // leafsize - 3.0 // duplicity - ) - ); - } - - return cellTree_(); -} - - template const Foam::PackedBoolList& Foam::Cloud::cellHasWallFaces() const @@ -364,22 +124,6 @@ const } - -template -Foam::label Foam::Cloud::getNewParticleID() const -{ - label id = particleCount_++; - - if (id == labelMax) - { - WarningIn("Cloud::getNewParticleID() const") - << "Particle counter has overflowed. This might cause problems" - << " when reconstructing particle tracks." << endl; - } - return id; -} - - template void Foam::Cloud::addParticle(ParticleType* pPtr) { @@ -399,14 +143,14 @@ void Foam::Cloud::cloudReset(const Cloud& c) { // Reset particle cound and particles only // - not changing the cloud object registry or reference to the polyMesh - particleCount_ = 0; + ParticleType::particleCount_ = 0; IDLList::operator=(c); } template -template -void Foam::Cloud::move(TrackingData& td, const scalar trackTime) +template +void Foam::Cloud::move(TrackData& td, const scalar trackTime) { const polyBoundaryMesh& pbm = pMesh().boundaryMesh(); const globalMeshData& pData = polyMesh_.globalData(); @@ -472,9 +216,9 @@ void Foam::Cloud::move(TrackingData& td, const scalar trackTime) { // If we are running in parallel and the particle is on a // boundary face - if (Pstream::parRun() && p.faceI_ >= pMesh().nInternalFaces()) + if (Pstream::parRun() && p.face() >= pMesh().nInternalFaces()) { - label patchI = pbm.whichPatch(p.faceI_); + label patchI = pbm.whichPatch(p.face()); // ... and the face is on a processor patch // prepare it for transfer @@ -571,7 +315,7 @@ void Foam::Cloud::move(TrackingData& td, const scalar trackTime) IDLList newParticles ( particleStream, - typename ParticleType::iNew(*this) + typename ParticleType::iNew(polyMesh_) ); label pI = 0; @@ -600,57 +344,64 @@ void Foam::Cloud::move(TrackingData& td, const scalar trackTime) template -void Foam::Cloud::autoMap(const mapPolyMesh& mapper) +template +void Foam::Cloud::autoMap +( + TrackData& td, + const mapPolyMesh& mapper +) { if (cloud::debug) { - Info<< "Cloud::autoMap(const morphFieldMapper& map) " - "for lagrangian cloud " << cloud::name() << endl; + Info<< "Cloud::autoMap(TrackData&, const mapPolyMesh&) " + << "for lagrangian cloud " << cloud::name() << endl; } const labelList& reverseCellMap = mapper.reverseCellMap(); const labelList& reverseFaceMap = mapper.reverseFaceMap(); // Reset stored data that relies on the mesh - cellTree_.clear(); +// polyMesh_.clearCellTree(); cellWallFacesPtr_.clear(); forAllIter(typename Cloud, *this, pIter) { - if (reverseCellMap[pIter().cellI_] >= 0) - { - pIter().cellI_ = reverseCellMap[pIter().cellI_]; + ParticleType& p = pIter(); - if (pIter().faceI_ >= 0 && reverseFaceMap[pIter().faceI_] >= 0) + if (reverseCellMap[p.cell()] >= 0) + { + p.cell() = reverseCellMap[p.cell()]; + + if (p.face() >= 0 && reverseFaceMap[p.face()] >= 0) { - pIter().faceI_ = reverseFaceMap[pIter().faceI_]; + p.face() = reverseFaceMap[p.face()]; } else { - pIter().faceI_ = -1; + p.face() = -1; } - pIter().initCellFacePt(); + p.initCellFacePt(); } else { - label trackStartCell = mapper.mergedCell(pIter().cellI_); + label trackStartCell = mapper.mergedCell(p.cell()); if (trackStartCell < 0) { trackStartCell = 0; } - vector p = pIter().position(); + vector pos = p.position(); - const_cast(pIter().position()) = + const_cast(p.position()) = polyMesh_.cellCentres()[trackStartCell]; - pIter().stepFraction() = 0; + p.stepFraction() = 0; - pIter().initCellFacePt(); + p.initCellFacePt(); - pIter().track(p); + p.track(pos, td); } } } diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H index 3715bcb75c..488fc63f3d 100644 --- a/src/lagrangian/basic/Cloud/Cloud.H +++ b/src/lagrangian/basic/Cloud/Cloud.H @@ -25,6 +25,7 @@ Class Foam::Cloud Description + Base cloud calls templated on particle type SourceFiles Cloud.C @@ -80,15 +81,9 @@ class Cloud const polyMesh& polyMesh_; - //- Overall count of particles ever created. Never decreases. - mutable label particleCount_; - //- Temporary storage for addressing. Used in findTris. mutable DynamicList