diff --git a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options index f1887d9aff..8ffc7f6d58 100644 --- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options +++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options @@ -5,11 +5,16 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/meshTools/lnInclude + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/fieldSources/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + EXE_LIBS = \ -lincompressibleTurbulenceModel \ -lincompressibleRASModels \ -lincompressibleTransportModels \ -lfiniteVolume \ - -lmeshTools + -lmeshTools \ + -lfieldSources \ + -lsampling diff --git a/applications/test/findSphereFeatureEdges-octree/Make/files b/applications/test/findSphereFeatureEdges-octree/Make/files new file mode 100644 index 0000000000..b9c5a5027f --- /dev/null +++ b/applications/test/findSphereFeatureEdges-octree/Make/files @@ -0,0 +1,3 @@ +Test-findSphereFeatureEdges-octree.C + +EXE = $(FOAM_USER_APPBIN)/Test-findSphereFeatureEdges-octree diff --git a/applications/test/findSphereFeatureEdges-octree/Make/options b/applications/test/findSphereFeatureEdges-octree/Make/options new file mode 100644 index 0000000000..6069034514 --- /dev/null +++ b/applications/test/findSphereFeatureEdges-octree/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/edgeMesh/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -ledgeMesh diff --git a/applications/test/findSphereFeatureEdges-octree/Test-findSphereFeatureEdges-octree.C b/applications/test/findSphereFeatureEdges-octree/Test-findSphereFeatureEdges-octree.C new file mode 100644 index 0000000000..4842756d48 --- /dev/null +++ b/applications/test/findSphereFeatureEdges-octree/Test-findSphereFeatureEdges-octree.C @@ -0,0 +1,134 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "fvMesh.H" +#include "IStringStream.H" +#include "indexedOctree.H" +#include "treeDataEdge.H" +#include "OFstream.H" +#include "extendedFeatureEdgeMesh.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Main program: + +int main(int argc, char *argv[]) +{ + +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + fileName sFeatFileName("unit_cube_rotated.extendedFeatureEdgeMesh"); + + extendedFeatureEdgeMesh efem + ( + IOobject + ( + sFeatFileName, + runTime.time().constant(), + "extendedFeatureEdgeMesh", + runTime.time(), + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + // Slightly extended bb. Slightly off-centred just so on symmetric + // geometry there are less face/edge aligned items. + treeBoundBox bb + ( + efem.points() + ); + + bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + + labelList allEdges(identity(efem.edges().size())); + + indexedOctree edgeTree + ( + treeDataEdge + ( + false, // cachebb + efem.edges(), // edges + efem.points(), // points + allEdges // selected edges + ), + bb, // bb + 8, // maxLevel + 10, // leafsize + 3.0 // duplicity + ); + + Info<< "Points: " << efem.points() << nl << endl; + Info<< "Edges: " << efem.edges() << nl << endl; + + Info<< "Find edge labels within sphere from point (0, 0, 0):" << endl; + + Info<< " Radius = 0 : " + << edgeTree.findSphere(point(0, 0, 0), 0) << endl; + + Info<< " Radius = 1 : " + << edgeTree.findSphere(point(0, 0, 0), 1) << endl; + + Info<< " Radius = root(1.5) : " + << edgeTree.findSphere(point(0, 0, 0), 1.5) << endl; + + Info<< " Radius = root(2) : " + << edgeTree.findSphere(point(0, 0, 0), 2) << endl; + + Info<< " Radius = root(0.5) : " + << edgeTree.findSphere(point(1, 0, 0), 0.5) << endl; + + Info<< " Radius = root(0.25) : " + << edgeTree.findSphere(point(0, 0, 0.5), 0.25) << endl; + + treeBoundBox tbb(point(0,0,0), point(0.1,0.1,0.1)); + Info<< " Box = " << tbb << " : " + << edgeTree.findBox(tbb) << endl; + + treeBoundBox tbb1(point(0,0,0), point(1,1,0.1)); + Info<< " Box = " << tbb1 << " : " + << edgeTree.findBox(tbb1) << endl; + + treeBoundBox tbb2(point(0.3,0,0), point(1,0.3,1)); + Info<< " Box = " << tbb2 << " : " + << edgeTree.findBox(tbb2) << endl; + + treeBoundBox tbb3(point(-0.2,0.5,0), point(0.3,0.9,1)); + Info<< " Box = " << tbb3 << " : " + << edgeTree.findBox(tbb3) << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/momentOfInertia/Test-momentOfInertia.C b/applications/test/momentOfInertia/Test-momentOfInertia.C index 0a825df73f..73aaee649b 100644 --- a/applications/test/momentOfInertia/Test-momentOfInertia.C +++ b/applications/test/momentOfInertia/Test-momentOfInertia.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,7 +35,7 @@ Description #include "polyMesh.H" #include "ListOps.H" #include "face.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "triFaceList.H" #include "OFstream.H" #include "meshTools.H" diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C index 01c9707958..b3d9d3f515 100644 --- a/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C +++ b/applications/utilities/mesh/manipulation/createBaffles/createBaffles.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -290,18 +290,6 @@ int main(int argc, char *argv[]) // Whether first use of face (modify) or consecutive (add) PackedBoolList modifiedFace(mesh.nFaces()); - // Never modify coupled faces - forAll(patches, patchI) - { - const polyPatch& pp = patches[patchI]; - if (pp.coupled()) - { - forAll(pp, i) - { - modifiedFace[pp.start()+i] = 1; - } - } - } label nModified = 0; forAll(newMasterPatches, i) diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C index a45d385894..d0a3d87276 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C +++ b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -254,5 +254,11 @@ Foam::UPstream::commsTypes Foam::UPstream::defaultCommsType commsTypeNames.read(debug::optimisationSwitches().lookup("commsType")) ); +// Number of polling cycles in processor updates +int Foam::UPstream::nPollProcInterfaces +( + debug::optimisationSwitch("nPollProcInterfaces", 0) +); + // ************************************************************************* // diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H index 6ae1dfb368..1ad4aa99cc 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H +++ b/src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -242,6 +242,8 @@ public: //- Default commsType static commsTypes defaultCommsType; + //- Number of polling cycles in processor updates + static int nPollProcInterfaces; // Constructors @@ -273,6 +275,9 @@ public: //- Wait until all requests (from start onwards) have finished. static void waitRequests(const label start = 0); + //- Wait until request i has finished. + static void waitRequest(const label i); + //- Non-blocking comms: has request i finished? static bool finishedRequest(const label i); diff --git a/src/OpenFOAM/global/argList/argList.C b/src/OpenFOAM/global/argList/argList.C index 5d8758e2d2..16f6157a22 100644 --- a/src/OpenFOAM/global/argList/argList.C +++ b/src/OpenFOAM/global/argList/argList.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -772,10 +772,11 @@ Foam::argList::argList Info<< "Roots : " << roots << nl; } Info<< "Pstream initialized with:" << nl - << " floatTransfer : " << Pstream::floatTransfer << nl - << " nProcsSimpleSum : " << Pstream::nProcsSimpleSum << nl - << " commsType : " - << Pstream::commsTypeNames[Pstream::defaultCommsType] + << " floatTransfer : " << Pstream::floatTransfer << nl + << " nProcsSimpleSum : " << Pstream::nProcsSimpleSum << nl + << " commsType : " + << Pstream::commsTypeNames[Pstream::defaultCommsType] << nl + << " polling iterations : " << Pstream::nPollProcInterfaces << endl; } } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H index 040c7f8da8..b8180a2b3e 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -60,6 +60,10 @@ class lduInterfaceField //- Reference to the coupled patch this field is defined for const lduInterface& interface_; + //- Update index used so that updateInterfaceMatrix is called only once + // during the construction of the matrix + bool updatedMatrix_; + // Private Member Functions @@ -69,7 +73,6 @@ class lduInterfaceField //- Disallow default bitwise assignment void operator=(const lduInterfaceField&); - public: //- Runtime type information @@ -81,7 +84,8 @@ public: //- Construct given coupled patch lduInterfaceField(const lduInterface& patch) : - interface_(patch) + interface_(patch), + updatedMatrix_(false) {} @@ -120,6 +124,24 @@ public: ) const {} + //- Whether matrix has been updated + bool updatedMatrix() const + { + return updatedMatrix_; + } + + //- Whether matrix has been updated + bool& updatedMatrix() + { + return updatedMatrix_; + } + + //- Is all data available + virtual bool ready() const + { + return true; + } + //- Update result field based on interface functionality virtual void updateInterfaceMatrix ( diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C index 45654d9731..07260b4a91 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -205,7 +205,15 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A) if (debug > 1) { WarningIn("lduMatrix::operator+=(const lduMatrix& A)") - << "Unknown matrix type combination" + << "Unknown matrix type combination" << nl + << " this :" + << " diagonal:" << diagonal() + << " symmetric:" << symmetric() + << " asymmetric:" << asymmetric() << nl + << " A :" + << " diagonal:" << A.diagonal() + << " symmetric:" << A.symmetric() + << " asymmetric:" << A.asymmetric() << endl; } } @@ -276,7 +284,15 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A) if (debug > 1) { WarningIn("lduMatrix::operator-=(const lduMatrix& A)") - << "Unknown matrix type combination" + << "Unknown matrix type combination" << nl + << " this :" + << " diagonal:" << diagonal() + << " symmetric:" << symmetric() + << " asymmetric:" << asymmetric() << nl + << " A :" + << " diagonal:" << A.diagonal() + << " symmetric:" << A.symmetric() + << " asymmetric:" << A.asymmetric() << endl; } } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C index c401c3e477..62d00b7626 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -87,7 +87,7 @@ void Foam::lduMatrix::initMatrixInterfaces } else { - FatalErrorIn("lduMatrix::initMatrixInterfaces") + FatalErrorIn("lduMatrix::initMatrixInterfaces(..)") << "Unsuported communications type " << Pstream::commsTypeNames[Pstream::defaultCommsType] << exit(FatalError); @@ -104,22 +104,8 @@ void Foam::lduMatrix::updateMatrixInterfaces const direction cmpt ) const { - if - ( - Pstream::defaultCommsType == Pstream::blocking - || Pstream::defaultCommsType == Pstream::nonBlocking - ) + if (Pstream::defaultCommsType == Pstream::blocking) { - // Block until all sends/receives have been finished - if - ( - Pstream::parRun() - && Pstream::defaultCommsType == Pstream::nonBlocking - ) - { - UPstream::waitRequests(); - } - forAll(interfaces, interfaceI) { if (interfaces.set(interfaceI)) @@ -136,6 +122,87 @@ void Foam::lduMatrix::updateMatrixInterfaces } } } + else if (Pstream::defaultCommsType == Pstream::nonBlocking) + { + // Try and consume interfaces as they become available + bool allUpdated = false; + + for (label i = 0; i < UPstream::nPollProcInterfaces; i++) + { + allUpdated = true; + + forAll(interfaces, interfaceI) + { + if (interfaces.set(interfaceI)) + { + if (!interfaces[interfaceI].updatedMatrix()) + { + if (interfaces[interfaceI].ready()) + { + interfaces[interfaceI].updateInterfaceMatrix + ( + psiif, + result, + *this, + coupleCoeffs[interfaceI], + cmpt, + Pstream::defaultCommsType + ); + } + else + { + allUpdated = false; + } + } + } + } + + if (allUpdated) + { + break; + } + } + + // Block for everything + if (Pstream::parRun()) + { + if (allUpdated) + { + // All received. Just remove all storage of requests + // Note that we don't know what starting number of requests + // was before start of sends and receives (since set from + // initMatrixInterfaces) so set to 0 and loose any in-flight + // requests. + UPstream::resetRequests(0); + } + else + { + // Block for all requests and remove storage + UPstream::waitRequests(); + } + } + + // Consume + forAll(interfaces, interfaceI) + { + if + ( + interfaces.set(interfaceI) + && !interfaces[interfaceI].updatedMatrix() + ) + { + interfaces[interfaceI].updateInterfaceMatrix + ( + psiif, + result, + *this, + coupleCoeffs[interfaceI], + cmpt, + Pstream::defaultCommsType + ); + } + } + } else if (Pstream::defaultCommsType == Pstream::scheduled) { const lduSchedule& patchSchedule = this->patchSchedule(); @@ -199,7 +266,7 @@ void Foam::lduMatrix::updateMatrixInterfaces } else { - FatalErrorIn("lduMatrix::updateMatrixInterfaces") + FatalErrorIn("lduMatrix::updateMatrixInterfaces(..)") << "Unsuported communications type " << Pstream::commsTypeNames[Pstream::defaultCommsType] << exit(FatalError); diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C index 3b6d1ee026..b8fa744ff6 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,11 +80,38 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate const Pstream::commsTypes commsType ) const { - procInterface_.compressedSend - ( - commsType, - procInterface_.interfaceInternalField(psiInternal)() - ); + procInterface_.interfaceInternalField(psiInternal, scalarSendBuf_); + + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) + { + // Fast path. + scalarReceiveBuf_.setSize(scalarSendBuf_.size()); + outstandingRecvRequest_ = UPstream::nRequests(); + IPstream::read + ( + Pstream::nonBlocking, + procInterface_.neighbProcNo(), + reinterpret_cast(scalarReceiveBuf_.begin()), + scalarReceiveBuf_.byteSize(), + procInterface_.tag() + ); + + outstandingSendRequest_ = UPstream::nRequests(); + OPstream::write + ( + Pstream::nonBlocking, + procInterface_.neighbProcNo(), + reinterpret_cast(scalarSendBuf_.begin()), + scalarSendBuf_.byteSize(), + procInterface_.tag() + ); + } + else + { + procInterface_.compressedSend(commsType, scalarSendBuf_); + } + + const_cast(*this).updatedMatrix() = false; } @@ -98,18 +125,54 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix const Pstream::commsTypes commsType ) const { - scalarField pnf - ( - procInterface_.compressedReceive(commsType, coeffs.size()) - ); - transformCoupleField(pnf, cmpt); + if (updatedMatrix()) + { + return; + } const labelUList& faceCells = procInterface_.faceCells(); - forAll(faceCells, elemI) + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) { - result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + UPstream::waitRequest(outstandingRecvRequest_); + } + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + + // Consume straight from scalarReceiveBuf_ + + // Transform according to the transformation tensor + transformCoupleField(scalarReceiveBuf_, cmpt); + + // Multiply the field by coefficients and add into the result + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI]; + } } + else + { + scalarField pnf + ( + procInterface_.compressedReceive(commsType, coeffs.size()) + ); + transformCoupleField(pnf, cmpt); + + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + } + } + + const_cast(*this).updatedMatrix() = true; } diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H index c1bdac0246..48d64baacf 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -65,6 +65,22 @@ class processorGAMGInterfaceField int rank_; + // Sending and receiving + + //- Outstanding request + mutable label outstandingSendRequest_; + + //- Outstanding request + mutable label outstandingRecvRequest_; + + //- Scalar send buffer + mutable Field scalarSendBuf_; + + //- Scalar receive buffer + mutable Field scalarReceiveBuf_; + + + // Private Member Functions //- Disallow default bitwise copy construct diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H index 8f25335415..ed6c7ca746 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterface.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -188,6 +188,14 @@ public: const UList& internalData ) const; + //- Get the interface internal field of the given field + template + void interfaceInternalField + ( + const UList& internalData, + List& + ) const; + //- Return the values of the given internal data adjacent to // the interface as a field virtual tmp interfaceInternalField diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C index 0d42d04cef..4890618ab1 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,14 +34,24 @@ Foam::tmp > Foam::GAMGInterface::interfaceInternalField ) const { tmp > tresult(new Field(size())); - Field& result = tresult(); + interfaceInternalField(iF, tresult()); + return tresult; +} + + +template +void Foam::GAMGInterface::interfaceInternalField +( + const UList& iF, + List& result +) const +{ + result.setSize(size()); forAll(result, elemI) { result[elemI] = iF[faceCells_[elemI]]; } - - return tresult; } diff --git a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H index a237335ccc..91e3926c7f 100644 --- a/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H +++ b/src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,7 +43,7 @@ SourceFiles #include "triFace.H" #include "edge.H" #include "pointField.H" -#include "tetPointRef.H" +#include "tetrahedron.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H index 181a2d1122..e1dfbfbf63 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -42,7 +42,7 @@ SourceFiles #include "polyMesh.H" #include "coupledPolyPatch.H" #include "syncTools.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "tetIndices.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H index 76f82fa78c..f0dea9aad8 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -38,7 +38,7 @@ SourceFiles #define tetIndices_H #include "label.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "triPointRef.H" #include "polyMesh.H" #include "triFace.H" diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index 84fd20affc..b5d12a042c 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -419,146 +419,148 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors << gAverage(half1Ctrs) << endl; } - switch (transform_) + if (half0Ctrs.size()) { - case ROTATIONAL: + switch (transform_) { - vector n0 = findFaceMaxRadius(half0Ctrs); - vector n1 = -findFaceMaxRadius(half1Ctrs); - n0 /= mag(n0) + VSMALL; - n1 /= mag(n1) + VSMALL; - - if (debug) + case ROTATIONAL: { - Pout<< "cyclicPolyPatch::getCentresAndAnchors :" - << " patch:" << name() - << " Specified rotation :" - << " n0:" << n0 << " n1:" << n1 << endl; - } + vector n0 = findFaceMaxRadius(half0Ctrs); + vector n1 = -findFaceMaxRadius(half1Ctrs); + n0 /= mag(n0) + VSMALL; + n1 /= mag(n1) + VSMALL; - // Extended tensor from two local coordinate systems calculated - // using normal and rotation axis - const tensor E0 - ( - rotationAxis_, - (n0 ^ rotationAxis_), - n0 - ); - const tensor E1 - ( - rotationAxis_, - (-n1 ^ rotationAxis_), - -n1 - ); - const tensor revT(E1.T() & E0); - - // Rotation - forAll(half0Ctrs, faceI) - { - half0Ctrs[faceI] = - Foam::transform - ( - revT, - half0Ctrs[faceI] - rotationCentre_ - ) - + rotationCentre_; - anchors0[faceI] = - Foam::transform - ( - revT, - anchors0[faceI] - rotationCentre_ - ) - + rotationCentre_; - } - - break; - } - case TRANSLATIONAL: - { - // Transform 0 points. - - if (debug) - { - Pout<< "cyclicPolyPatch::getCentresAndAnchors :" - << " patch:" << name() - << "Specified translation : " << separationVector_ - << endl; - } - - // Note: getCentresAndAnchors gets called on the slave side - // so separationVector is owner-slave points. - - half0Ctrs -= separationVector_; - anchors0 -= separationVector_; - break; - } - default: - { - // Assumes that cyclic is planar. This is also the initial - // condition for patches without faces. - - // Determine the face with max area on both halves. These - // two faces are used to determine the transformation tensors - label max0I = findMaxArea(pp0.points(), pp0); - vector n0 = pp0[max0I].normal(pp0.points()); - n0 /= mag(n0) + VSMALL; - - label max1I = findMaxArea(pp1.points(), pp1); - vector n1 = pp1[max1I].normal(pp1.points()); - n1 /= mag(n1) + VSMALL; - - if (mag(n0 & n1) < 1-matchTolerance()) - { if (debug) { Pout<< "cyclicPolyPatch::getCentresAndAnchors :" << " patch:" << name() - << " Detected rotation :" + << " Specified rotation :" << " n0:" << n0 << " n1:" << n1 << endl; } - // Rotation (around origin) - const tensor revT(rotationTensor(n0, -n1)); + // Extended tensor from two local coordinate systems calculated + // using normal and rotation axis + const tensor E0 + ( + rotationAxis_, + (n0 ^ rotationAxis_), + n0 + ); + const tensor E1 + ( + rotationAxis_, + (-n1 ^ rotationAxis_), + -n1 + ); + const tensor revT(E1.T() & E0); // Rotation forAll(half0Ctrs, faceI) { - half0Ctrs[faceI] = Foam::transform - ( - revT, - half0Ctrs[faceI] - ); - anchors0[faceI] = Foam::transform - ( - revT, - anchors0[faceI] - ); + half0Ctrs[faceI] = + Foam::transform + ( + revT, + half0Ctrs[faceI] - rotationCentre_ + ) + + rotationCentre_; + anchors0[faceI] = + Foam::transform + ( + revT, + anchors0[faceI] - rotationCentre_ + ) + + rotationCentre_; } - } - else - { - // Parallel translation. Get average of all used points. - const point ctr0(sum(pp0.localPoints())/pp0.nPoints()); - const point ctr1(sum(pp1.localPoints())/pp1.nPoints()); + break; + } + case TRANSLATIONAL: + { + // Transform 0 points. if (debug) { Pout<< "cyclicPolyPatch::getCentresAndAnchors :" << " patch:" << name() - << " Detected translation :" - << " n0:" << n0 << " n1:" << n1 - << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl; + << "Specified translation : " << separationVector_ + << endl; } - half0Ctrs += ctr1 - ctr0; - anchors0 += ctr1 - ctr0; + // Note: getCentresAndAnchors gets called on the slave side + // so separationVector is owner-slave points. + + half0Ctrs -= separationVector_; + anchors0 -= separationVector_; + break; + } + default: + { + // Assumes that cyclic is planar. This is also the initial + // condition for patches without faces. + + // Determine the face with max area on both halves. These + // two faces are used to determine the transformation tensors + label max0I = findMaxArea(pp0.points(), pp0); + vector n0 = pp0[max0I].normal(pp0.points()); + n0 /= mag(n0) + VSMALL; + + label max1I = findMaxArea(pp1.points(), pp1); + vector n1 = pp1[max1I].normal(pp1.points()); + n1 /= mag(n1) + VSMALL; + + if (mag(n0 & n1) < 1-matchTolerance()) + { + if (debug) + { + Pout<< "cyclicPolyPatch::getCentresAndAnchors :" + << " patch:" << name() + << " Detected rotation :" + << " n0:" << n0 << " n1:" << n1 << endl; + } + + // Rotation (around origin) + const tensor revT(rotationTensor(n0, -n1)); + + // Rotation + forAll(half0Ctrs, faceI) + { + half0Ctrs[faceI] = Foam::transform + ( + revT, + half0Ctrs[faceI] + ); + anchors0[faceI] = Foam::transform + ( + revT, + anchors0[faceI] + ); + } + } + else + { + // Parallel translation. Get average of all used points. + + const point ctr0(sum(pp0.localPoints())/pp0.nPoints()); + const point ctr1(sum(pp1.localPoints())/pp1.nPoints()); + + if (debug) + { + Pout<< "cyclicPolyPatch::getCentresAndAnchors :" + << " patch:" << name() + << " Detected translation :" + << " n0:" << n0 << " n1:" << n1 + << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl; + } + + half0Ctrs += ctr1 - ctr0; + anchors0 += ctr1 - ctr0; + } + break; } - break; } } - // Calculate typical distance per face tols = matchTolerance()*calcFaceTol(pp1, pp1.points(), half1Ctrs); } @@ -1270,7 +1272,7 @@ bool Foam::cyclicPolyPatch::order rotation.setSize(pp.size()); rotation = 0; - if (pp.empty() || transform_ == NOORDERING) + if (transform_ == NOORDERING) { // No faces, nothing to change. return false; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C index f5f189dbc4..bf103c2254 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,7 @@ License #include "primitiveMesh.H" #include "pyramidPointFaceRef.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "ListOps.H" #include "unitConversion.H" #include "SortableList.H" diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H deleted file mode 100644 index 56382a55c7..0000000000 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H +++ /dev/null @@ -1,52 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ 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 3 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, see . - -Typedef - Foam::tetPointRef - -Description - -\*---------------------------------------------------------------------------*/ - -#ifndef tetPointRef_H -#define tetPointRef_H - -#include "point.H" -#include "tetrahedron.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -typedef tetrahedron tetPointRef; - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -#endif - -// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H similarity index 99% rename from src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H rename to src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H index b23f7540dc..2e1ecf288d 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H @@ -35,9 +35,9 @@ SourceFiles #ifndef tetPoints_H #define tetPoints_H +#include "tetrahedron.H" #include "FixedList.H" #include "treeBoundBox.H" -#include "tetPointRef.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H index a2ccbc250f..1a76022205 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) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -54,6 +54,8 @@ namespace Foam class Istream; class Ostream; +class tetPoints; +class plane; // Forward declaration of friend functions and operators @@ -73,6 +75,7 @@ inline Ostream& operator<< const tetrahedron& ); +typedef tetrahedron tetPointRef; /*---------------------------------------------------------------------------*\ class tetrahedron Declaration @@ -81,10 +84,71 @@ inline Ostream& operator<< template class tetrahedron { +public: + + // Classes for use in sliceWithPlane. What to do with decomposition + // of tet. + + //- Dummy + class dummyOp + { + public: + inline void operator()(const tetPoints&); + }; + + //- Sum resulting volumes + class sumVolOp + { + public: + scalar vol_; + + inline sumVolOp(); + + inline void operator()(const tetPoints&); + }; + + //- Store resulting tets + class storeOp + { + FixedList& tets_; + label& nTets_; + + public: + inline storeOp(FixedList&, label&); + + inline void operator()(const tetPoints&); + }; + +private: + // Private data PointRef a_, b_, c_, d_; + inline static point planeIntersection + ( + const FixedList&, + const tetPoints&, + const label, + const label + ); + + template + inline static void decomposePrism + ( + const FixedList& points, + TetOp& op + ); + + template + inline static void tetSliceWithPlane + ( + const plane& pl, + const tetPoints& tet, + AboveTetOp& aboveOp, + BelowTetOp& belowOp + ); + public: @@ -184,6 +248,16 @@ public: //- Return true if point is inside tetrahedron inline bool inside(const point& pt) const; + //- Decompose tet into tets above and below plane + template + inline void sliceWithPlane + ( + const plane& pl, + AboveTetOp& aboveOp, + BelowTetOp& belowOp + ) const; + + //- Return (min)containment sphere, i.e. the smallest sphere with // all points inside. Returns pointHit with: // - hit : if sphere is equal to circumsphere diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 40b38f8067..8a69e88e04 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) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,8 @@ License #include "triangle.H" #include "IOstreams.H" -#include "triPointRef.H" +#include "tetPoints.H" +#include "plane.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -492,6 +493,483 @@ bool Foam::tetrahedron::inside(const point& pt) const } +template +inline void Foam::tetrahedron::dummyOp::operator() +( + const tetPoints& +) +{} + + +template +inline Foam::tetrahedron::sumVolOp::sumVolOp() +: + vol_(0.0) +{} + + +template +inline void Foam::tetrahedron::sumVolOp::operator() +( + const tetPoints& tet +) +{ + vol_ += tet.tet().mag(); +} + + +template +inline Foam::tetrahedron::storeOp::storeOp +( + FixedList& tets, + label& nTets +) +: + tets_(tets), + nTets_(nTets) +{} + + +template +inline void Foam::tetrahedron::storeOp::operator() +( + const tetPoints& tet +) +{ + tets_[nTets_++] = tet; +} + + +template +inline Foam::point Foam::tetrahedron::planeIntersection +( + const FixedList& d, + const tetPoints& t, + const label negI, + const label posI +) +{ + return + (d[posI]*t[negI] - d[negI]*t[posI]) + / (-d[negI]+d[posI]); +} + + +template +template +inline void Foam::tetrahedron::decomposePrism +( + const FixedList& points, + TetOp& op +) +{ + op(tetPoints(points[1], points[3], points[2], points[0])); + op(tetPoints(points[1], points[2], points[3], points[4])); + op(tetPoints(points[4], points[2], points[3], points[5])); +} + + +template +template +inline void Foam::tetrahedron:: +tetSliceWithPlane +( + const plane& pl, + const tetPoints& tet, + AboveTetOp& aboveOp, + BelowTetOp& belowOp +) +{ + // Distance to plane + FixedList d; + label nPos = 0; + forAll(tet, i) + { + d[i] = ((tet[i]-pl.refPoint()) & pl.normal()); + if (d[i] > 0) + { + nPos++; + } + } + + if (nPos == 4) + { + aboveOp(tet); + } + else if (nPos == 3) + { + // Sliced into below tet and above prism. Prism gets split into + // two tets. + + // Find the below tet + label i0 = -1; + forAll(d, i) + { + if (d[i] <= 0) + { + i0 = i; + break; + } + } + + label i1 = d.fcIndex(i0); + label i2 = d.fcIndex(i1); + label i3 = d.fcIndex(i2); + + point p01 = planeIntersection(d, tet, i0, i1); + point p02 = planeIntersection(d, tet, i0, i2); + point p03 = planeIntersection(d, tet, i0, i3); + + // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad + // ,, 1 : ,, inwards pointing triad + // ,, 2 : ,, outwards pointing triad + // ,, 3 : ,, inwards pointing triad + + //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl; + + if (i0 == 0 || i0 == 2) + { + tetPoints t(tet[i0], p01, p02, p03); + //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; + //checkTet(t, "nPos 3, belowTet i0==0 or 2"); + belowOp(t); + + // Prism + FixedList p; + p[0] = tet[i1]; + p[1] = tet[i3]; + p[2] = tet[i2]; + p[3] = p01; + p[4] = p03; + p[5] = p02; + //Pout<< " aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + else + { + tetPoints t(p01, p02, p03, tet[i0]); + //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; + //checkTet(t, "nPos 3, belowTet i0==1 or 3"); + belowOp(t); + + // Prism + FixedList p; + p[0] = tet[i3]; + p[1] = tet[i1]; + p[2] = tet[i2]; + p[3] = p03; + p[4] = p01; + p[5] = p02; + //Pout<< " aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + } + else if (nPos == 2) + { + // Tet cut into two prisms. Determine the positive one. + label pos0 = -1; + label pos1 = -1; + label neg0 = -1; + label neg1 = -1; + forAll(d, i) + { + if (d[i] > 0) + { + if (pos0 == -1) + { + pos0 = i; + } + else + { + pos1 = i; + } + } + else + { + if (neg0 == -1) + { + neg0 = i; + } + else + { + neg1 = i; + } + } + } + + //Pout<< "Split 2pos tet " << tet << " d:" << d + // << " around pos0:" << pos0 << " pos1:" << pos1 + // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl; + + const edge posEdge(pos0, pos1); + + if (posEdge == edge(0, 1)) + { + point p02 = planeIntersection(d, tet, 0, 2); + point p03 = planeIntersection(d, tet, 0, 3); + point p12 = planeIntersection(d, tet, 1, 2); + point p13 = planeIntersection(d, tet, 1, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[0]; + p[1] = p02; + p[2] = p03; + p[3] = tet[1]; + p[4] = p12; + p[5] = p13; + //Pout<< " 01 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[2]; + p[1] = p02; + p[2] = p12; + p[3] = tet[3]; + p[4] = p03; + p[5] = p13; + //Pout<< " 01 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else if (posEdge == edge(1, 2)) + { + point p01 = planeIntersection(d, tet, 0, 1); + point p13 = planeIntersection(d, tet, 1, 3); + point p02 = planeIntersection(d, tet, 0, 2); + point p23 = planeIntersection(d, tet, 2, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[1]; + p[1] = p01; + p[2] = p13; + p[3] = tet[2]; + p[4] = p02; + p[5] = p23; + //Pout<< " 12 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[3]; + p[1] = p23; + p[2] = p13; + p[3] = tet[0]; + p[4] = p02; + p[5] = p01; + //Pout<< " 12 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else if (posEdge == edge(2, 0)) + { + point p01 = planeIntersection(d, tet, 0, 1); + point p03 = planeIntersection(d, tet, 0, 3); + point p12 = planeIntersection(d, tet, 1, 2); + point p23 = planeIntersection(d, tet, 2, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[2]; + p[1] = p12; + p[2] = p23; + p[3] = tet[0]; + p[4] = p01; + p[5] = p03; + //Pout<< " 20 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[1]; + p[1] = p12; + p[2] = p01; + p[3] = tet[3]; + p[4] = p23; + p[5] = p03; + //Pout<< " 20 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else if (posEdge == edge(0, 3)) + { + point p01 = planeIntersection(d, tet, 0, 1); + point p02 = planeIntersection(d, tet, 0, 2); + point p13 = planeIntersection(d, tet, 1, 3); + point p23 = planeIntersection(d, tet, 2, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[3]; + p[1] = p23; + p[2] = p13; + p[3] = tet[0]; + p[4] = p02; + p[5] = p01; + //Pout<< " 03 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[2]; + p[1] = p23; + p[2] = p02; + p[3] = tet[1]; + p[4] = p13; + p[5] = p01; + //Pout<< " 03 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else if (posEdge == edge(1, 3)) + { + point p01 = planeIntersection(d, tet, 0, 1); + point p12 = planeIntersection(d, tet, 1, 2); + point p03 = planeIntersection(d, tet, 0, 3); + point p23 = planeIntersection(d, tet, 2, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[1]; + p[1] = p12; + p[2] = p01; + p[3] = tet[3]; + p[4] = p23; + p[5] = p03; + //Pout<< " 13 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[2]; + p[1] = p12; + p[2] = p23; + p[3] = tet[0]; + p[4] = p01; + p[5] = p03; + //Pout<< " 13 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else if (posEdge == edge(2, 3)) + { + point p02 = planeIntersection(d, tet, 0, 2); + point p12 = planeIntersection(d, tet, 1, 2); + point p03 = planeIntersection(d, tet, 0, 3); + point p13 = planeIntersection(d, tet, 1, 3); + // Split the resulting prism + { + FixedList p; + p[0] = tet[2]; + p[1] = p02; + p[2] = p12; + p[3] = tet[3]; + p[4] = p03; + p[5] = p13; + //Pout<< " 23 aboveprism:" << p << endl; + decomposePrism(p, aboveOp); + } + { + FixedList p; + p[0] = tet[0]; + p[1] = p02; + p[2] = p03; + p[3] = tet[1]; + p[4] = p12; + p[5] = p13; + //Pout<< " 23 belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else + { + FatalErrorIn("tetSliceWithPlane(..)") + << "Missed edge:" << posEdge + << abort(FatalError); + } + } + else if (nPos == 1) + { + // Find the positive tet + label i0 = -1; + forAll(d, i) + { + if (d[i] > 0) + { + i0 = i; + break; + } + } + + label i1 = d.fcIndex(i0); + label i2 = d.fcIndex(i1); + label i3 = d.fcIndex(i2); + + point p01 = planeIntersection(d, tet, i0, i1); + point p02 = planeIntersection(d, tet, i0, i2); + point p03 = planeIntersection(d, tet, i0, i3); + + //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl; + + if (i0 == 0 || i0 == 2) + { + tetPoints t(tet[i0], p01, p02, p03); + //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; + //checkTet(t, "nPos 1, aboveTets i0==0 or 2"); + aboveOp(t); + + // Prism + FixedList p; + p[0] = tet[i1]; + p[1] = tet[i3]; + p[2] = tet[i2]; + p[3] = p01; + p[4] = p03; + p[5] = p02; + //Pout<< " belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + else + { + tetPoints t(p01, p02, p03, tet[i0]); + //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; + //checkTet(t, "nPos 1, aboveTets i0==1 or 3"); + aboveOp(t); + + // Prism + FixedList p; + p[0] = tet[i3]; + p[1] = tet[i1]; + p[2] = tet[i2]; + p[3] = p03; + p[4] = p01; + p[5] = p02; + //Pout<< " belowprism:" << p << endl; + decomposePrism(p, belowOp); + } + } + else // nPos == 0 + { + belowOp(tet); + } +} + + +template +template +inline void Foam::tetrahedron::sliceWithPlane +( + const plane& pl, + AboveTetOp& aboveOp, + BelowTetOp& belowOp +) const +{ + tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp); +} + + // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H index f128f5ad51..60106598b1 100644 --- a/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H +++ b/src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/Pstream/dummy/UPstream.C b/src/Pstream/dummy/UPstream.C index f208ac880e..5cd8bbcbcb 100644 --- a/src/Pstream/dummy/UPstream.C +++ b/src/Pstream/dummy/UPstream.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,6 +74,10 @@ void Foam::UPstream::waitRequests(const label start) {} +void Foam::UPstream::waitRequest(const label i) +{} + + bool Foam::UPstream::finishedRequest(const label i) { notImplemented("UPstream::finishedRequest()"); diff --git a/src/Pstream/mpi/UPstream.C b/src/Pstream/mpi/UPstream.C index c85f341d8f..1605eca8ce 100644 --- a/src/Pstream/mpi/UPstream.C +++ b/src/Pstream/mpi/UPstream.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -516,11 +516,55 @@ void Foam::UPstream::waitRequests(const label start) } +void Foam::UPstream::waitRequest(const label i) +{ + if (debug) + { + Pout<< "UPstream::waitRequest : starting wait for request:" << i + << endl; + } + + if (i >= PstreamGlobals::outstandingRequests_.size()) + { + FatalErrorIn + ( + "UPstream::waitRequest(const label)" + ) << "There are " << PstreamGlobals::outstandingRequests_.size() + << " outstanding send requests and you are asking for i=" << i + << nl + << "Maybe you are mixing blocking/non-blocking comms?" + << Foam::abort(FatalError); + } + + int flag; + if + ( + MPI_Wait + ( + &PstreamGlobals::outstandingRequests_[i], + MPI_STATUS_IGNORE + ) + ) + { + FatalErrorIn + ( + "UPstream::waitRequest()" + ) << "MPI_Wait returned with error" << Foam::endl; + } + + if (debug) + { + Pout<< "UPstream::waitRequest : finished wait for request:" << i + << endl; + } +} + + bool Foam::UPstream::finishedRequest(const label i) { if (debug) { - Pout<< "UPstream::waitRequests : starting wait for request:" << i + Pout<< "UPstream::waitRequests : checking finishedRequest:" << i << endl; } @@ -546,7 +590,7 @@ bool Foam::UPstream::finishedRequest(const label i) if (debug) { - Pout<< "UPstream::waitRequests : finished wait for request:" << i + Pout<< "UPstream::waitRequests : finished finishedRequest:" << i << endl; } @@ -554,6 +598,4 @@ bool Foam::UPstream::finishedRequest(const label i) } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - // ************************************************************************* // diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index afac722830..1670de1c92 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,7 +26,7 @@ License #include "polyMeshGeometry.H" #include "polyMeshTetDecomposition.H" #include "pyramidPointFaceRef.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "syncTools.H" #include "unitConversion.H" diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C index 1ec37f1a2f..9f621b24ee 100644 --- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C +++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C @@ -756,6 +756,66 @@ void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType } +void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges +( + const point& sample, + const scalar searchRadiusSqr, + List& info +) const +{ + const PtrList >& edgeTrees = edgeTreesByType(); + + info.setSize(edgeTrees.size()); + + labelList sliceStarts(edgeTrees.size()); + + sliceStarts[0] = externalStart_; + sliceStarts[1] = internalStart_; + sliceStarts[2] = flatStart_; + sliceStarts[3] = openStart_; + sliceStarts[4] = multipleStart_; + + DynamicList dynEdgeHit; + + // Loop over all the feature edge types + forAll(edgeTrees, i) + { + // Pick up all the edges that intersect the search sphere + labelList elems = edgeTrees[i].findSphere + ( + sample, + searchRadiusSqr + ); + + forAll(elems, elemI) + { + label index = elems[elemI]; + label edgeI = edgeTrees[i].shapes().edgeLabels()[index]; + const edge& e = edges()[edgeI]; + + pointHit hitPoint = e.line(points()).nearestDist(sample); + + label hitIndex = index + sliceStarts[i]; + + pointIndexHit nearHit; + + if (!hitPoint.hit()) + { + nearHit = pointIndexHit(true, hitPoint.missPoint(), hitIndex); + } + else + { + nearHit = pointIndexHit(true, hitPoint.hitPoint(), hitIndex); + } + + dynEdgeHit.append(nearHit); + } + } + + info.transfer(dynEdgeHit); +} + + const Foam::indexedOctree& Foam::extendedFeatureEdgeMesh::edgeTree() const { diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H index 82aef34ab0..3268e719c5 100644 --- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H +++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H @@ -271,6 +271,14 @@ public: List& info ) const; + //- Find all the feature edges within searchDistSqr of sample + void allNearestFeatureEdges + ( + const point& sample, + const scalar searchRadiusSqr, + List& info + ) const; + // Access //- Return the index of the start of the convex feature points diff --git a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C index acf50983ef..f530fd1ac9 100644 --- a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -198,6 +198,17 @@ tmp > slicedFvPatchField::patchInternalField() const } +template +void slicedFvPatchField::patchInternalField(Field&) const +{ + notImplemented + ( + "slicedFvPatchField::" + "patchInternalField(Field&)" + ); +} + + template tmp > slicedFvPatchField::patchNeighbourField ( diff --git a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H index 9784c51323..a40d212cd7 100644 --- a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -148,6 +148,9 @@ public: //- Return internal field next to patch as patch field virtual tmp > patchInternalField() const; + //- Return internal field next to patch as patch field + virtual void patchInternalField(Field&) const; + //- Return neighbour coupled given internal cell data virtual tmp > patchNeighbourField ( diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C index 4bf8e67cd7..a44ea4cc5c 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,7 +43,12 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(p, iF), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + sendBuf_(0), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(0), + scalarReceiveBuf_(0) {} @@ -56,7 +61,12 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(p, iF, f), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + sendBuf_(0), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(0), + scalarReceiveBuf_(0) {} @@ -71,7 +81,12 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(ptf, p, iF, mapper), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + sendBuf_(0), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(0), + scalarReceiveBuf_(0) { if (!isA(this->patch())) { @@ -91,6 +106,12 @@ processorFvPatchField::processorFvPatchField << " in file " << this->dimensionedInternalField().objectPath() << exit(FatalIOError); } + if (debug && !ptf.ready()) + { + FatalErrorIn("processorFvPatchField::processorFvPatchField(..)") + << "On patch " << procPatch_.name() << " outstanding request." + << abort(FatalError); + } } @@ -103,7 +124,12 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(p, iF, dict), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + sendBuf_(0), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(0), + scalarReceiveBuf_(0) { if (!isA(p)) { @@ -134,8 +160,20 @@ processorFvPatchField::processorFvPatchField : processorLduInterfaceField(), coupledFvPatchField(ptf), - procPatch_(refCast(ptf.patch())) -{} + procPatch_(refCast(ptf.patch())), + sendBuf_(ptf.sendBuf_.xfer()), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(ptf.scalarSendBuf_.xfer()), + scalarReceiveBuf_(ptf.scalarReceiveBuf_.xfer()) +{ + if (debug && !ptf.ready()) + { + FatalErrorIn("processorFvPatchField::processorFvPatchField(..)") + << "On patch " << procPatch_.name() << " outstanding request." + << abort(FatalError); + } +} template @@ -146,8 +184,20 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(ptf, iF), - procPatch_(refCast(ptf.patch())) -{} + procPatch_(refCast(ptf.patch())), + sendBuf_(0), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + scalarSendBuf_(0), + scalarReceiveBuf_(0) +{ + if (debug && !ptf.ready()) + { + FatalErrorIn("processorFvPatchField::processorFvPatchField(..)") + << "On patch " << procPatch_.name() << " outstanding request." + << abort(FatalError); + } +} // * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // @@ -162,6 +212,13 @@ processorFvPatchField::~processorFvPatchField() template tmp > processorFvPatchField::patchNeighbourField() const { + if (debug && !this->ready()) + { + FatalErrorIn("processorFvPatchField::patchNeighbourField()") + << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } return *this; } @@ -174,7 +231,36 @@ void processorFvPatchField::initEvaluate { if (Pstream::parRun()) { - procPatch_.compressedSend(commsType, this->patchInternalField()()); + this->patchInternalField(sendBuf_); + + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) + { + // Fast path. Receive into *this + this->setSize(sendBuf_.size()); + outstandingRecvRequest_ = UPstream::nRequests(); + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(this->begin()), + this->byteSize(), + procPatch_.tag() + ); + + outstandingSendRequest_ = UPstream::nRequests(); + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(sendBuf_.begin()), + this->byteSize(), + procPatch_.tag() + ); + } + else + { + procPatch_.compressedSend(commsType, sendBuf_); + } } } @@ -187,7 +273,25 @@ void processorFvPatchField::evaluate { if (Pstream::parRun()) { - procPatch_.compressedReceive(commsType, *this); + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) + { + // Fast path. Received into *this + + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + UPstream::waitRequest(outstandingRecvRequest_); + } + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + } + else + { + procPatch_.compressedReceive(commsType, *this); + } if (doTransform()) { @@ -215,11 +319,49 @@ void processorFvPatchField::initInterfaceMatrixUpdate const Pstream::commsTypes commsType ) const { - procPatch_.compressedSend - ( - commsType, - this->patch().patchInternalField(psiInternal)() - ); + this->patch().patchInternalField(psiInternal, scalarSendBuf_); + + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) + { + // Fast path. + if (debug && !this->ready()) + { + FatalErrorIn + ( + "processorFvPatchField::initInterfaceMatrixUpdate(..)" + ) << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } + + + scalarReceiveBuf_.setSize(scalarSendBuf_.size()); + outstandingRecvRequest_ = UPstream::nRequests(); + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarReceiveBuf_.begin()), + scalarReceiveBuf_.byteSize(), + procPatch_.tag() + ); + + outstandingSendRequest_ = UPstream::nRequests(); + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarSendBuf_.begin()), + scalarSendBuf_.byteSize(), + procPatch_.tag() + ); + } + else + { + procPatch_.compressedSend(commsType, scalarSendBuf_); + } + + const_cast&>(*this).updatedMatrix() = false; } @@ -234,22 +376,92 @@ void processorFvPatchField::updateInterfaceMatrix const Pstream::commsTypes commsType ) const { - scalarField pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); - - // Transform according to the transformation tensor - transformCoupleField(pnf, cmpt); - - // Multiply the field by coefficients and add into the result + if (this->updatedMatrix()) + { + return; + } const labelUList& faceCells = this->patch().faceCells(); - forAll(faceCells, elemI) + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) { - result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + UPstream::waitRequest(outstandingRecvRequest_); + } + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + + // Consume straight from scalarReceiveBuf_ + + // Transform according to the transformation tensor + transformCoupleField(scalarReceiveBuf_, cmpt); + + // Multiply the field by coefficients and add into the result + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI]; + } } + else + { + scalarField pnf + ( + procPatch_.compressedReceive(commsType, this->size())() + ); + + // Transform according to the transformation tensor + transformCoupleField(pnf, cmpt); + + // Multiply the field by coefficients and add into the result + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + } + } + + const_cast&>(*this).updatedMatrix() = true; +} + + +template +bool processorFvPatchField::ready() const +{ + if + ( + outstandingSendRequest_ >= 0 + && outstandingSendRequest_ < Pstream::nRequests() + ) + { + bool finished = UPstream::finishedRequest(outstandingSendRequest_); + if (!finished) + { + return false; + } + } + outstandingSendRequest_ = -1; + + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + bool finished = UPstream::finishedRequest(outstandingRecvRequest_); + if (!finished) + { + return false; + } + } + outstandingRecvRequest_ = -1; + + return true; } diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H index 0130909d04..9541e80ea5 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -59,6 +59,22 @@ class processorFvPatchField //- Local reference cast into the processor patch const processorFvPatch& procPatch_; + // Sending and receiving + + //- Send buffer. + mutable Field sendBuf_; + + //- Outstanding request + mutable label outstandingSendRequest_; + + //- Outstanding request + mutable label outstandingRecvRequest_; + + //- Scalar send buffer + mutable Field scalarSendBuf_; + + //- Scalar receive buffer + mutable Field scalarReceiveBuf_; public: @@ -179,6 +195,9 @@ public: const Pstream::commsTypes commsType ) const; + //- Is all data available + virtual bool ready() const; + //- Update result field based on interface functionality virtual void updateInterfaceMatrix ( @@ -221,6 +240,7 @@ public: { return pTraits::rank; } + }; diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C index b94e8049f9..c15b1c1be7 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,11 +43,53 @@ void processorFvPatchField::initInterfaceMatrixUpdate const Pstream::commsTypes commsType ) const { - procPatch_.compressedSend + this->patch().patchInternalField(psiInternal, scalarSendBuf_); + + if ( - commsType, - patch().patchInternalField(psiInternal)() - ); + Pstream::defaultCommsType == Pstream::nonBlocking + && !Pstream::floatTransfer + ) + { + // Fast path. + if (debug && !this->ready()) + { + FatalErrorIn + ( + "processorFvPatchField::initInterfaceMatrixUpdate(..)" + ) << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } + + + scalarReceiveBuf_.setSize(scalarSendBuf_.size()); + outstandingRecvRequest_ = UPstream::nRequests(); + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarReceiveBuf_.begin()), + scalarReceiveBuf_.byteSize(), + procPatch_.tag() + ); + + outstandingSendRequest_ = UPstream::nRequests(); + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarSendBuf_.begin()), + scalarSendBuf_.byteSize(), + procPatch_.tag() + ); + } + else + { + procPatch_.compressedSend(commsType, scalarSendBuf_); + } + + const_cast&>(*this).updatedMatrix() = false; } @@ -62,17 +104,48 @@ void processorFvPatchField::updateInterfaceMatrix const Pstream::commsTypes commsType ) const { - scalarField pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); - - const labelUList& faceCells = patch().faceCells(); - - forAll(faceCells, facei) + if (this->updatedMatrix()) { - result[faceCells[facei]] -= coeffs[facei]*pnf[facei]; + return; } + + const labelUList& faceCells = this->patch().faceCells(); + + if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer) + { + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + UPstream::waitRequest(outstandingRecvRequest_); + } + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + + + // Consume straight from scalarReceiveBuf_ + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI]; + } + } + else + { + scalarField pnf + ( + procPatch_.compressedReceive(commsType, this->size())() + ); + + forAll(faceCells, elemI) + { + result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + } + } + + const_cast&>(*this).updatedMatrix() = true; } diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C index 63038d3dd3..9266ab2749 100644 --- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -190,6 +190,13 @@ Foam::fvPatchField::patchInternalField() const } +template +void Foam::fvPatchField::patchInternalField(Field& pif) const +{ + patch_.patchInternalField(internalField_, pif); +} + + template void Foam::fvPatchField::autoMap ( diff --git a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H index 4d9d2e0f88..369246d3b0 100644 --- a/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -352,6 +352,9 @@ public: //- Return internal field next to patch as patch field virtual tmp > patchInternalField() const; + //- Return internal field next to patch as patch field + virtual void patchInternalField(Field&) const; + //- Return patchField on the opposite patch of a coupled patch virtual tmp > patchNeighbourField() const { diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H index be7ee1bda4..28a07343e5 100644 --- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -239,6 +239,10 @@ public: template tmp > patchInternalField(const UList&) const; + //- Return given internal field next to patch as patch field + template + void patchInternalField(const UList&, Field&) const; + //- Return the corresponding patchField of the named field template const typename GeometricField::PatchFieldType& patchField diff --git a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchTemplates.C b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchTemplates.C index 5d47411fa0..11aef36b05 100644 --- a/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchTemplates.C +++ b/src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchTemplates.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,6 +47,24 @@ Foam::tmp > Foam::fvPatch::patchInternalField } +template +void Foam::fvPatch::patchInternalField +( + const UList& f, + Field& pif +) const +{ + pif.setSize(size()); + + const labelUList& faceCells = this->faceCells(); + + forAll(pif, facei) + { + pif[facei] = f[faceCells[facei]]; + } +} + + template const typename GeometricField::PatchFieldType& Foam::fvPatch::patchField ( diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C index 13828ff7cf..9d75e32336 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,7 +25,6 @@ License #include "cellPointWeight.H" #include "polyMesh.H" -#include "tetPointRef.H" #include "polyMeshTetDecomposition.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // diff --git a/src/lagrangian/basic/Cloud/Cloud.H b/src/lagrangian/basic/Cloud/Cloud.H index 56b4759f1a..dd7358d2c7 100644 --- a/src/lagrangian/basic/Cloud/Cloud.H +++ b/src/lagrangian/basic/Cloud/Cloud.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,9 +41,6 @@ SourceFiles #include "IOField.H" #include "CompactIOField.H" #include "polyMesh.H" -#include "indexedOctree.H" -#include "treeDataCell.H" -#include "tetPointRef.H" #include "PackedBoolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 08913b4145..8ab17472d0 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,11 +35,10 @@ Description #include "vector.H" #include "Cloud.H" #include "IDLList.H" -#include "labelList.H" #include "pointField.H" #include "faceList.H" #include "OFstream.H" -#include "tetPointRef.H" +#include "tetrahedron.H" #include "FixedList.H" #include "polyMeshTetDecomposition.H" diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index 67d4fec174..f70381471b 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -1144,6 +1144,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement ); // Reset if differs + // 1. take max forAll(syncPatchNLayers, i) { if (syncPatchNLayers[i] != patchNLayers[i]) @@ -1174,6 +1175,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement ); // Reset if differs + // 2. take min forAll(syncPatchNLayers, i) { if (syncPatchNLayers[i] != patchNLayers[i]) @@ -2726,7 +2728,7 @@ void Foam::autoLayerDriver::addLayers Info<< "Writing shrunk mesh to " << meshRefiner_.timeName() << endl; // See comment in autoSnapDriver why we should not remove meshPhi - // using mesh.clearPout(). + // using mesh.clearOut(). meshRefiner_.write ( diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C index d394582236..62ba5f2917 100644 --- a/src/meshTools/indexedOctree/treeDataEdge.C +++ b/src/meshTools/indexedOctree/treeDataEdge.C @@ -128,14 +128,37 @@ bool Foam::treeDataEdge::overlaps const treeBoundBox& cubeBb ) const { - if (cacheBb_) + const edge& e = edges_[edgeLabels_[index]]; + + const point& start = points_[e.start()]; + const point& end = points_[e.end()]; + + point intersect; + + return cubeBb.intersects(start, end, intersect); +} + + +// Check if any point on shape is inside sphere. +bool Foam::treeDataEdge::overlaps +( + const label index, + const point& centre, + const scalar radiusSqr +) const +{ + const edge& e = edges_[edgeLabels_[index]]; + + const pointHit nearHit = e.line(points_).nearestDist(centre); + + const scalar distSqr = sqr(nearHit.distance()); + + if (distSqr <= radiusSqr) { - return cubeBb.overlaps(bbs_[index]); - } - else - { - return cubeBb.overlaps(calcBb(edgeLabels_[index])); + return true; } + + return false; } diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H index c6af44e9ed..9a8deadb73 100644 --- a/src/meshTools/indexedOctree/treeDataEdge.H +++ b/src/meshTools/indexedOctree/treeDataEdge.H @@ -148,6 +148,14 @@ public: const treeBoundBox& sampleBb ) const; + //- Does (bb of) shape at index overlap bb + bool overlaps + ( + const label index, + const point& centre, + const scalar radiusSqr + ) const; + //- Calculates nearest (to sample) point in shape. // Returns actual point and distance (squared) void findNearest diff --git a/src/meshTools/momentOfInertia/momentOfInertia.H b/src/meshTools/momentOfInertia/momentOfInertia.H index a4d03ad808..f21f1f2983 100644 --- a/src/meshTools/momentOfInertia/momentOfInertia.H +++ b/src/meshTools/momentOfInertia/momentOfInertia.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,7 +37,6 @@ SourceFiles #ifndef momentOfInertia_H #define momentOfInertia_H -#include "tetPointRef.H" #include "triFaceList.H" #include "triSurface.H" #include "polyMesh.H" diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C index d36730a0ee..16d311815c 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,436 +24,27 @@ License \*---------------------------------------------------------------------------*/ #include "tetOverlapVolume.H" +#include "tetrahedron.H" +#include "tetPoints.H" +#include "polyMesh.H" +#include "OFstream.H" +#include "treeBoundBox.H" +#include "indexedOctree.H" +#include "treeDataCell.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::tetOverlapVolume, 0); + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::tetOverlapVolume::tetOverlapVolume() {} -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::tetOverlapVolume::~tetOverlapVolume() -{} - - // * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // -Foam::point Foam::tetOverlapVolume::planeIntersection -( - const FixedList& d, - const tetPoints& t, - const label negI, - const label posI -) -{ - return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI]+d[posI]); -} - - -template -inline void Foam::tetOverlapVolume::decomposePrism -( - const FixedList& points, - TetOp& op -) -{ - op(tetPoints(points[1], points[3], points[2], points[0])); - op(tetPoints(points[1], points[2], points[3], points[4])); - op(tetPoints(points[4], points[2], points[3], points[5])); -} - - -template -inline void Foam::tetOverlapVolume::tetSliceWithPlane -( - const tetPoints& tet, - const plane& pl, - - AboveTetOp& aboveOp, - BelowTetOp& belowOp -) -{ - // Distance to plane - FixedList d; - label nPos = 0; - forAll(tet, i) - { - d[i] = ((tet[i]-pl.refPoint()) & pl.normal()); - if (d[i] > 0) - { - nPos++; - } - } - - if (nPos == 4) - { - aboveOp(tet); - } - else if (nPos == 3) - { - // Sliced into below tet and above prism. Prism gets split into - // two tets. - - // Find the below tet - label i0 = -1; - forAll(d, i) - { - if (d[i] <= 0) - { - i0 = i; - break; - } - } - - label i1 = d.fcIndex(i0); - label i2 = d.fcIndex(i1); - label i3 = d.fcIndex(i2); - - point p01 = planeIntersection(d, tet, i0, i1); - point p02 = planeIntersection(d, tet, i0, i2); - point p03 = planeIntersection(d, tet, i0, i3); - - // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad - // ,, 1 : ,, inwards pointing triad - // ,, 2 : ,, outwards pointing triad - // ,, 3 : ,, inwards pointing triad - - //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl; - - if (i0 == 0 || i0 == 2) - { - tetPoints t(tet[i0], p01, p02, p03); - //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 3, belowTet i0==0 or 2"); - belowOp(t); - - // Prism - FixedList p; - p[0] = tet[i1]; - p[1] = tet[i3]; - p[2] = tet[i2]; - p[3] = p01; - p[4] = p03; - p[5] = p02; - //Pout<< " aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - else - { - tetPoints t(p01, p02, p03, tet[i0]); - //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 3, belowTet i0==1 or 3"); - belowOp(t); - - // Prism - FixedList p; - p[0] = tet[i3]; - p[1] = tet[i1]; - p[2] = tet[i2]; - p[3] = p03; - p[4] = p01; - p[5] = p02; - //Pout<< " aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - } - else if (nPos == 2) - { - // Tet cut into two prisms. Determine the positive one. - label pos0 = -1; - label pos1 = -1; - label neg0 = -1; - label neg1 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - if (pos0 == -1) - { - pos0 = i; - } - else - { - pos1 = i; - } - } - else - { - if (neg0 == -1) - { - neg0 = i; - } - else - { - neg1 = i; - } - } - } - - //Pout<< "Split 2pos tet " << tet << " d:" << d - // << " around pos0:" << pos0 << " pos1:" << pos1 - // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl; - - const edge posEdge(pos0, pos1); - - if (posEdge == edge(0, 1)) - { - point p02 = planeIntersection(d, tet, 0, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p12 = planeIntersection(d, tet, 1, 2); - point p13 = planeIntersection(d, tet, 1, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[0]; - p[1] = p02; - p[2] = p03; - p[3] = tet[1]; - p[4] = p12; - p[5] = p13; - //Pout<< " 01 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p02; - p[2] = p12; - p[3] = tet[3]; - p[4] = p03; - p[5] = p13; - //Pout<< " 01 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(1, 2)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p13 = planeIntersection(d, tet, 1, 3); - point p02 = planeIntersection(d, tet, 0, 2); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[1]; - p[1] = p01; - p[2] = p13; - p[3] = tet[2]; - p[4] = p02; - p[5] = p23; - //Pout<< " 12 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[3]; - p[1] = p23; - p[2] = p13; - p[3] = tet[0]; - p[4] = p02; - p[5] = p01; - //Pout<< " 12 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(2, 0)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p03 = planeIntersection(d, tet, 0, 3); - point p12 = planeIntersection(d, tet, 1, 2); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[2]; - p[1] = p12; - p[2] = p23; - p[3] = tet[0]; - p[4] = p01; - p[5] = p03; - //Pout<< " 20 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[1]; - p[1] = p12; - p[2] = p01; - p[3] = tet[3]; - p[4] = p23; - p[5] = p03; - //Pout<< " 20 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(0, 3)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p02 = planeIntersection(d, tet, 0, 2); - point p13 = planeIntersection(d, tet, 1, 3); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[3]; - p[1] = p23; - p[2] = p13; - p[3] = tet[0]; - p[4] = p02; - p[5] = p01; - //Pout<< " 03 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p23; - p[2] = p02; - p[3] = tet[1]; - p[4] = p13; - p[5] = p01; - //Pout<< " 03 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(1, 3)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p12 = planeIntersection(d, tet, 1, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[1]; - p[1] = p12; - p[2] = p01; - p[3] = tet[3]; - p[4] = p23; - p[5] = p03; - //Pout<< " 13 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p12; - p[2] = p23; - p[3] = tet[0]; - p[4] = p01; - p[5] = p03; - //Pout<< " 13 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(2, 3)) - { - point p02 = planeIntersection(d, tet, 0, 2); - point p12 = planeIntersection(d, tet, 1, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p13 = planeIntersection(d, tet, 1, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[2]; - p[1] = p02; - p[2] = p12; - p[3] = tet[3]; - p[4] = p03; - p[5] = p13; - //Pout<< " 23 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[0]; - p[1] = p02; - p[2] = p03; - p[3] = tet[1]; - p[4] = p12; - p[5] = p13; - //Pout<< " 23 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else - { - FatalErrorIn("tetSliceWithPlane(..)") << "Missed edge:" << posEdge - << abort(FatalError); - } - } - else if (nPos == 1) - { - // Find the positive tet - label i0 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - i0 = i; - break; - } - } - - label i1 = d.fcIndex(i0); - label i2 = d.fcIndex(i1); - label i3 = d.fcIndex(i2); - - point p01 = planeIntersection(d, tet, i0, i1); - point p02 = planeIntersection(d, tet, i0, i2); - point p03 = planeIntersection(d, tet, i0, i3); - - //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl; - - if (i0 == 0 || i0 == 2) - { - tetPoints t(tet[i0], p01, p02, p03); - //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 1, aboveTets i0==0 or 2"); - aboveOp(t); - - // Prism - FixedList p; - p[0] = tet[i1]; - p[1] = tet[i3]; - p[2] = tet[i2]; - p[3] = p01; - p[4] = p03; - p[5] = p02; - //Pout<< " belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - else - { - tetPoints t(p01, p02, p03, tet[i0]); - //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 1, aboveTets i0==1 or 3"); - aboveOp(t); - - // Prism - FixedList p; - p[0] = tet[i3]; - p[1] = tet[i1]; - p[2] = tet[i2]; - p[3] = p03; - p[4] = p01; - p[5] = p02; - //Pout<< " belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else // nPos == 0 - { - belowOp(tet); - } -} - - void Foam::tetOverlapVolume::tetTetOverlap ( const tetPoints& tetA, @@ -462,16 +53,15 @@ void Foam::tetOverlapVolume::tetTetOverlap label& nInside, FixedList& outsideTets, label& nOutside -) +) const { // Work storage FixedList cutInsideTets; label nCutInside = 0; - storeTetOp inside(insideTets, nInside); - storeTetOp cutInside(cutInsideTets, nCutInside); - dummyTetOp outside; - + tetPointRef::storeOp inside(insideTets, nInside); + tetPointRef::storeOp cutInside(cutInsideTets, nCutInside); + tetPointRef::dummyOp outside; // Cut tetA with all inwards pointing faces of tetB. Any tets remaining @@ -484,7 +74,7 @@ void Foam::tetOverlapVolume::tetTetOverlap // Cut and insert subtets into cutInsideTets (either by getting // an index from freeSlots or by appending to insideTets) or // insert into outsideTets - tetSliceWithPlane(tetA, pl0, cutInside, outside); + tetA.tet().sliceWithPlane(pl0, cutInside, outside); } if (nCutInside == 0) @@ -501,7 +91,7 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nCutInside; i++) { - tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside); + cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside); } if (nInside == 0) @@ -518,7 +108,7 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nInside; i++) { - tetSliceWithPlane(insideTets[i], pl2, cutInside, outside); + insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside); } if (nCutInside == 0) @@ -536,28 +126,27 @@ void Foam::tetOverlapVolume::tetTetOverlap for (label i = 0; i < nCutInside; i++) { - tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside); + cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside); } } } -inline Foam::scalar -Foam::tetOverlapVolume::tetTetOverlapVol +Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol ( const tetPoints& tetA, const tetPoints& tetB -) +) const { FixedList insideTets; label nInside = 0; FixedList cutInsideTets; label nCutInside = 0; - storeTetOp inside(insideTets, nInside); - storeTetOp cutInside(cutInsideTets, nCutInside); - sumTetVolOp volInside; - dummyTetOp outside; + tetPointRef::storeOp inside(insideTets, nInside); + tetPointRef::storeOp cutInside(cutInsideTets, nCutInside); + tetPointRef::sumVolOp volInside; + tetPointRef::dummyOp outside; // face0 plane pl0(tetB[1], tetB[3], tetB[2]); @@ -605,12 +194,12 @@ Foam::tetOverlapVolume::tetTetOverlapVol } -inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb +Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb ( const pointField& points, const face& f, const point& fc -) +) const { treeBoundBox bb(fc, fc); forAll(f, fp) @@ -750,8 +339,8 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp Foam::labelList Foam::tetOverlapVolume::overlappingCells ( - const fvMesh& fromMesh, - const fvMesh& toMesh, + const polyMesh& fromMesh, + const polyMesh& toMesh, const label iTo ) const { @@ -766,30 +355,4 @@ Foam::labelList Foam::tetOverlapVolume::overlappingCells } -/* - - forAll(cellsA, i) - { - label cellAI = cellsA[i]; - treeBoundBox bbA - ( - pointField(meshA.points(), meshA.cellPoints()[cellAI]) - ); - - scalar v = cellCellOverlapVolumeMinDecomp - ( - meshA, - cellAI, - bbA, - meshB, - cellBI, - bbB - ); - - overlapVol += v; - nOverlapTests++; - } - -*/ - // ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H index d0891c14d1..d392958e81 100644 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H +++ b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -36,21 +36,17 @@ SourceFiles #ifndef tetOverlapVolume_H #define tetOverlapVolume_H -#include "tetrahedron.H" -#include "fvMesh.H" -#include "plane.H" -#include "tetPointRef.H" -#include "OFstream.H" -#include "meshTools.H" -#include "indexedOctree.H" -#include "treeDataCell.H" -#include "tetPoints.H" -#include "tetCell.H" -#include "EdgeMap.H" +#include "FixedList.H" +#include "labelList.H" +#include "treeBoundBox.H" namespace Foam { +class primitiveMesh; +class polyMesh; +class tetPoints; + /*---------------------------------------------------------------------------*\ Class tetOverlapVolume Declaration \*---------------------------------------------------------------------------*/ @@ -59,82 +55,6 @@ class tetOverlapVolume { // Private member functions - //- Plane intersection - inline point planeIntersection - ( - const FixedList& d, - const tetPoints& t, - const label negI, - const label posI - ); - - - //- Decompose prism - template inline void decomposePrism - ( - const FixedList& points, - TetOp& op - ); - - - //- Helping cľasses - class dummyTetOp - { - public: - - inline void operator()(const tetPoints&){} - }; - - - class sumTetVolOp - { - public: - scalar vol_; - - inline sumTetVolOp() - : - vol_(0.0) - {} - - inline void operator()(const tetPoints& tet) - { - vol_ += tet.tet().mag(); - } - }; - - - class storeTetOp - { - FixedList& tets_; - label& nTets_; - - public: - - inline storeTetOp(FixedList& tets, label& nTets) - : - tets_(tets), - nTets_(nTets) - {} - - inline void operator()(const tetPoints& tet) - { - tets_[nTets_++] = tet; - } - }; - - - //- Slice. Split tet into subtets above and below plane - template - inline void tetSliceWithPlane - ( - const tetPoints& tet, - const plane& pl, - - AboveTetOp& aboveOp, - BelowTetOp& belowOp - ); - - //- Tet overlap void tetTetOverlap ( @@ -144,31 +64,28 @@ class tetOverlapVolume label& nInside, FixedList& outsideTets, label& nOutside - ); - + ) const; //- Tet Overlap Vol - inline scalar tetTetOverlapVol + scalar tetTetOverlapVol ( const tetPoints& tetA, const tetPoints& tetB - ); - + ) const; //- Return a const treeBoundBox - inline const treeBoundBox pyrBb + treeBoundBox pyrBb ( const pointField& points, const face& f, const point& fc - ); - + ) const; public: //- Runtime type information - TypeName("tetOverlapVolume"); + ClassName("tetOverlapVolume"); // Constructors @@ -177,18 +94,14 @@ public: tetOverlapVolume(); - //- Destructor - virtual ~tetOverlapVolume(); - - // Public members //- Return a list of cells in meshA which overlaps with cellBI in // meshB labelList overlappingCells ( - const fvMesh& meshA, - const fvMesh& meshB, + const polyMesh& meshA, + const polyMesh& meshB, const label cellBI ) const; diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C deleted file mode 100644 index d82cf30a50..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C +++ /dev/null @@ -1,341 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 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 3 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, see . - -Description - Calculation of shape function product for a tetrahedron - -\*---------------------------------------------------------------------------*/ - -#include "tetrahedron.H" -#include "triPointRef.H" -#include "scalarField.H" - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -// (Probably very inefficient) minimum containment sphere calculation. -// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf: -// Sphere ctr is smallest one of -// - tet circumcentre -// - triangle circumcentre -// - edge mids -template -Foam::pointHit Foam::tetrahedron::containmentSphere -( - const scalar tol -) const -{ - const scalar fac = 1 + tol; - - // Halve order of tolerance for comparisons of sqr. - const scalar facSqr = Foam::sqrt(fac); - - - // 1. Circumcentre itself. - - pointHit pHit(circumCentre()); - pHit.setHit(); - scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_); - - - // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and - // check if 4th point is inside. - - { - point ctr = triPointRef(a_, b_, c_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(a_, b_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(a_, c_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - a_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - { - point ctr = triPointRef(b_, c_, d_).circumCentre(); - scalar radiusSqr = magSqr(ctr - b_); - - if - ( - radiusSqr < minRadiusSqr - && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - - // 3. Try midpoints of edges - - // mid of edge A-B - { - point ctr = 0.5*(a_ + b_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(c_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge A-C - { - point ctr = 0.5*(a_ + c_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(b_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge A-D - { - point ctr = 0.5*(a_ + d_); - scalar radiusSqr = magSqr(a_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(b_-ctr) <= testRadSrq - && magSqr(c_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge B-C - { - point ctr = 0.5*(b_ + c_); - scalar radiusSqr = magSqr(b_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(d_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge B-D - { - point ctr = 0.5*(b_ + d_); - scalar radiusSqr = magSqr(b_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(c_-ctr) <= testRadSrq) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - // mid of edge C-D - { - point ctr = 0.5*(c_ + d_); - scalar radiusSqr = magSqr(c_ - ctr); - scalar testRadSrq = facSqr*radiusSqr; - - if - ( - radiusSqr < minRadiusSqr - && magSqr(a_-ctr) <= testRadSrq - && magSqr(b_-ctr) <= testRadSrq - ) - { - pHit.setMiss(false); - pHit.setPoint(ctr); - minRadiusSqr = radiusSqr; - } - } - - - pHit.setDistance(sqrt(minRadiusSqr)); - - return pHit; -} - - -template -void Foam::tetrahedron::gradNiSquared -( - scalarField& buffer -) const -{ - // Change of sign between face area vector and gradient - // does not matter because of square - - // Warning: Added a mag to produce positive coefficients even if - // the tetrahedron is twisted inside out. This is pretty - // dangerous, but essential for mesh motion. - scalar magVol = Foam::mag(mag()); - - buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol; - buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol; - buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol; - buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol; -} - - -template -void Foam::tetrahedron::gradNiDotGradNj -( - scalarField& buffer -) const -{ - // Warning. Ordering of edges needs to be the same for a tetrahedron - // class, a tetrahedron cell shape model and a tetCell - - // Warning: Added a mag to produce positive coefficients even if - // the tetrahedron is twisted inside out. This is pretty - // dangerous, but essential for mesh motion. - - // Double change of sign between face area vector and gradient - - scalar magVol = Foam::mag(mag()); - vector sa = Sa(); - vector sb = Sb(); - vector sc = Sc(); - vector sd = Sd(); - - buffer[0] = (1.0/9.0)*(sa & sb)/magVol; - buffer[1] = (1.0/9.0)*(sa & sc)/magVol; - buffer[2] = (1.0/9.0)*(sa & sd)/magVol; - buffer[3] = (1.0/9.0)*(sd & sb)/magVol; - buffer[4] = (1.0/9.0)*(sb & sc)/magVol; - buffer[5] = (1.0/9.0)*(sd & sc)/magVol; -} - - -template -void Foam::tetrahedron::gradNiGradNi -( - tensorField& buffer -) const -{ - // Change of sign between face area vector and gradient - // does not matter because of square - - scalar magVol = Foam::mag(mag()); - - buffer[0] = (1.0/9.0)*sqr(Sa())/magVol; - buffer[1] = (1.0/9.0)*sqr(Sb())/magVol; - buffer[2] = (1.0/9.0)*sqr(Sc())/magVol; - buffer[3] = (1.0/9.0)*sqr(Sd())/magVol; -} - - -template -void Foam::tetrahedron::gradNiGradNj -( - tensorField& buffer -) const -{ - // Warning. Ordering of edges needs to be the same for a tetrahedron - // class, a tetrahedron cell shape model and a tetCell - - // Double change of sign between face area vector and gradient - - scalar magVol = Foam::mag(mag()); - vector sa = Sa(); - vector sb = Sb(); - vector sc = Sc(); - vector sd = Sd(); - - buffer[0] = (1.0/9.0)*(sa * sb)/magVol; - buffer[1] = (1.0/9.0)*(sa * sc)/magVol; - buffer[2] = (1.0/9.0)*(sa * sd)/magVol; - buffer[3] = (1.0/9.0)*(sd * sb)/magVol; - buffer[4] = (1.0/9.0)*(sb * sc)/magVol; - buffer[5] = (1.0/9.0)*(sd * sc)/magVol; -} - - -// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H deleted file mode 100644 index 13e7acd4eb..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H +++ /dev/null @@ -1,315 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 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 3 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, see . - -Class - Foam::tetrahedron - -Description - A tetrahedron primitive. - - Ordering of edges needs to be the same for a tetrahedron - class, a tetrahedron cell shape model and a tetCell. - -SourceFiles - tetrahedronI.H - tetrahedron.C - -\*---------------------------------------------------------------------------*/ - -#ifndef tetrahedron_H -#define tetrahedron_H - -#include "point.H" -#include "primitiveFieldsFwd.H" -#include "pointHit.H" -#include "cachedRandom.H" -#include "Random.H" -#include "FixedList.H" -#include "UList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -class Istream; -class Ostream; -class tetPoints; -class plane; - -// Forward declaration of friend functions and operators - -template class tetrahedron; - -template -inline Istream& operator>> -( - Istream&, - tetrahedron& -); - -template -inline Ostream& operator<< -( - Ostream&, - const tetrahedron& -); - - -/*---------------------------------------------------------------------------*\ - class tetrahedron Declaration -\*---------------------------------------------------------------------------*/ - -template -class tetrahedron -{ -public: - - // Classes for use in sliceWithPlane. What to do with decomposition - // of tet. - - //- Dummy - class dummyOp - { - public: - inline void operator()(const tetPoints&); - }; - - //- Sum resulting volumes - class sumVolOp - { - public: - scalar vol_; - - inline sumVolOp(); - - inline void operator()(const tetPoints&); - }; - - //- Store resulting tets - class storeOp - { - FixedList& tets_; - label& nTets_; - - public: - inline storeOp(FixedList&, label&); - - inline void operator()(const tetPoints&); - }; - -private: - - // Private data - - PointRef a_, b_, c_, d_; - - inline static point planeIntersection - ( - const FixedList&, - const tetPoints&, - const label, - const label - ); - - template - inline static void decomposePrism - ( - const FixedList& points, - TetOp& op - ); - - template - inline static void tetSliceWithPlane - ( - const plane& pl, - const tetPoints& tet, - AboveTetOp& aboveOp, - BelowTetOp& belowOp - ); - - -public: - - // Member constants - - enum - { - nVertices = 4, // Number of vertices in tetrahedron - nEdges = 6 // Number of edges in tetrahedron - }; - - - // Constructors - - //- Construct from points - inline tetrahedron - ( - const Point& a, - const Point& b, - const Point& c, - const Point& d - ); - - //- Construct from four points in the list of points - inline tetrahedron - ( - const UList&, - const FixedList& indices - ); - - //- Construct from Istream - inline tetrahedron(Istream&); - - - // Member Functions - - // Access - - //- Return vertices - inline const Point& a() const; - - inline const Point& b() const; - - inline const Point& c() const; - - inline const Point& d() const; - - - // Properties - - //- Return face normal - inline vector Sa() const; - - inline vector Sb() const; - - inline vector Sc() const; - - inline vector Sd() const; - - //- Return centre (centroid) - inline Point centre() const; - - //- Return volume - inline scalar mag() const; - - //- Return circum-centre - inline Point circumCentre() const; - - //- Return circum-radius - inline scalar circumRadius() const; - - //- Return quality: Ratio of tetrahedron and circum-sphere - // volume, scaled so that a regular tetrahedron has a - // quality of 1 - inline scalar quality() const; - - //- Return a random point in the tetrahedron from a - // uniform distribution - inline Point randomPoint(Random& rndGen) const; - - //- Return a random point in the tetrahedron from a - // uniform distribution - inline Point randomPoint(cachedRandom& rndGen) const; - - //- Calculate the barycentric coordinates of the given - // point, in the same order as a, b, c, d. Returns the - // determinant of the solution. - inline scalar barycentric - ( - const point& pt, - List& bary - ) const; - - //- Return nearest point to p on tetrahedron - inline pointHit nearestPoint(const point& p) const; - - //- Return true if point is inside tetrahedron - inline bool inside(const point& pt) const; - - //- Decompose tet into tets above and below plane - template - inline void sliceWithPlane - ( - const plane& pl, - AboveTetOp& aboveOp, - BelowTetOp& belowOp - ) const; - - - //- Return (min)containment sphere, i.e. the smallest sphere with - // all points inside. Returns pointHit with: - // - hit : if sphere is equal to circumsphere - // (biggest sphere) - // - point : centre of sphere - // - distance : radius of sphere - // - eligiblemiss: false - // Tol (small compared to 1, e.g. 1E-9) is used to determine - // whether point is inside: mag(pt - ctr) < (1+tol)*radius. - pointHit containmentSphere(const scalar tol) const; - - //- Fill buffer with shape function products - void gradNiSquared(scalarField& buffer) const; - - void gradNiDotGradNj(scalarField& buffer) const; - - void gradNiGradNi(tensorField& buffer) const; - - void gradNiGradNj(tensorField& buffer) const; - - - // IOstream operators - - friend Istream& operator>> - ( - Istream&, - tetrahedron& - ); - - friend Ostream& operator<< - ( - Ostream&, - const tetrahedron& - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "tetrahedronI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "tetrahedron.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H b/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H deleted file mode 100644 index 922ae0057e..0000000000 --- a/src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedronI.H +++ /dev/null @@ -1,1012 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2012-2012 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 3 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, see . - -\*---------------------------------------------------------------------------*/ - -#include "triangle.H" -#include "IOstreams.H" -#include "triPointRef.H" -#include "tetPoints.H" -#include "plane.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -inline Foam::tetrahedron::tetrahedron -( - const Point& a, - const Point& b, - const Point& c, - const Point& d -) -: - a_(a), - b_(b), - c_(c), - d_(d) -{} - - -template -inline Foam::tetrahedron::tetrahedron -( - const UList& points, - const FixedList& indices -) -: - a_(points[indices[0]]), - b_(points[indices[1]]), - c_(points[indices[2]]), - d_(points[indices[3]]) -{} - - -template -inline Foam::tetrahedron::tetrahedron(Istream& is) -{ - is >> *this; -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -inline const Point& Foam::tetrahedron::a() const -{ - return a_; -} - - -template -inline const Point& Foam::tetrahedron::b() const -{ - return b_; -} - - -template -inline const Point& Foam::tetrahedron::c() const -{ - return c_; -} - - -template -inline const Point& Foam::tetrahedron::d() const -{ - return d_; -} - - -template -inline Foam::vector Foam::tetrahedron::Sa() const -{ - return triangle(b_, c_, d_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sb() const -{ - return triangle(a_, d_, c_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sc() const -{ - return triangle(a_, b_, d_).normal(); -} - - -template -inline Foam::vector Foam::tetrahedron::Sd() const -{ - return triangle(a_, c_, b_).normal(); -} - - -template -inline Point Foam::tetrahedron::centre() const -{ - return 0.25*(a_ + b_ + c_ + d_); -} - - -template -inline Foam::scalar Foam::tetrahedron::mag() const -{ - return (1.0/6.0)*(((b_ - a_) ^ (c_ - a_)) & (d_ - a_)); -} - - -template -inline Point Foam::tetrahedron::circumCentre() const -{ - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; - - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); - - vector ba = b ^ a; - vector ca = c ^ a; - - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); - - if (Foam::mag(denom) < ROOTVSMALL) - { - // Degenerate tetrahedron, returning centre instead of circumCentre. - - return centre(); - } - - return a_ + 0.5*(a + num/denom); -} - - -template -inline Foam::scalar Foam::tetrahedron::circumRadius() const -{ - vector a = b_ - a_; - vector b = c_ - a_; - vector c = d_ - a_; - - scalar lambda = magSqr(c) - (a & c); - scalar mu = magSqr(b) - (a & b); - - vector ba = b ^ a; - vector ca = c ^ a; - - vector num = lambda*ba - mu*ca; - scalar denom = (c & ba); - - if (Foam::mag(denom) < ROOTVSMALL) - { - // Degenerate tetrahedron, returning GREAT for circumRadius. - return GREAT; - } - - return Foam::mag(0.5*(a + num/denom)); -} - - -template -inline Foam::scalar Foam::tetrahedron::quality() const -{ - return - mag() - /( - 8.0/(9.0*sqrt(3.0)) - *pow3(min(circumRadius(), GREAT)) - + ROOTVSMALL - ); -} - - -template -inline Point Foam::tetrahedron::randomPoint -( - Random& rndGen -) const -{ - // Adapted from - // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html - - scalar s = rndGen.scalar01(); - scalar t = rndGen.scalar01(); - scalar u = rndGen.scalar01(); - - if (s + t > 1.0) - { - s = 1.0 - s; - t = 1.0 - t; - } - - if (t + u > 1.0) - { - scalar tmp = u; - u = 1.0 - s - t; - t = 1.0 - tmp; - } - else if (s + t + u > 1.0) - { - scalar tmp = u; - u = s + t + u - 1.0; - s = 1.0 - t - tmp; - } - - return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_; -} - - -template -inline Point Foam::tetrahedron::randomPoint -( - cachedRandom& rndGen -) const -{ - // Adapted from - // http://vcg.isti.cnr.it/activities/geometryegraphics/pointintetraedro.html - - scalar s = rndGen.sample01(); - scalar t = rndGen.sample01(); - scalar u = rndGen.sample01(); - - if (s + t > 1.0) - { - s = 1.0 - s; - t = 1.0 - t; - } - - if (t + u > 1.0) - { - scalar tmp = u; - u = 1.0 - s - t; - t = 1.0 - tmp; - } - else if (s + t + u > 1.0) - { - scalar tmp = u; - u = s + t + u - 1.0; - s = 1.0 - t - tmp; - } - - return (1 - s - t - u)*a_ + s*b_ + t*c_ + u*d_; -} - - -template -Foam::scalar Foam::tetrahedron::barycentric -( - const point& pt, - List& bary -) const -{ - // From: - // http://en.wikipedia.org/wiki/Barycentric_coordinate_system_(mathematics) - - vector e0(a_ - d_); - vector e1(b_ - d_); - vector e2(c_ - d_); - - tensor t - ( - e0.x(), e1.x(), e2.x(), - e0.y(), e1.y(), e2.y(), - e0.z(), e1.z(), e2.z() - ); - - scalar detT = det(t); - - if (Foam::mag(detT) < SMALL) - { - // Degenerate tetrahedron, returning 1/4 barycentric coordinates. - - bary = List(4, 0.25); - - return detT; - } - - vector res = inv(t, detT) & (pt - d_); - - bary.setSize(4); - - bary[0] = res.x(); - bary[1] = res.y(); - bary[2] = res.z(); - bary[3] = (1.0 - res.x() - res.y() - res.z()); - - return detT; -} - - -template -inline Foam::pointHit Foam::tetrahedron::nearestPoint -( - const point& p -) const -{ - // Adapted from: - // Real-time collision detection, Christer Ericson, 2005, p142-144 - - // Assuming initially that the point is inside all of the - // halfspaces, so closest to itself. - - point closestPt = p; - - scalar minOutsideDistance = VGREAT; - - bool inside = true; - - if (((p - b_) & Sa()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(b_, c_, d_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sb()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, d_, c_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sc()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, b_, d_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - if (((p - a_) & Sd()) >= 0) - { - // p is outside halfspace plane of tri - pointHit info = triangle(a_, c_, b_).nearestPoint(p); - - inside = false; - - if (info.distance() < minOutsideDistance) - { - closestPt = info.rawPoint(); - - minOutsideDistance = info.distance(); - } - } - - // If the point is inside, then the distance to the closest point - // is zero - if (inside) - { - minOutsideDistance = 0; - } - - return pointHit - ( - inside, - closestPt, - minOutsideDistance, - !inside - ); -} - - -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; -} - - -template -inline void Foam::tetrahedron::dummyOp::operator() -( - const tetPoints& -) -{} - - -template -inline Foam::tetrahedron::sumVolOp::sumVolOp() -: - vol_(0.0) -{} - - -template -inline void Foam::tetrahedron::sumVolOp::operator() -( - const tetPoints& tet -) -{ - vol_ += tet.tet().mag(); -} - - -template -inline Foam::tetrahedron::storeOp::storeOp -( - FixedList& tets, - label& nTets -) -: - tets_(tets), - nTets_(nTets) -{} - - -template -inline void Foam::tetrahedron::storeOp::operator() -( - const tetPoints& tet -) -{ - tets_[nTets_++] = tet; -} - - -template -inline Foam::point Foam::tetrahedron::planeIntersection -( - const FixedList& d, - const tetPoints& t, - const label negI, - const label posI -) -{ - return - (d[posI]*t[negI] - d[negI]*t[posI]) - / (-d[negI]+d[posI]); -} - - -template -template -inline void Foam::tetrahedron::decomposePrism -( - const FixedList& points, - TetOp& op -) -{ - op(tetPoints(points[1], points[3], points[2], points[0])); - op(tetPoints(points[1], points[2], points[3], points[4])); - op(tetPoints(points[4], points[2], points[3], points[5])); -} - - -template -template -inline void Foam::tetrahedron:: -tetSliceWithPlane -( - const plane& pl, - const tetPoints& tet, - AboveTetOp& aboveOp, - BelowTetOp& belowOp -) -{ - // Distance to plane - FixedList d; - label nPos = 0; - forAll(tet, i) - { - d[i] = ((tet[i]-pl.refPoint()) & pl.normal()); - if (d[i] > 0) - { - nPos++; - } - } - - if (nPos == 4) - { - aboveOp(tet); - } - else if (nPos == 3) - { - // Sliced into below tet and above prism. Prism gets split into - // two tets. - - // Find the below tet - label i0 = -1; - forAll(d, i) - { - if (d[i] <= 0) - { - i0 = i; - break; - } - } - - label i1 = d.fcIndex(i0); - label i2 = d.fcIndex(i1); - label i3 = d.fcIndex(i2); - - point p01 = planeIntersection(d, tet, i0, i1); - point p02 = planeIntersection(d, tet, i0, i2); - point p03 = planeIntersection(d, tet, i0, i3); - - // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad - // ,, 1 : ,, inwards pointing triad - // ,, 2 : ,, outwards pointing triad - // ,, 3 : ,, inwards pointing triad - - //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl; - - if (i0 == 0 || i0 == 2) - { - tetPoints t(tet[i0], p01, p02, p03); - //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 3, belowTet i0==0 or 2"); - belowOp(t); - - // Prism - FixedList p; - p[0] = tet[i1]; - p[1] = tet[i3]; - p[2] = tet[i2]; - p[3] = p01; - p[4] = p03; - p[5] = p02; - //Pout<< " aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - else - { - tetPoints t(p01, p02, p03, tet[i0]); - //Pout<< " belowtet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 3, belowTet i0==1 or 3"); - belowOp(t); - - // Prism - FixedList p; - p[0] = tet[i3]; - p[1] = tet[i1]; - p[2] = tet[i2]; - p[3] = p03; - p[4] = p01; - p[5] = p02; - //Pout<< " aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - } - else if (nPos == 2) - { - // Tet cut into two prisms. Determine the positive one. - label pos0 = -1; - label pos1 = -1; - label neg0 = -1; - label neg1 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - if (pos0 == -1) - { - pos0 = i; - } - else - { - pos1 = i; - } - } - else - { - if (neg0 == -1) - { - neg0 = i; - } - else - { - neg1 = i; - } - } - } - - //Pout<< "Split 2pos tet " << tet << " d:" << d - // << " around pos0:" << pos0 << " pos1:" << pos1 - // << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl; - - const edge posEdge(pos0, pos1); - - if (posEdge == edge(0, 1)) - { - point p02 = planeIntersection(d, tet, 0, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p12 = planeIntersection(d, tet, 1, 2); - point p13 = planeIntersection(d, tet, 1, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[0]; - p[1] = p02; - p[2] = p03; - p[3] = tet[1]; - p[4] = p12; - p[5] = p13; - //Pout<< " 01 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p02; - p[2] = p12; - p[3] = tet[3]; - p[4] = p03; - p[5] = p13; - //Pout<< " 01 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(1, 2)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p13 = planeIntersection(d, tet, 1, 3); - point p02 = planeIntersection(d, tet, 0, 2); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[1]; - p[1] = p01; - p[2] = p13; - p[3] = tet[2]; - p[4] = p02; - p[5] = p23; - //Pout<< " 12 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[3]; - p[1] = p23; - p[2] = p13; - p[3] = tet[0]; - p[4] = p02; - p[5] = p01; - //Pout<< " 12 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(2, 0)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p03 = planeIntersection(d, tet, 0, 3); - point p12 = planeIntersection(d, tet, 1, 2); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[2]; - p[1] = p12; - p[2] = p23; - p[3] = tet[0]; - p[4] = p01; - p[5] = p03; - //Pout<< " 20 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[1]; - p[1] = p12; - p[2] = p01; - p[3] = tet[3]; - p[4] = p23; - p[5] = p03; - //Pout<< " 20 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(0, 3)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p02 = planeIntersection(d, tet, 0, 2); - point p13 = planeIntersection(d, tet, 1, 3); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[3]; - p[1] = p23; - p[2] = p13; - p[3] = tet[0]; - p[4] = p02; - p[5] = p01; - //Pout<< " 03 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p23; - p[2] = p02; - p[3] = tet[1]; - p[4] = p13; - p[5] = p01; - //Pout<< " 03 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(1, 3)) - { - point p01 = planeIntersection(d, tet, 0, 1); - point p12 = planeIntersection(d, tet, 1, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p23 = planeIntersection(d, tet, 2, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[1]; - p[1] = p12; - p[2] = p01; - p[3] = tet[3]; - p[4] = p23; - p[5] = p03; - //Pout<< " 13 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[2]; - p[1] = p12; - p[2] = p23; - p[3] = tet[0]; - p[4] = p01; - p[5] = p03; - //Pout<< " 13 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else if (posEdge == edge(2, 3)) - { - point p02 = planeIntersection(d, tet, 0, 2); - point p12 = planeIntersection(d, tet, 1, 2); - point p03 = planeIntersection(d, tet, 0, 3); - point p13 = planeIntersection(d, tet, 1, 3); - // Split the resulting prism - { - FixedList p; - p[0] = tet[2]; - p[1] = p02; - p[2] = p12; - p[3] = tet[3]; - p[4] = p03; - p[5] = p13; - //Pout<< " 23 aboveprism:" << p << endl; - decomposePrism(p, aboveOp); - } - { - FixedList p; - p[0] = tet[0]; - p[1] = p02; - p[2] = p03; - p[3] = tet[1]; - p[4] = p12; - p[5] = p13; - //Pout<< " 23 belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else - { - FatalErrorIn("tetSliceWithPlane(..)") - << "Missed edge:" << posEdge - << abort(FatalError); - } - } - else if (nPos == 1) - { - // Find the positive tet - label i0 = -1; - forAll(d, i) - { - if (d[i] > 0) - { - i0 = i; - break; - } - } - - label i1 = d.fcIndex(i0); - label i2 = d.fcIndex(i1); - label i3 = d.fcIndex(i2); - - point p01 = planeIntersection(d, tet, i0, i1); - point p02 = planeIntersection(d, tet, i0, i2); - point p03 = planeIntersection(d, tet, i0, i3); - - //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl; - - if (i0 == 0 || i0 == 2) - { - tetPoints t(tet[i0], p01, p02, p03); - //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 1, aboveTets i0==0 or 2"); - aboveOp(t); - - // Prism - FixedList p; - p[0] = tet[i1]; - p[1] = tet[i3]; - p[2] = tet[i2]; - p[3] = p01; - p[4] = p03; - p[5] = p02; - //Pout<< " belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - else - { - tetPoints t(p01, p02, p03, tet[i0]); - //Pout<< " abovetet:" << t << " around i0:" << i0 << endl; - //checkTet(t, "nPos 1, aboveTets i0==1 or 3"); - aboveOp(t); - - // Prism - FixedList p; - p[0] = tet[i3]; - p[1] = tet[i1]; - p[2] = tet[i2]; - p[3] = p03; - p[4] = p01; - p[5] = p02; - //Pout<< " belowprism:" << p << endl; - decomposePrism(p, belowOp); - } - } - else // nPos == 0 - { - belowOp(tet); - } -} - - -template -template -inline void Foam::tetrahedron::sliceWithPlane -( - const plane& pl, - AboveTetOp& aboveOp, - BelowTetOp& belowOp -) const -{ - tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp); -} - - -// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // - -template -inline Foam::Istream& Foam::operator>> -( - Istream& is, - tetrahedron& t -) -{ - is.readBegin("tetrahedron"); - is >> t.a_ >> t.b_ >> t.c_ >> t.d_; - is.readEnd("tetrahedron"); - - is.check("Istream& operator>>(Istream&, tetrahedron&)"); - - return is; -} - - -template -inline Foam::Ostream& Foam::operator<< -( - Ostream& os, - const tetrahedron& t -) -{ - os << nl - << token::BEGIN_LIST - << t.a_ << token::SPACE - << t.b_ << token::SPACE - << t.c_ << token::SPACE - << t.d_ - << token::END_LIST; - - return os; -} - - -// ************************************************************************* // diff --git a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSolution b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSolution index 8ae8ee05ae..97d24d6130 100644 --- a/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSolution +++ b/tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSolution @@ -68,11 +68,11 @@ PIMPLE cAlpha 2; } -relaxation +relaxationFactors { - U 1; - k 1; - epsilon 1; + "U.*" 1; + "k.*" 1; + "epsilon.*" 1; } // ************************************************************************* //