From 23cf5d94cd937aae0a912f01da023f43a9a09d2c Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 3 Nov 2009 14:33:01 +0000 Subject: [PATCH 01/29] code tidying --- .../vtkPV3Foam/vtkPV3FoamFaceField.H | 21 ++++++++++--------- .../vtkPV3Foam/vtkPV3FoamMeshLagrangian.C | 4 ++-- .../vtkPV3Foam/vtkPV3FoamMeshPatch.C | 13 ++++-------- .../vtkPV3Foam/vtkPV3FoamMeshSet.C | 12 +++++------ .../vtkPV3Foam/vtkPV3FoamMeshVolume.C | 8 +++---- .../vtkPV3Foam/vtkPV3FoamMeshZone.C | 10 ++++----- .../vtkPV3Foam/vtkPV3FoamPatchField.H | 18 ++++++++-------- .../vtkPV3Foam/vtkPV3FoamPointFields.H | 6 +++--- .../vtkPV3Foam/vtkPV3FoamUpdateInfo.C | 4 +--- .../vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H | 3 ++- .../vtkPV3Foam/vtkPV3FoamUtils.C | 1 + 11 files changed, 45 insertions(+), 55 deletions(-) diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H index dba8d97d37..14c0d5f9b9 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H @@ -35,6 +35,7 @@ InClass #include "vtkFloatArray.h" #include "vtkMultiBlockDataSet.h" #include "vtkPolyData.h" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // template @@ -53,11 +54,11 @@ void Foam::vtkPV3Foam::convertFaceField const labelList& faceOwner = mesh.faceOwner(); const labelList& faceNeigh = mesh.faceNeighbour(); - vtkFloatArray *cellData = vtkFloatArray::New(); - cellData->SetNumberOfTuples( faceLabels.size() ); - cellData->SetNumberOfComponents( nComp ); - cellData->Allocate( nComp*faceLabels.size() ); - cellData->SetName( tf.name().c_str() ); + vtkFloatArray* cellData = vtkFloatArray::New(); + cellData->SetNumberOfTuples(faceLabels.size()); + cellData->SetNumberOfComponents(nComp); + cellData->Allocate(nComp*faceLabels.size()); + cellData->SetName(tf.name().c_str()); if (debug) { @@ -123,11 +124,11 @@ void Foam::vtkPV3Foam::convertFaceField const labelList& faceOwner = mesh.faceOwner(); const labelList& faceNeigh = mesh.faceNeighbour(); - vtkFloatArray *cellData = vtkFloatArray::New(); - cellData->SetNumberOfTuples( fSet.size() ); - cellData->SetNumberOfComponents( nComp ); - cellData->Allocate( nComp*fSet.size() ); - cellData->SetName( tf.name().c_str() ); + vtkFloatArray* cellData = vtkFloatArray::New(); + cellData->SetNumberOfTuples(fSet.size()); + cellData->SetNumberOfComponents(nComp); + cellData->Allocate(nComp*fSet.size()); + cellData->SetName(tf.name().c_str()); if (debug) { diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshLagrangian.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshLagrangian.C index e653b05adb..0d46db1d49 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshLagrangian.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshLagrangian.C @@ -80,8 +80,8 @@ vtkPolyData* Foam::vtkPV3Foam::lagrangianVTKMesh vtkPoints* vtkpoints = vtkPoints::New(); vtkCellArray* vtkcells = vtkCellArray::New(); - vtkpoints->Allocate( parcels.size() ); - vtkcells->Allocate( parcels.size() ); + vtkpoints->Allocate(parcels.size()); + vtkcells->Allocate(parcels.size()); vtkIdType particleId = 0; forAllConstIter(Cloud, parcels, iter) diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshPatch.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshPatch.C index d7c0f2f013..ace25f527f 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshPatch.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshPatch.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "vtkPV3Foam.H" @@ -40,10 +38,7 @@ Description // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh -( - const polyPatch& p -) +vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh(const polyPatch& p) { vtkPolyData* vtkmesh = vtkPolyData::New(); @@ -56,8 +51,8 @@ vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh // Convert Foam mesh vertices to VTK const Foam::pointField& points = p.localPoints(); - vtkPoints *vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( points.size() ); + vtkPoints* vtkpoints = vtkPoints::New(); + vtkpoints->Allocate(points.size()); forAll(points, i) { vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]); @@ -71,7 +66,7 @@ vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh const faceList& faces = p.localFaces(); vtkCellArray* vtkcells = vtkCellArray::New(); - vtkcells->Allocate( faces.size() ); + vtkcells->Allocate(faces.size()); forAll(faces, faceI) { const face& f = faces[faceI]; diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshSet.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshSet.C index f363beae45..431df1f9c5 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshSet.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshSet.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "vtkPV3Foam.H" @@ -71,8 +69,8 @@ vtkPolyData* Foam::vtkPV3Foam::faceSetVTKMesh // Convert Foam mesh vertices to VTK const pointField& points = p.localPoints(); - vtkPoints *vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( points.size() ); + vtkPoints* vtkpoints = vtkPoints::New(); + vtkpoints->Allocate(points.size()); forAll(points, i) { vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]); @@ -84,7 +82,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceSetVTKMesh const faceList& faces = p.localFaces(); vtkCellArray* vtkcells = vtkCellArray::New(); - vtkcells->Allocate( faces.size() ); + vtkcells->Allocate(faces.size()); forAll(faces, faceI) { @@ -127,8 +125,8 @@ vtkPolyData* Foam::vtkPV3Foam::pointSetVTKMesh const pointField& meshPoints = mesh.points(); - vtkPoints *vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( pSet.size() ); + vtkPoints* vtkpoints = vtkPoints::New(); + vtkpoints->Allocate(pSet.size()); forAllConstIter(pointSet, pSet, iter) { diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C index e1fd59ccd0..b39dfd453f 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "vtkPV3Foam.H" @@ -136,8 +134,8 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh } // Convert Foam mesh vertices to VTK - vtkPoints *vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( mesh.nPoints() + nAddPoints ); + vtkPoints* vtkpoints = vtkPoints::New(); + vtkpoints->Allocate(mesh.nPoints() + nAddPoints); const Foam::pointField& points = mesh.points(); @@ -152,7 +150,7 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh Info<< "... converting cells" << endl; } - vtkmesh->Allocate( mesh.nCells() + nAddCells ); + vtkmesh->Allocate(mesh.nCells() + nAddCells); // Set counters for additional points and additional cells label addPointI = 0, addCellI = 0; diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C index 0a21310a63..9b8bcb956d 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C @@ -22,8 +22,6 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ #include "vtkPV3Foam.H" @@ -69,7 +67,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceZoneVTKMesh const pointField& points = p.localPoints(); vtkPoints* vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( points.size() ); + vtkpoints->Allocate(points.size()); forAll(points, i) { vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]); @@ -83,7 +81,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceZoneVTKMesh const faceList& faces = p.localFaces(); vtkCellArray* vtkcells = vtkCellArray::New(); - vtkcells->Allocate( faces.size() ); + vtkcells->Allocate(faces.size()); forAll(faces, faceI) { @@ -126,8 +124,8 @@ vtkPolyData* Foam::vtkPV3Foam::pointZoneVTKMesh const pointField& meshPoints = mesh.points(); - vtkPoints *vtkpoints = vtkPoints::New(); - vtkpoints->Allocate( pointLabels.size() ); + vtkPoints* vtkpoints = vtkPoints::New(); + vtkpoints->Allocate(pointLabels.size()); forAll(pointLabels, pointI) { diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPatchField.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPatchField.H index 7b7de4d022..bfa2d1ac46 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPatchField.H +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPatchField.H @@ -52,10 +52,10 @@ void Foam::vtkPV3Foam::convertPatchField const label nComp = pTraits::nComponents; vtkFloatArray* cellData = vtkFloatArray::New(); - cellData->SetNumberOfTuples( ptf.size() ); - cellData->SetNumberOfComponents( nComp ); - cellData->Allocate( nComp*ptf.size() ); - cellData->SetName( name.c_str() ); + cellData->SetNumberOfTuples(ptf.size()); + cellData->SetNumberOfComponents(nComp); + cellData->Allocate(nComp*ptf.size()); + cellData->SetName(name.c_str()); float vec[nComp]; forAll(ptf, i) @@ -91,11 +91,11 @@ void Foam::vtkPV3Foam::convertPatchPointField { const label nComp = pTraits::nComponents; - vtkFloatArray *pointData = vtkFloatArray::New(); - pointData->SetNumberOfTuples( pptf.size() ); - pointData->SetNumberOfComponents( nComp ); - pointData->Allocate( nComp*pptf.size() ); - pointData->SetName( name.c_str() ); + vtkFloatArray* pointData = vtkFloatArray::New(); + pointData->SetNumberOfTuples(pptf.size()); + pointData->SetNumberOfComponents(nComp); + pointData->Allocate(nComp*pptf.size()); + pointData->SetName(name.c_str()); float vec[nComp]; forAll(pptf, i) diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H index 3a609b2bde..5e3ae1df42 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H @@ -110,9 +110,9 @@ void Foam::vtkPV3Foam::convertPointFields ++partId ) { - const word patchName = getPartName(partId); + const word patchName = getPartName(partId); const label datasetNo = partDataset_[partId]; - const label patchId = patches.findPatchID(patchName); + const label patchId = patches.findPatchID(patchName); if (!partStatus_[partId] || datasetNo < 0 || patchId < 0) { @@ -187,7 +187,7 @@ void Foam::vtkPV3Foam::convertPointField nPoints = ptf.size(); } - vtkFloatArray *pointData = vtkFloatArray::New(); + vtkFloatArray* pointData = vtkFloatArray::New(); pointData->SetNumberOfTuples(nPoints + addPointCellLabels.size()); pointData->SetNumberOfComponents(nComp); pointData->Allocate(nComp*(nPoints + addPointCellLabels.size())); diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C index 59e445e341..5bb2143c2a 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C @@ -140,7 +140,6 @@ void Foam::vtkPV3Foam::updateInfoInternalMesh() Info<< " Foam::vtkPV3Foam::updateInfoInternalMesh" << endl; } - } @@ -440,14 +439,13 @@ void Foam::vtkPV3Foam::updateInfoLagrangianFields() << endl; } - vtkDataArraySelection *fieldSelection = + vtkDataArraySelection* fieldSelection = reader_->GetLagrangianFieldSelection(); // preserve the enabled selections stringList enabledEntries = getSelectedArrayEntries(fieldSelection); fieldSelection->RemoveAllArrays(); - // // TODO - currently only get fields from ONE cloud // have to decide if the second set of fields get mixed in // or dealt with separately diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H index 7e28460a46..4d35525efe 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H @@ -35,7 +35,7 @@ InClass template class patchType, class meshType> void Foam::vtkPV3Foam::updateInfoFields ( - vtkDataArraySelection *select + vtkDataArraySelection* select ) { if (debug) @@ -112,4 +112,5 @@ void Foam::vtkPV3Foam::updateInfoFields // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif + // ************************************************************************* // diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUtils.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUtils.C index a9d12c303a..9683a15914 100644 --- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUtils.C +++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUtils.C @@ -69,6 +69,7 @@ namespace Foam } // End namespace Foam + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::vtkPV3Foam::AddToBlock From e7347dbd624989eea8d9b8f4614128f0c8f711bb Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 9 Nov 2009 10:15:07 +0000 Subject: [PATCH 02/29] removed unnecessary interpolation of h - not used --- .../solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C index bd1989f433..757ad95e25 100644 --- a/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C +++ b/applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C @@ -90,7 +90,6 @@ int main(int argc, char *argv[]) // --- PISO loop for (int corr=0; corr Date: Wed, 11 Nov 2009 17:19:18 +0000 Subject: [PATCH 03/29] incorrect reading of optional wedgePlane --- .../surfaceDisplacementPointPatchVectorField.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C index 72d0541e7e..f6d91a9edc 100644 --- a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C +++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C @@ -329,7 +329,7 @@ surfaceDisplacementPointPatchVectorField surfacesDict_(dict.subDict("geometry")), projectMode_(projectModeNames_.read(dict.lookup("projectMode"))), projectDir_(dict.lookup("projectDirection")), - wedgePlane_(dict.lookupOrDefault(dict.lookup("wedgePlane"), -1)), + wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)), frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null)) { if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0) From d7e880289647c3bad384aaa0907239002aa7ea79 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 11 Nov 2009 17:19:52 +0000 Subject: [PATCH 04/29] removed octreeDataFaceList; created generic PrimitivePatch search one --- src/dynamicMesh/Make/files | 1 - src/dynamicMesh/boundaryMesh/boundaryMesh.C | 138 +++-- .../boundaryMesh/octreeDataFaceList.C | 573 ------------------ .../boundaryMesh/octreeDataFaceList.H | 220 ------- src/meshTools/Make/files | 1 + .../indexedOctree/treeDataPrimitivePatch.C | 567 +++++++++++++++++ .../indexedOctree/treeDataPrimitivePatch.H | 210 +++++++ .../treeDataPrimitivePatchName.C | 33 + 8 files changed, 888 insertions(+), 855 deletions(-) delete mode 100644 src/dynamicMesh/boundaryMesh/octreeDataFaceList.C delete mode 100644 src/dynamicMesh/boundaryMesh/octreeDataFaceList.H create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatch.C create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatch.H create mode 100644 src/meshTools/indexedOctree/treeDataPrimitivePatchName.C diff --git a/src/dynamicMesh/Make/files b/src/dynamicMesh/Make/files index e59bbe3800..8201b7236d 100644 --- a/src/dynamicMesh/Make/files +++ b/src/dynamicMesh/Make/files @@ -49,7 +49,6 @@ perfectInterface/perfectInterface.C boundaryMesh/boundaryMesh.C boundaryPatch/boundaryPatch.C -boundaryMesh/octreeDataFaceList.C setUpdater/setUpdater.C meshModifiers = meshCut/meshModifiers diff --git a/src/dynamicMesh/boundaryMesh/boundaryMesh.C b/src/dynamicMesh/boundaryMesh/boundaryMesh.C index 494cda9295..b71e0a569e 100644 --- a/src/dynamicMesh/boundaryMesh/boundaryMesh.C +++ b/src/dynamicMesh/boundaryMesh/boundaryMesh.C @@ -29,11 +29,12 @@ License #include "polyMesh.H" #include "repatchPolyTopoChanger.H" #include "faceList.H" -#include "octree.H" -#include "octreeDataFaceList.H" +#include "indexedOctree.H" +#include "treeDataPrimitivePatch.H" #include "triSurface.H" #include "SortableList.H" #include "OFstream.H" +#include "uindirectPrimitivePatch.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -892,6 +893,17 @@ Foam::labelList Foam::boundaryMesh::getNearest << endl; } + uindirectPrimitivePatch leftPatch + ( + UIndirectList(mesh(), leftFaces), + mesh().points() + ); + uindirectPrimitivePatch rightPatch + ( + UIndirectList(mesh(), rightFaces), + mesh().points() + ); + // Overall bb treeBoundBox overallBb(mesh().localPoints()); @@ -911,29 +923,35 @@ Foam::labelList Foam::boundaryMesh::getNearest bbMax.z() += 2*tol; // Create the octrees - octree leftTree + indexedOctree + < + treeDataPrimitivePatch + > leftTree ( - overallBb, - octreeDataFaceList + treeDataPrimitivePatch ( - mesh(), - leftFaces + false, // cacheBb + leftPatch ), - 1, - 20, - 10 + overallBb, + 10, // maxLevel + 10, // leafSize + 3.0 // duplicity ); - octree rightTree + indexedOctree + < + treeDataPrimitivePatch + > rightTree ( - overallBb, - octreeDataFaceList + treeDataPrimitivePatch ( - mesh(), - rightFaces + false, // cacheBb + rightPatch ), - 1, - 20, - 10 + overallBb, + 10, // maxLevel + 10, // leafSize + 3.0 // duplicity ); if (debug) @@ -953,7 +971,7 @@ Foam::labelList Foam::boundaryMesh::getNearest treeBoundBox tightest; - const scalar searchDim = mag(searchSpan); + const scalar searchDimSqr = magSqr(searchSpan); forAll(nearestBFaceI, patchFaceI) { @@ -982,50 +1000,25 @@ Foam::labelList Foam::boundaryMesh::getNearest } // Search right tree - tightest.min() = ctr - searchSpan; - tightest.max() = ctr + searchSpan; - scalar rightDist = searchDim; - label rightI = rightTree.findNearest(ctr, tightest, rightDist); - + pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr); // Search left tree. Note: could start from rightDist bounding box // instead of starting from top. - tightest.min() = ctr - searchSpan; - tightest.max() = ctr + searchSpan; - scalar leftDist = searchDim; - label leftI = leftTree.findNearest(ctr, tightest, leftDist); + pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr); - - if (rightI == -1) + if (rightInfo.hit()) { - // No face found in right tree. - - if (leftI == -1) - { - // No face found in left tree. - nearestBFaceI[patchFaceI] = -1; - } - else - { - // Found in left but not in right. Choose left regardless if - // correct sign. Note: do we want this? - nearestBFaceI[patchFaceI] = leftFaces[leftI]; - } - } - else - { - if (leftI == -1) - { - // Found in right but not in left. Choose right regardless if - // correct sign. Note: do we want this? - nearestBFaceI[patchFaceI] = rightFaces[rightI]; - } - else + if (leftInfo.hit()) { // Found in both trees. Compare normals. + label rightFaceI = rightFaces[rightInfo.index()]; + label leftFaceI = leftFaces[leftInfo.index()]; - scalar rightSign = n & ns[rightFaces[rightI]]; - scalar leftSign = n & ns[leftFaces[leftI]]; + label rightDist = mag(rightInfo.hitPoint()-ctr); + label leftDist = mag(leftInfo.hitPoint()-ctr); + + scalar rightSign = n & ns[rightFaceI]; + scalar leftSign = n & ns[leftFaceI]; if ( @@ -1036,11 +1029,11 @@ Foam::labelList Foam::boundaryMesh::getNearest // Both same sign. Choose nearest. if (rightDist < leftDist) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } else @@ -1059,11 +1052,11 @@ Foam::labelList Foam::boundaryMesh::getNearest // Different sign and nearby. Choosing matching normal if (rightSign > 0) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } else @@ -1071,15 +1064,38 @@ Foam::labelList Foam::boundaryMesh::getNearest // Different sign but faraway. Choosing nearest. if (rightDist < leftDist) { - nearestBFaceI[patchFaceI] = rightFaces[rightI]; + nearestBFaceI[patchFaceI] = rightFaceI; } else { - nearestBFaceI[patchFaceI] = leftFaces[leftI]; + nearestBFaceI[patchFaceI] = leftFaceI; } } } } + else + { + // Found in right but not in left. Choose right regardless if + // correct sign. Note: do we want this? + label rightFaceI = rightFaces[rightInfo.index()]; + nearestBFaceI[patchFaceI] = rightFaceI; + } + } + else + { + // No face found in right tree. + + if (leftInfo.hit()) + { + // Found in left but not in right. Choose left regardless if + // correct sign. Note: do we want this? + nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()]; + } + else + { + // No face found in left tree. + nearestBFaceI[patchFaceI] = -1; + } } } diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C deleted file mode 100644 index c44c455c74..0000000000 --- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.C +++ /dev/null @@ -1,573 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - -\*---------------------------------------------------------------------------*/ - -#include "octreeDataFaceList.H" -#include "octree.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -defineTypeNameAndDebug(Foam::octreeDataFaceList, 0); - -Foam::scalar Foam::octreeDataFaceList::tol = 1E-6; - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::octreeDataFaceList::calcBb() -{ - allBb_.setSize(faceLabels_.size()); - allBb_ = treeBoundBox - ( - vector(GREAT, GREAT, GREAT), - vector(-GREAT, -GREAT, -GREAT) - ); - - forAll (faceLabels_, faceLabelI) - { - label faceI = faceLabels_[faceLabelI]; - - // Update bb of face - treeBoundBox& myBb = allBb_[faceLabelI]; - - const face& f = mesh_.localFaces()[faceI]; - - forAll(f, fp) - { - const point& coord = mesh_.localPoints()[f[fp]]; - - myBb.min() = min(myBb.min(), coord); - myBb.max() = max(myBb.max(), coord); - } - } -} - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from all faces in bMesh -Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh) -: - mesh_(mesh), - faceLabels_(mesh_.size()), - allBb_(mesh_.size()) -{ - forAll(faceLabels_, faceI) - { - faceLabels_[faceI] = faceI; - } - - // Generate tight fitting bounding box - calcBb(); -} - - -// Construct from selected faces in bMesh -Foam::octreeDataFaceList::octreeDataFaceList -( - const bMesh& mesh, - const labelList& faceLabels -) -: - mesh_(mesh), - faceLabels_(faceLabels), - allBb_(faceLabels.size()) -{ - // Generate tight fitting bounding box - calcBb(); -} - - - -// Construct as copy -Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes) -: - mesh_(shapes.mesh()), - faceLabels_(shapes.faceLabels()), - allBb_(shapes.allBb()) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::octreeDataFaceList::~octreeDataFaceList() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -Foam::label Foam::octreeDataFaceList::getSampleType -( - const octree& oc, - const point& sample -) const -{ - // Need to determine whether sample is 'inside' or 'outside' - // Done by finding nearest face. This gives back a face which is - // guaranteed to contain nearest point. This point can be - // - in interior of face: compare to face normal - // - on edge of face: compare to edge normal - // - on point of face: compare to point normal - // Unfortunately the octree does not give us back the intersection point - // or where on the face it has hit so we have to recreate all that - // information. - - - // Find nearest face to sample - treeBoundBox tightest(treeBoundBox::greatBox); - - scalar tightestDist = GREAT; - - label index = oc.findNearest(sample, tightest, tightestDist); - - if (index == -1) - { - FatalErrorIn - ( - "octreeDataFaceList::getSampleType" - "(octree&, const point&)" - ) << "Could not find " << sample << " in octree." - << abort(FatalError); - } - - label faceI = faceLabels_[index]; - - // Get actual intersection point on face - - if (debug & 2) - { - Pout<< "getSampleType : sample:" << sample - << " nearest face:" << faceI; - } - - const face& f = mesh_.localFaces()[faceI]; - - const pointField& points = mesh_.localPoints(); - - pointHit curHit = f.nearestPoint(sample, points); - - // - // 1] Check whether sample is above face - // - - if (curHit.hit()) - { - // Simple case. Compare to face normal. - - if (debug & 2) - { - Pout<< " -> face hit:" << curHit.hitPoint() - << " comparing to face normal " << mesh_.faceNormals()[faceI] - << endl; - } - return octree::getVolType - ( - mesh_.faceNormals()[faceI], - sample - curHit.hitPoint() - ); - } - - if (debug & 2) - { - Pout<< " -> face miss:" << curHit.missPoint(); - } - - // - // 2] Check whether intersection is on one of the face vertices or - // face centre - // - - // typical dimension as sqrt of face area. - scalar typDim = sqrt(mag(f.normal(points))) + VSMALL; - - forAll(f, fp) - { - if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol) - { - // Face intersection point equals face vertex fp - - if (debug & 2) - { - Pout<< " -> face point hit :" << points[f[fp]] - << " point normal:" << mesh_.pointNormals()[f[fp]] - << " distance:" - << mag(points[f[fp]] - curHit.missPoint())/typDim - << endl; - } - return octree::getVolType - ( - mesh_.pointNormals()[f[fp]], - sample - curHit.missPoint() - ); - } - } - - // Get face centre - point ctr(f.centre(points)); - - if ((mag(ctr - curHit.missPoint())/typDim) < tol) - { - // Face intersection point equals face centre. Normal at face centre - // is already average of face normals - - if (debug & 2) - { - Pout<< " -> centre hit:" << ctr - << " distance:" - << mag(ctr - curHit.missPoint())/typDim - << endl; - } - - return octree::getVolType - ( - mesh_.faceNormals()[faceI], - sample - curHit.missPoint() - ); - } - - - // - // 3] Get the 'real' edge the face intersection is on - // - - const labelList& myEdges = mesh_.faceEdges()[faceI]; - - forAll(myEdges, myEdgeI) - { - const edge& e = mesh_.edges()[myEdges[myEdgeI]]; - - pointHit edgeHit = - line - ( - points[e.start()], - points[e.end()] - ).nearestDist(sample); - - point edgePoint; - if (edgeHit.hit()) - { - edgePoint = edgeHit.hitPoint(); - } - else - { - edgePoint = edgeHit.missPoint(); - } - - - if ((mag(edgePoint - curHit.missPoint())/typDim) < tol) - { - // Face intersection point lies on edge e - - // Calculate edge normal (wrong: uses face normals instead of - // triangle normals) - const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]]; - - vector edgeNormal(vector::zero); - - forAll(myFaces, myFaceI) - { - edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]]; - } - - if (debug & 2) - { - Pout<< " -> real edge hit point:" << edgePoint - << " comparing to edge normal:" << edgeNormal - << endl; - } - - // Found face intersection point on this edge. Compare to edge - // normal - return octree::getVolType - ( - edgeNormal, - sample - curHit.missPoint() - ); - } - } - - - // - // 4] Get the internal edge (vertex - faceCentre) the face intersection - // is on - // - - forAll(f, fp) - { - pointHit edgeHit = - line - ( - points[f[fp]], - ctr - ).nearestDist(sample); - - point edgePoint; - if (edgeHit.hit()) - { - edgePoint = edgeHit.hitPoint(); - } - else - { - edgePoint = edgeHit.missPoint(); - } - - if ((mag(edgePoint - curHit.missPoint())/typDim) < tol) - { - // Face intersection point lies on edge between two face triangles - - // Calculate edge normal as average of the two triangle normals - label fpPrev = f.rcIndex(fp); - label fpNext = f.fcIndex(fp); - - vector e = points[f[fp]] - ctr; - vector ePrev = points[f[fpPrev]] - ctr; - vector eNext = points[f[fpNext]] - ctr; - - vector nLeft = ePrev ^ e; - nLeft /= mag(nLeft) + VSMALL; - - vector nRight = e ^ eNext; - nRight /= mag(nRight) + VSMALL; - - if (debug & 2) - { - Pout<< " -> internal edge hit point:" << edgePoint - << " comparing to edge normal " - << 0.5*(nLeft + nRight) - << endl; - } - - // Found face intersection point on this edge. Compare to edge - // normal - return octree::getVolType - ( - 0.5*(nLeft + nRight), - sample - curHit.missPoint() - ); - } - } - - if (debug & 2) - { - Pout<< "Did not find sample " << sample - << " anywhere related to nearest face " << faceI << endl - << "Face:"; - - forAll(f, fp) - { - Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]] - << endl; - } - } - - // Can't determine status of sample with respect to nearest face. - // Either - // - tolerances are wrong. (if e.g. face has zero area) - // - or (more likely) surface is not closed. - - return octree::UNKNOWN; -} - - -bool Foam::octreeDataFaceList::overlaps -( - const label index, - const treeBoundBox& sampleBb -) const -{ - return sampleBb.overlaps(allBb_[index]); -} - - -bool Foam::octreeDataFaceList::contains -( - const label, - const point& -) const -{ - notImplemented - ( - "octreeDataFaceList::contains(const label, const point&)" - ); - return false; -} - - -bool Foam::octreeDataFaceList::intersects -( - const label index, - const point& start, - const point& end, - point& intersectionPoint -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - const vector dir(end - start); - - // Disable picking up intersections behind us. - scalar oldTol = intersection::setPlanarTol(0.0); - - pointHit inter = - f.ray - ( - start, - dir, - mesh_.localPoints(), - intersection::HALF_RAY, - intersection::VECTOR - ); - - intersection::setPlanarTol(oldTol); - - if (inter.hit() && inter.distance() <= mag(dir)) - { - intersectionPoint = inter.hitPoint(); - - return true; - } - else - { - return false; - } -} - - -bool Foam::octreeDataFaceList::findTightest -( - const label index, - const point& sample, - treeBoundBox& tightest -) const -{ - // Get nearest and furthest away vertex - point myNear, myFar; - allBb_[index].calcExtremities(sample, myNear, myFar); - - const point dist = myFar - sample; - scalar myFarDist = mag(dist); - - point tightestNear, tightestFar; - tightest.calcExtremities(sample, tightestNear, tightestFar); - - scalar tightestFarDist = mag(tightestFar - sample); - - if (tightestFarDist < myFarDist) - { - // Keep current tightest. - return false; - } - else - { - // Construct bb around sample and myFar - const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z())); - - tightest.min() = sample - dist2; - tightest.max() = sample + dist2; - - return true; - } -} - - -// Determine numerical value of sign of sample compared to shape at index -Foam::scalar Foam::octreeDataFaceList::calcSign -( - const label index, - const point& sample, - vector& -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - point ctr = f.centre(mesh_.localPoints()); - - vector vec = sample - ctr; - - vec /= mag(vec) + VSMALL; - - return mesh_.faceNormals()[faceI] & vec; -} - - -// Calculate nearest point on/in shapei -Foam::scalar Foam::octreeDataFaceList::calcNearest -( - const label index, - const point& sample, - point& nearest -) const -{ - label faceI = faceLabels_[index]; - - const face& f = mesh_.localFaces()[faceI]; - - pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints()); - - if (nearHit.hit()) - { - nearest = nearHit.hitPoint(); - } - else - { - nearest = nearHit.missPoint(); - } - - if (debug & 1) - { - point ctr = f.centre(mesh_.localPoints()); - - scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest); - - Pout<< "octreeDataFaceList::calcNearest : " - << "sample:" << sample - << " faceI:" << faceI - << " ctr:" << ctr - << " sign:" << sign - << " nearest point:" << nearest - << " distance to face:" << nearHit.distance() - << endl; - } - return nearHit.distance(); -} - - -void Foam::octreeDataFaceList::write -( - Ostream& os, - const label index -) const -{ - os << faceLabels_[index] << " " << allBb_[index]; -} - - -// ************************************************************************* // diff --git a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H b/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H deleted file mode 100644 index addd71d79b..0000000000 --- a/src/dynamicMesh/boundaryMesh/octreeDataFaceList.H +++ /dev/null @@ -1,220 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::octreeDataFaceList - -Description - Holds data for octree to work on list of faces on a bMesh - (= PrimitivePatch which holds faces, not references them) - Same as octreeDataFace except for that. - -SourceFiles - octreeDataFaceList.C - -\*---------------------------------------------------------------------------*/ - -#ifndef octreeDataFaceList_H -#define octreeDataFaceList_H - -#include "treeBoundBoxList.H" -#include "faceList.H" -#include "point.H" -#include "className.H" -#include "bMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -template class octree; - -/*---------------------------------------------------------------------------*\ - Class octreeDataFaceList Declaration -\*---------------------------------------------------------------------------*/ - -class octreeDataFaceList -{ - // Static data - - //- tolerance on linear dimensions - static scalar tol; - - - // Static function - - static inline label nexti(label max, label i) - { - return (i + 1) % max; - } - - - // Private data - - //- the mesh - const bMesh& mesh_; - - //- labels (in mesh indexing) of faces - labelList faceLabels_; - - //- bbs for all above faces - treeBoundBoxList allBb_; - - - // Private Member Functions - - //- Set allBb to tight fitting bounding box - void calcBb(); - -public: - - // Declare name of the class and its debug switch - ClassName("octreeDataFaceList"); - - // Constructors - - //- Construct from all faces in bMesh. - octreeDataFaceList(const bMesh& mesh); - - //- Construct from selected faces in bMesh. - octreeDataFaceList(const bMesh& mesh, const labelList& faceLabels); - - //- Construct as copy - octreeDataFaceList(const octreeDataFaceList&); - - - // Destructor - - ~octreeDataFaceList(); - - - // Member Functions - - // Access - - const bMesh& mesh() const - { - return mesh_; - } - - const labelList& faceLabels() const - { - return faceLabels_; - } - - const treeBoundBoxList& allBb() const - { - return allBb_; - } - - label size() const - { - return allBb_.size(); - } - - // Search - - //- Get type of sample - label getSampleType - ( - const octree&, - const point& - ) const; - - //- Does (bb of) shape at index overlap bb - bool overlaps - ( - const label index, - const treeBoundBox& sampleBb - ) const; - - //- Does shape at index contain sample - bool contains - ( - const label index, - const point& sample - ) const; - - //- Segment (from start to end) intersection with shape - // at index. If intersects returns true and sets intersectionPoint - bool intersects - ( - const label index, - const point& start, - const point& end, - point& intersectionPoint - ) const; - - //- Sets newTightest to bounding box (and returns true) if - // nearer to sample than tightest bounding box. Otherwise - // returns false. - bool findTightest - ( - const label index, - const point& sample, - treeBoundBox& tightest - ) const; - - //- Given index get unit normal and calculate (numerical) sign - // of sample. - // Used to determine accuracy of calcNearest or inside/outside. - scalar calcSign - ( - const label index, - const point& sample, - vector& n - ) const; - - //- Calculates nearest (to sample) point in shape. - // Returns point and mag(nearest - sample). Returns GREAT if - // sample does not project onto (triangle decomposition) of face. - scalar calcNearest - ( - const label index, - const point& sample, - point& nearest - ) const; - - - // Edit - - // Write - - //- Write shape at index - void write(Ostream& os, const label index) const; -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - -#endif - -// ************************************************************************* // diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 272f3e067c..ada844b5a5 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -50,6 +50,7 @@ indexedOctree/treeDataCell.C indexedOctree/treeDataEdge.C indexedOctree/treeDataFace.C indexedOctree/treeDataPoint.C +indexedOctree/treeDataPrimitivePatchName.C indexedOctree/treeDataTriSurface.C searchableSurface = searchableSurface diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C new file mode 100644 index 0000000000..97b7821b1e --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C @@ -0,0 +1,567 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "treeDataPrimitivePatch.H" +#include "indexedOctree.H" +#include "triangleFuncs.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::scalar +Foam::treeDataPrimitivePatch:: +tolSqr = sqr(1E-6); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::treeBoundBox +Foam::treeDataPrimitivePatch:: +calcBb +( + const pointField& points, + const face& f +) +{ + treeBoundBox bb(points[f[0]], points[f[0]]); + + for (label fp = 1; fp < f.size(); fp++) + { + const point& p = points[f[fp]]; + + bb.min() = min(bb.min(), p); + bb.max() = max(bb.max(), p); + } + return bb; +} + + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +void Foam::treeDataPrimitivePatch:: +update() +{ + if (cacheBb_) + { + bbs_.setSize(patch_.size()); + + forAll(patch_, i) + { + bbs_[i] = calcBb(patch_.points(), patch_[i]); + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::treeDataPrimitivePatch:: +treeDataPrimitivePatch +( + const bool cacheBb, + const PrimitivePatch& patch +) +: + patch_(patch), + cacheBb_(cacheBb) +{ + update(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::pointField +Foam::treeDataPrimitivePatch:: +points() const +{ + pointField cc(patch_.size()); + + forAll(patch_, i) + { + cc[i] = patch_[i].centre(patch_.points()); + } + + return cc; +} + + +//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. +// Only makes sense for closed surfaces. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +Foam::label +Foam::treeDataPrimitivePatch:: +getVolumeType +( + const indexedOctree + < + treeDataPrimitivePatch + < + Face, + FaceList, + PointField, + PointType + > + >& oc, + const point& sample +) const +{ + // Need to determine whether sample is 'inside' or 'outside' + // Done by finding nearest face. This gives back a face which is + // guaranteed to contain nearest point. This point can be + // - in interior of face: compare to face normal + // - on edge of face: compare to edge normal + // - on point of face: compare to point normal + // Unfortunately the octree does not give us back the intersection point + // or where on the face it has hit so we have to recreate all that + // information. + + + // Find nearest face to sample + pointIndexHit info = oc.findNearest(sample, sqr(GREAT)); + + if (info.index() == -1) + { + FatalErrorIn + ( + "treeDataPrimitivePatch::getSampleType" + "(indexedOctree&, const point&)" + ) << "Could not find " << sample << " in octree." + << abort(FatalError); + } + + + // Get actual intersection point on face + label faceI = info.index(); + + if (debug & 2) + { + Pout<< "getSampleType : sample:" << sample + << " nearest face:" << faceI; + } + + const pointField& points = patch_.localPoints(); + const face& f = patch_.localFaces()[faceI]; + + // Retest to classify where on face info is. Note: could be improved. We + // already have point. + + pointHit curHit = f.nearestPoint(sample, points); + const vector area = f.normal(points); + const point& curPt = curHit.rawPoint(); + + // + // 1] Check whether sample is above face + // + + if (curHit.hit()) + { + // Nearest point inside face. Compare to face normal. + + if (debug & 2) + { + Pout<< " -> face hit:" << curPt + << " comparing to face normal " << area << endl; + } + return indexedOctree::getSide + ( + area, + sample - curPt + ); + } + + if (debug & 2) + { + Pout<< " -> face miss:" << curPt; + } + + // + // 2] Check whether intersection is on one of the face vertices or + // face centre + // + + const scalar typDimSqr = mag(area) + VSMALL; + + forAll(f, fp) + { + if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point equals face vertex fp + + // Calculate point normal (wrong: uses face normals instead of + // triangle normals) + + return indexedOctree::getSide + ( + patch_.pointNormals()[f[fp]], + sample - curPt + ); + } + } + + const point fc(f.centre(points)); + + if ((magSqr(fc - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point equals face centre. Normal at face centre + // is already average of face normals + + if (debug & 2) + { + Pout<< " -> centre hit:" << fc + << " distance:" << magSqr(fc - curPt)/typDimSqr << endl; + } + + return indexedOctree::getSide + ( + area, + sample - curPt + ); + } + + + + // + // 3] Get the 'real' edge the face intersection is on + // + + const labelList& fEdges = patch_.faceEdges()[faceI]; + + forAll(fEdges, fEdgeI) + { + label edgeI = fEdges[fEdgeI]; + const edge& e = patch_.edges()[edgeI]; + + pointHit edgeHit = e.line(points).nearestDist(sample); + + if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point lies on edge e + + // Calculate edge normal (wrong: uses face normals instead of + // triangle normals) + const labelList& eFaces = patch_.edgeFaces()[edgeI]; + + vector edgeNormal(vector::zero); + + forAll(eFaces, i) + { + edgeNormal += patch_.faceNormal()[eFaces[i]]; + } + + if (debug & 2) + { + Pout<< " -> real edge hit point:" << edgeHit.rawPoint() + << " comparing to edge normal:" << edgeNormal + << endl; + } + + // Found face intersection point on this edge. Compare to edge + // normal + return indexedOctree::getSide + ( + edgeNormal, + sample - curPt + ); + } + } + + + // + // 4] Get the internal edge the face intersection is on + // + + forAll(f, fp) + { + pointHit edgeHit = linePointRef + ( + points[f[fp]], + fc + ).nearestDist(sample); + + if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr) + { + // Face intersection point lies on edge between two face triangles + + // Calculate edge normal as average of the two triangle normals + vector e = points[f[fp]] - fc; + vector ePrev = points[f[f.rcIndex(fp)]] - fc; + vector eNext = points[f[f.fcIndex(fp)]] - fc; + + vector nLeft = ePrev ^ e; + nLeft /= mag(nLeft) + VSMALL; + + vector nRight = e ^ eNext; + nRight /= mag(nRight) + VSMALL; + + if (debug & 2) + { + Pout<< " -> internal edge hit point:" << edgeHit.rawPoint() + << " comparing to edge normal " + << 0.5*(nLeft + nRight) + << endl; + } + + // Found face intersection point on this edge. Compare to edge + // normal + return indexedOctree::getSide + ( + 0.5*(nLeft + nRight), + sample - curPt + ); + } + } + + if (debug & 2) + { + Pout<< "Did not find sample " << sample + << " anywhere related to nearest face " << faceI << endl + << "Face:"; + + forAll(f, fp) + { + Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]] + << endl; + } + } + + // Can't determine status of sample with respect to nearest face. + // Either + // - tolerances are wrong. (if e.g. face has zero area) + // - or (more likely) surface is not closed. + + return indexedOctree::UNKNOWN; +} + + +// Check if any point on shape is inside cubeBb. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +bool +Foam::treeDataPrimitivePatch:: +overlaps +( + const label index, + const treeBoundBox& cubeBb +) const +{ + // 1. Quick rejection: bb does not intersect face bb at all + if (cacheBb_) + { + if (!cubeBb.overlaps(bbs_[index])) + { + return false; + } + } + else + { + if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index]))) + { + return false; + } + } + + + // 2. Check if one or more face points inside + + const pointField& points = patch_.points(); + const face& f = patch_[index]; + + forAll(f, fp) + { + if (cubeBb.contains(points[f[fp]])) + { + return true; + } + } + + // 3. Difficult case: all points are outside but connecting edges might + // go through cube. Use triangle-bounding box intersection. + const point fc = f.centre(points); + + forAll(f, fp) + { + bool triIntersects = triangleFuncs::intersectBb + ( + points[f[fp]], + points[f[f.fcIndex(fp)]], + fc, + cubeBb + ); + + if (triIntersects) + { + return true; + } + } + return false; +} + + +// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex, +// nearestPoint. +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +void +Foam::treeDataPrimitivePatch:: +findNearest +( + const labelList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& minIndex, + point& nearestPoint +) const +{ + const pointField& points = patch_.points(); + + forAll(indices, i) + { + label index = indices[i]; + + const face& f = patch_[index]; + + pointHit nearHit = f.nearestPoint(sample, points); + scalar distSqr = sqr(nearHit.distance()); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + minIndex = index; + nearestPoint = nearHit.rawPoint(); + } + } +} + + +template +< + class Face, + template class FaceList, + class PointField, + class PointType +> +bool +Foam::treeDataPrimitivePatch:: +intersects +( + const label index, + const point& start, + const point& end, + point& intersectionPoint +) const +{ + // Do quick rejection test + if (cacheBb_) + { + const treeBoundBox& faceBb = bbs_[index]; + + if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0) + { + // start and end in same block outside of faceBb. + return false; + } + } + + const pointField& points = patch_.points(); + const face& f = patch_[index]; + const point fc = f.centre(points); + const vector dir(end - start); + + pointHit inter = patch_[index].intersection + ( + start, + dir, + fc, + points, + intersection::HALF_RAY + ); + + if (inter.hit() && inter.distance() <= 1) + { + // Note: no extra test on whether intersection is in front of us + // since using half_ray + intersectionPoint = inter.hitPoint(); + return true; + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.H b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H new file mode 100644 index 0000000000..9fd4da2fd5 --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.H @@ -0,0 +1,210 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::treeDataPrimitivePatch + +Description + Encapsulation of data needed to search on PrimitivePatches + +SourceFiles + treeDataPrimitivePatch.C + +\*---------------------------------------------------------------------------*/ + +#ifndef treeDataPrimitivePatch_H +#define treeDataPrimitivePatch_H + +#include "PrimitivePatch.H" +//#include "indexedOctree.H" +#include "treeBoundBoxList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +template class indexedOctree; + + +/*---------------------------------------------------------------------------*\ + Class treeDataPrimitivePatchName Declaration +\*---------------------------------------------------------------------------*/ + +TemplateName(treeDataPrimitivePatch); + + +/*---------------------------------------------------------------------------*\ + Class treeDataPrimitivePatch Declaration +\*---------------------------------------------------------------------------*/ + +template +< + class Face, + template class FaceList, + class PointField, + class PointType=point +> +class treeDataPrimitivePatch +: + public treeDataPrimitivePatchName +{ + // Static data + + //- tolerance on linear dimensions + static scalar tolSqr; + + // Private data + + //- Underlying geometry + const PrimitivePatch& patch_; + + //- Whether to precalculate and store face bounding box + const bool cacheBb_; + + //- face bounding boxes (valid only if cacheBb_) + treeBoundBoxList bbs_; + + + // Private Member Functions + + //- Calculate face bounding box + static treeBoundBox calcBb(const pointField&, const face&); + + //- Initialise all member data + void update(); + +public: + + // Constructors + + //- Construct from patch. + treeDataPrimitivePatch + ( + const bool cacheBb, + const PrimitivePatch& + ); + + + // Member Functions + + // Access + + label size() const + { + return patch_.size(); + } + + //- Get representative point cloud for all shapes inside + // (one point per shape) + pointField points() const; + + + // Search + + //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface. + // Only makes sense for closed surfaces. + label getVolumeType + ( + const indexedOctree + < + treeDataPrimitivePatch + < + Face, + FaceList, + PointField, + PointType + > + >&, + const point& + ) const; + + //- Does (bb of) shape at index overlap bb + bool overlaps + ( + const label index, + const treeBoundBox& sampleBb + ) const; + + //- Calculates nearest (to sample) point in shape. + // Returns actual point and distance (squared) + void findNearest + ( + const labelList& indices, + const point& sample, + + scalar& nearestDistSqr, + label& nearestIndex, + point& nearestPoint + ) const; + + //- Calculates nearest (to line) point in shape. + // Returns point and distance (squared) + void findNearest + ( + const labelList& indices, + const linePointRef& ln, + + treeBoundBox& tightest, + label& minIndex, + point& linePoint, + point& nearestPoint + ) const + { + notImplemented + ( + "treeDataPrimitivePatch::findNearest" + "(const labelList&, const linePointRef&, ..)" + ); + } + + //- Calculate intersection of shape with ray. Sets result + // accordingly + bool intersects + ( + const label index, + const point& start, + const point& end, + point& result + ) const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "treeDataPrimitivePatch.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C b/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C new file mode 100644 index 0000000000..59cf44b8a6 --- /dev/null +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatchName.C @@ -0,0 +1,33 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "treeDataPrimitivePatch.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::treeDataPrimitivePatchName, 0); + +// ************************************************************************* // From 100f9c25456a94d2f5435033279bf21e798a0a44 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 11 Nov 2009 19:14:33 +0000 Subject: [PATCH 05/29] stitching points if no points have been merged --- src/triSurface/triSurface/stitchTriangles.C | 8 +++++--- src/triSurface/triSurface/triSurface.C | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/triSurface/triSurface/stitchTriangles.C b/src/triSurface/triSurface/stitchTriangles.C index 463eb41ba9..7a441ac679 100644 --- a/src/triSurface/triSurface/stitchTriangles.C +++ b/src/triSurface/triSurface/stitchTriangles.C @@ -46,10 +46,7 @@ bool triSurface::stitchTriangles pointField newPoints; bool hasMerged = mergePoints(rawPoints, tol, verbose, pointMap, newPoints); - pointField& ps = storedPoints(); - // Set the coordinates to the merged ones - ps = newPoints; if (hasMerged) { @@ -59,6 +56,11 @@ bool triSurface::stitchTriangles << " points down to " << newPoints.size() << endl; } + pointField& ps = storedPoints(); + + // Set the coordinates to the merged ones + ps = newPoints; + // Reset the triangle point labels to the unique points array label newTriangleI = 0; forAll(*this, i) diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C index 65258557e8..83674b6b90 100644 --- a/src/triSurface/triSurface/triSurface.C +++ b/src/triSurface/triSurface/triSurface.C @@ -814,7 +814,7 @@ void Foam::triSurface::scalePoints(const scalar scaleFactor) void Foam::triSurface::cleanup(const bool verbose) { // Merge points (already done for STL, TRI) - stitchTriangles(pointField(points()), SMALL, verbose); + stitchTriangles(points(), SMALL, verbose); // Merging points might have changed geometric factors clearOut(); From a54bbffbbfc612034fa9b66539fdd7f37042f93a Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 12 Nov 2009 13:28:18 +0000 Subject: [PATCH 06/29] typo --- .../alphatJayatillekeWallFunctionFvPatchScalarField.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C index b5022e2756..5c73b8b146 100644 --- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C +++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C @@ -239,7 +239,7 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs() // Molecular Prandtl number scalar Pr = muw[faceI]/alphaw[faceI]; - // Molecular-to-turbulenbt Prandtl number ratio + // Molecular-to-turbulent Prandtl number ratio scalar Prat = Pr/Prt_; // Thermal sublayer thickness From 963a8cb8d52db5da8a5484c9c5b3a581297b20ec Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 12 Nov 2009 13:28:55 +0000 Subject: [PATCH 07/29] cleaner code --- .../readAdditionalSolutionControls.H | 21 +++---------------- 1 file changed, 3 insertions(+), 18 deletions(-) diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H index 8159205cae..340623fbfd 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H @@ -1,20 +1,5 @@ dictionary additional = mesh.solutionDict().subDict("additional"); -bool dpdt = true; -if (additional.found("dpdt")) -{ - additional.lookup("dpdt") >> dpdt; -} - -bool eWork = true; -if (additional.found("eWork")) -{ - additional.lookup("eWork") >> eWork; -} - -bool hWork = true; -if (additional.found("hWork")) -{ - additional.lookup("hWork") >> hWork; -} - +bool dpdt = additional.lookupOrDefault("dpdt", true); +bool eWork = additional.lookupOrDefault("eWork", true); +bool hWork = additional.lookupOrDefault("hWork", true); From 9e57cef139841dab800bc68b2ace7619840bf3b1 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 12 Nov 2009 18:23:22 +0000 Subject: [PATCH 08/29] fixed reverseFaceMap --- src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C index 65965f15d9..7e38d09831 100644 --- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C +++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C @@ -146,6 +146,10 @@ void Foam::singleCellFvMesh::agglomerateMesh // From new patch face back to agglomeration patchFaceMap_.setSize(oldPatches.size()); + // From fine face to coarse face (or -1) + reverseFaceMap_.setSize(mesh.nFaces()); + reverseFaceMap_.labelList::operator=(-1); + // Face counter coarseI = 0; From 6e793380d9697cb4cc1644f64ab474028e3854af Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 12 Nov 2009 18:25:08 +0000 Subject: [PATCH 09/29] rmeoved printing --- src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C index 3df03c8c70..37562abc35 100644 --- a/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C +++ b/src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C @@ -330,8 +330,6 @@ void Foam::fvMeshAdder::MapVolFields if (fieldsToAdd.found(fld.name())) { - Pout<< "Mapping field " << fld.name() << endl; - const GeometricField& fldToAdd = *fieldsToAdd[fld.name()]; @@ -339,7 +337,7 @@ void Foam::fvMeshAdder::MapVolFields } else { - WarningIn("fvMeshAdder::MapVolFields") + WarningIn("fvMeshAdder::MapVolFields(..)") << "Not mapping field " << fld.name() << " since not present on mesh to add" << endl; @@ -642,15 +640,13 @@ void Foam::fvMeshAdder::MapSurfaceFields if (fieldsToAdd.found(fld.name())) { - Pout<< "Mapping field " << fld.name() << endl; - const fldType& fldToAdd = *fieldsToAdd[fld.name()]; MapSurfaceField(meshMap, fld, fldToAdd); } else { - WarningIn("fvMeshAdder::MapSurfaceFields") + WarningIn("fvMeshAdder::MapSurfaceFields(..)") << "Not mapping field " << fld.name() << " since not present on mesh to add" << endl; From 1d4533a41eb2f95e4a11a6d704f1be3e1352a0a6 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 12 Nov 2009 18:25:37 +0000 Subject: [PATCH 10/29] >80 line length --- .../Lists/CompactListList/CompactListList.H | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H b/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H index f1875dde70..b2aad34d22 100644 --- a/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H +++ b/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H @@ -61,8 +61,16 @@ namespace Foam template class CompactListList; -template Istream& operator>>(Istream&, CompactListList&); -template Ostream& operator<<(Ostream&, const CompactListList&); +template Istream& operator>> +( + Istream&, + CompactListList& +); +template Ostream& operator<< +( + Ostream&, + const CompactListList& +); /*---------------------------------------------------------------------------*\ From b7845edeba957575b3ed69462f93db1ce123a4cb Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 13 Nov 2009 16:09:27 +0000 Subject: [PATCH 11/29] updated to new thermo --- applications/test/readCHEMKINIII/readCHEMKINIII.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/test/readCHEMKINIII/readCHEMKINIII.C b/applications/test/readCHEMKINIII/readCHEMKINIII.C index a028e9fd2b..50721d64ea 100644 --- a/applications/test/readCHEMKINIII/readCHEMKINIII.C +++ b/applications/test/readCHEMKINIII/readCHEMKINIII.C @@ -57,7 +57,7 @@ int main(int argc, char *argv[]) //Info<< ck.specieThermo() << endl; //Info<< ck.reactions() << endl; - PtrList reactions = ck.reactions(); + const SLPtrList& reactions = ck.reactions(); { OFstream reactionStream("reactions"); @@ -70,17 +70,17 @@ int main(int argc, char *argv[]) label nReactions(readLabel(reactionStream)); reactionStream.readBeginList(args.executable().c_str()); - PtrList testReactions(nReactions); + PtrList testReactions(nReactions); forAll (testReactions, i) { testReactions.set ( i, - chemkinReader::reaction::New + gasReaction::New ( ck.species(), - ck.specieThermo(), + ck.speciesThermo(), reactionStream ) ); From 8349d891109d2e4c1c5e8165678d82e18fa96c6d Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 13 Nov 2009 16:09:39 +0000 Subject: [PATCH 12/29] construct from dictionary --- .../GeometricField/GeometricField.C | 77 +++++++++++++++---- .../GeometricField/GeometricField.H | 11 +++ 2 files changed, 75 insertions(+), 13 deletions(-) diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C index 705ee5e2b4..25a90ba64c 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C @@ -52,20 +52,11 @@ Foam::tmp typename Foam::GeometricField:: GeometricBoundaryField > -Foam::GeometricField::readField(Istream& is) +Foam::GeometricField::readField +( + const dictionary& fieldDict +) { - if (is.version() < 2.0) - { - FatalIOErrorIn - ( - "GeometricField::readField(Istream&)", - is - ) << "IO versions < 2.0 are not supported." - << exit(FatalIOError); - } - - dictionary fieldDict(is); - DimensionedField::readField(fieldDict, "internalField"); tmp tboundaryField @@ -96,6 +87,28 @@ Foam::GeometricField::readField(Istream& is) } +template class PatchField, class GeoMesh> +Foam::tmp +< + typename Foam::GeometricField:: + GeometricBoundaryField +> +Foam::GeometricField::readField(Istream& is) +{ + if (is.version() < 2.0) + { + FatalIOErrorIn + ( + "GeometricField::readField(Istream&)", + is + ) << "IO versions < 2.0 are not supported." + << exit(FatalIOError); + } + + return readField(dictionary(is)); +} + + template class PatchField, class GeoMesh> bool Foam::GeometricField::readIfPresent() { @@ -404,6 +417,44 @@ Foam::GeometricField::GeometricField } +template class PatchField, class GeoMesh> +Foam::GeometricField::GeometricField +( + const IOobject& io, + const Mesh& mesh, + const dictionary& dict +) +: + DimensionedField(io, mesh, dimless), + timeIndex_(this->time().timeIndex()), + field0Ptr_(NULL), + fieldPrevIterPtr_(NULL), + boundaryField_(*this, readField(dict)) +{ + // Check compatibility between field and mesh + + if (this->size() != GeoMesh::size(this->mesh())) + { + FatalErrorIn + ( + "GeometricField::GeometricField" + "(const IOobject&, const Mesh&, const dictionary&)" + ) << " number of field elements = " << this->size() + << " number of mesh elements = " << GeoMesh::size(this->mesh()) + << exit(FatalIOError); + } + + readOldTimeIfPresent(); + + if (debug) + { + Info<< "Finishing dictionary-construct of " + "GeometricField" + << endl << this->info() << endl; + } +} + + // construct as copy template class PatchField, class GeoMesh> Foam::GeometricField::GeometricField diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H index 71e2bbf34c..1abca4cc87 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H +++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H @@ -240,6 +240,9 @@ private: // Private member functions + //- Read the field from the dictionary + tmp readField(const dictionary&); + //- Read the field from the given stream tmp readField(Istream& is); @@ -327,6 +330,14 @@ public: Istream& ); + //- Construct from dictionary + GeometricField + ( + const IOobject&, + const Mesh&, + const dictionary& + ); + //- Construct as copy GeometricField ( From 03d777f594b1092c6d773e8531d38c1cbc946c26 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 13 Nov 2009 16:10:11 +0000 Subject: [PATCH 13/29] message size not set; debug printing --- src/Pstream/mpi/UIPread.C | 57 +++++++++++++++++++++++++++++++++++++- src/Pstream/mpi/UOPwrite.C | 36 +++++++++++++++++++++++- src/Pstream/mpi/UPstream.C | 45 ++++++++++++++++++++++++++++++ 3 files changed, 136 insertions(+), 2 deletions(-) diff --git a/src/Pstream/mpi/UIPread.C b/src/Pstream/mpi/UIPread.C index 14360e0cc0..fb6a9ec7d7 100644 --- a/src/Pstream/mpi/UIPread.C +++ b/src/Pstream/mpi/UIPread.C @@ -66,6 +66,14 @@ Foam::UIPstream::UIPstream label wantedSize = externalBuf_.capacity(); + if (debug) + { + Pout<< "UIPstream::UIPstream : read from:" << fromProcNo + << " tag:" << tag << " wanted size:" << wantedSize + << Foam::endl; + } + + // If the buffer size is not specified, probe the incomming message // and set it if (!wantedSize) @@ -75,6 +83,12 @@ Foam::UIPstream::UIPstream externalBuf_.setCapacity(messageSize_); wantedSize = messageSize_; + + if (debug) + { + Pout<< "UIPstream::UIPstream : probed size:" << wantedSize + << Foam::endl; + } } messageSize_ = UIPstream::read @@ -127,6 +141,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers) if (commsType() == UPstream::nonBlocking) { // Message is already received into externalBuf + messageSize_ = buffers.recvBuf_[fromProcNo].size(); } else { @@ -134,6 +149,14 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers) label wantedSize = externalBuf_.capacity(); + if (debug) + { + Pout<< "UIPstream::UIPstream PstreamBuffers :" + << " read from:" << fromProcNo + << " tag:" << tag_ << " wanted size:" << wantedSize + << Foam::endl; + } + // If the buffer size is not specified, probe the incomming message // and set it if (!wantedSize) @@ -143,6 +166,12 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers) externalBuf_.setCapacity(messageSize_); wantedSize = messageSize_; + + if (debug) + { + Pout<< "UIPstream::UIPstream PstreamBuffers : probed size:" + << wantedSize << Foam::endl; + } } messageSize_ = UIPstream::read @@ -180,6 +209,14 @@ Foam::label Foam::UIPstream::read const int tag ) { + if (debug) + { + Pout<< "UIPstream::read : starting read from:" << fromProcNo + << " tag:" << tag << " wanted size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << Foam::endl; + } + if (commsType == blocking || commsType == scheduled) { MPI_Status status; @@ -214,6 +251,14 @@ Foam::label Foam::UIPstream::read label messageSize; MPI_Get_count(&status, MPI_BYTE, &messageSize); + if (debug) + { + Pout<< "UIPstream::read : finished read from:" << fromProcNo + << " tag:" << tag << " read size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << Foam::endl; + } + if (messageSize > bufSize) { FatalErrorIn @@ -256,6 +301,15 @@ Foam::label Foam::UIPstream::read return 0; } + if (debug) + { + Pout<< "UIPstream::read : started read from:" << fromProcNo + << " tag:" << tag << " read size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << " request:" << PstreamGlobals::outstandingRequests_.size() + << Foam::endl; + } + PstreamGlobals::outstandingRequests_.append(request); // Assume the message is completely received. @@ -267,7 +321,8 @@ Foam::label Foam::UIPstream::read ( "UIPstream::read" "(const int fromProcNo, char* buf, std::streamsize bufSize)" - ) << "Unsupported communications type " << commsType + ) << "Unsupported communications type " + << commsType << Foam::abort(FatalError); return 0; diff --git a/src/Pstream/mpi/UOPwrite.C b/src/Pstream/mpi/UOPwrite.C index f8439ba517..9f8726d9a0 100644 --- a/src/Pstream/mpi/UOPwrite.C +++ b/src/Pstream/mpi/UOPwrite.C @@ -43,6 +43,14 @@ bool Foam::UOPstream::write const int tag ) { + if (debug) + { + Pout<< "UIPstream::write : starting write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << Foam::endl; + } + bool transferFailed = true; if (commsType == blocking) @@ -56,6 +64,14 @@ bool Foam::UOPstream::write tag, MPI_COMM_WORLD ); + + if (debug) + { + Pout<< "UIPstream::write : finished write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << Foam::endl; + } } else if (commsType == scheduled) { @@ -68,6 +84,14 @@ bool Foam::UOPstream::write tag, MPI_COMM_WORLD ); + + if (debug) + { + Pout<< "UIPstream::write : finished write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << Foam::endl; + } } else if (commsType == nonBlocking) { @@ -84,6 +108,15 @@ bool Foam::UOPstream::write &request ); + if (debug) + { + Pout<< "UIPstream::write : started write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << UPstream::commsTypeNames[commsType] + << " request:" << PstreamGlobals::outstandingRequests_.size() + << Foam::endl; + } + PstreamGlobals::outstandingRequests_.append(request); } else @@ -93,7 +126,8 @@ bool Foam::UOPstream::write "UOPstream::write" "(const int fromProcNo, char* buf, std::streamsize bufSize" ", const int)" - ) << "Unsupported communications type " << commsType + ) << "Unsupported communications type " + << UPstream::commsTypeNames[commsType] << Foam::abort(FatalError); } diff --git a/src/Pstream/mpi/UPstream.C b/src/Pstream/mpi/UPstream.C index b5a76e36e9..3623f5da58 100644 --- a/src/Pstream/mpi/UPstream.C +++ b/src/Pstream/mpi/UPstream.C @@ -70,6 +70,12 @@ bool Foam::UPstream::init(int& argc, char**& argv) MPI_Comm_size(MPI_COMM_WORLD, &numprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myProcNo_); + if (debug) + { + Pout<< "UPstream::init : initialised with numProcs:" << numprocs + << " myProcNo:" << myProcNo_ << endl; + } + if (numprocs <= 1) { FatalErrorIn("UPstream::init(int& argc, char**& argv)") @@ -124,6 +130,11 @@ bool Foam::UPstream::init(int& argc, char**& argv) void Foam::UPstream::exit(int errnum) { + if (debug) + { + Pout<< "UPstream::exit." << endl; + } + # ifndef SGIMPI int size; char* buff; @@ -164,6 +175,11 @@ void Foam::UPstream::abort() void Foam::reduce(scalar& Value, const sumOp& bop) { + if (Pstream::debug) + { + Pout<< "UPstream::reduce : value:" << Value << endl; + } + if (!UPstream::parRun()) { return; @@ -433,11 +449,23 @@ void Foam::reduce(scalar& Value, const sumOp& bop) } */ } + + if (Pstream::debug) + { + Pout<< "UPstream::reduce : reduced value:" << Value << endl; + } } void Foam::UPstream::waitRequests() { + if (debug) + { + Pout<< "UPstream::waitRequests : starting wait for " + << PstreamGlobals::outstandingRequests_.size() + << " outstanding requests." << endl; + } + if (PstreamGlobals::outstandingRequests_.size()) { if @@ -458,11 +486,22 @@ void Foam::UPstream::waitRequests() PstreamGlobals::outstandingRequests_.clear(); } + + if (debug) + { + Pout<< "UPstream::waitRequests : finished wait." << endl; + } } bool Foam::UPstream::finishedRequest(const label i) { + if (debug) + { + Pout<< "UPstream::waitRequests : starting wait for request:" << i + << endl; + } + if (i >= PstreamGlobals::outstandingRequests_.size()) { FatalErrorIn @@ -483,6 +522,12 @@ bool Foam::UPstream::finishedRequest(const label i) MPI_STATUS_IGNORE ); + if (debug) + { + Pout<< "UPstream::waitRequests : finished wait for request:" << i + << endl; + } + return flag != 0; } From 4a5f0fbebec395d92c190911b30d0bbf462de1bc Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 13 Nov 2009 16:10:25 +0000 Subject: [PATCH 14/29] extraneous pointer cast --- .../constraint/jumpCyclic/jumpCyclicFvPatchField.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C index 4f4ec60dd0..dcd3440e36 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C @@ -154,7 +154,7 @@ void jumpCyclicFvPatchField::updateInterfaceMatrix label sizeby2 = this->size()/2; const unallocLabelList& faceCells = this->cyclicPatch().faceCells(); - if (long(&psiInternal) == long(&this->internalField())) + if (&psiInternal == &this->internalField()) { tmp > tjf = jump(); const Field& jf = tjf(); From 4f1305316da6b1bc39671311935861252ac8e950 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 13 Nov 2009 16:10:54 +0000 Subject: [PATCH 15/29] fix multiple field transfer --- .../fvMeshDistribute/fvMeshDistribute.C | 549 ++++++++++-------- .../fvMeshDistribute/fvMeshDistribute.H | 44 +- .../fvMeshDistributeTemplates.C | 60 +- 3 files changed, 352 insertions(+), 301 deletions(-) diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C index fea10f410b..5d4144a32c 100644 --- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C +++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C @@ -25,8 +25,6 @@ License \*----------------------------------------------------------------------------*/ #include "fvMeshDistribute.H" -#include "ProcessorTopology.H" -#include "commSchedule.H" #include "PstreamCombineReduceOps.H" #include "fvMeshAdder.H" #include "faceCoupleInfo.H" @@ -39,6 +37,8 @@ License #include "mergePoints.H" #include "mapDistributePolyMesh.H" #include "surfaceFields.H" +#include "syncTools.H" +#include "CompactListList.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -47,71 +47,6 @@ defineTypeNameAndDebug(Foam::fvMeshDistribute, 0); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -//Foam::List Foam::fvMeshDistribute::getSchedule -//( -// const labelList& distribution -//) -//{ -// labelList nCellsPerProc(countCells(distribution)); -// -// if (debug) -// { -// Pout<< "getSchedule : Wanted distribution:" << nCellsPerProc << endl; -// } -// -// // Processors I need to send data to -// labelListList mySendProcs(Pstream::nProcs()); -// -// // Count -// label nSendProcs = 0; -// forAll(nCellsPerProc, sendProcI) -// { -// if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0) -// { -// nSendProcs++; -// } -// } -// -// // Fill -// mySendProcs[Pstream::myProcNo()].setSize(nSendProcs); -// nSendProcs = 0; -// forAll(nCellsPerProc, sendProcI) -// { -// if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0) -// { -// mySendProcs[Pstream::myProcNo()][nSendProcs++] = sendProcI; -// } -// } -// -// // Synchronise -// Pstream::gatherList(mySendProcs); -// Pstream::scatterList(mySendProcs); -// -// // Combine into list (same on all procs) giving sending and receiving -// // processor -// label nComms = 0; -// forAll(mySendProcs, procI) -// { -// nComms += mySendProcs[procI].size(); -// } -// -// List schedule(nComms); -// nComms = 0; -// -// forAll(mySendProcs, procI) -// { -// const labelList& sendProcs = mySendProcs[procI]; -// -// forAll(sendProcs, i) -// { -// schedule[nComms++] = labelPair(procI, sendProcs[i]); -// } -// } -// -// return schedule; -//} - - Foam::labelList Foam::fvMeshDistribute::select ( const bool selectEqual, @@ -144,14 +79,29 @@ Foam::labelList Foam::fvMeshDistribute::select // Check all procs have same names and in exactly same order. -void Foam::fvMeshDistribute::checkEqualWordList(const wordList& lst) +void Foam::fvMeshDistribute::checkEqualWordList +( + const string& msg, + const wordList& lst +) { - wordList myObjects(lst); + List allNames(Pstream::nProcs()); + allNames[Pstream::myProcNo()] = lst; + Pstream::gatherList(allNames); + Pstream::scatterList(allNames); - // Check that all meshes have same objects - Pstream::listCombineGather(myObjects, checkEqualType()); - // Below scatter only needed to balance sends and receives. - Pstream::listCombineScatter(myObjects); + for (label procI = 1; procI < Pstream::nProcs(); procI++) + { + if (allNames[procI] != allNames[0]) + { + FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)") + << "When checking for equal " << msg.c_str() << " :" << endl + << "processor0 has:" << allNames[0] << endl + << "processor" << procI << " has:" << allNames[procI] << endl + << msg.c_str() << " need to be synchronised on all processors." + << exit(FatalError); + } + } } @@ -664,21 +614,6 @@ Foam::autoPtr Foam::fvMeshDistribute::mergeSharedPoints << " newPointI:" << newPointI << abort(FatalError); } } - //- old: use pointToMaster map. - //forAll(constructMap, i) - //{ - // label oldPointI = constructMap[i]; - // - // // See if merged into other point - // Map