diff --git a/applications/solvers/compressible/sonicFoam/UEqn.H b/applications/solvers/compressible/sonicFoam/UEqn.H new file mode 100644 index 0000000000..25506783ae --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/UEqn.H @@ -0,0 +1,8 @@ +fvVectorMatrix UEqn +( + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) +); + +solve(UEqn == -fvc::grad(p)); diff --git a/applications/solvers/compressible/sonicFoam/compressibleContinuityErrs.H b/applications/solvers/compressible/sonicFoam/compressibleContinuityErrs.H deleted file mode 100644 index 128d99c946..0000000000 --- a/applications/solvers/compressible/sonicFoam/compressibleContinuityErrs.H +++ /dev/null @@ -1,12 +0,0 @@ -{ -# include "rhoEqn.H" -} -{ - scalar sumLocalContErr = (sum(mag(rho - psi*p))/sum(rho)).value(); - scalar globalContErr = (sum(rho - psi*p)/sum(rho)).value(); - cumulativeContErr += globalContErr; - - Info<< "time step continuity errors : sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr << endl; -} diff --git a/applications/solvers/compressible/sonicFoam/createFields.H b/applications/solvers/compressible/sonicFoam/createFields.H index bbb5d10526..5d03dd2bb7 100644 --- a/applications/solvers/compressible/sonicFoam/createFields.H +++ b/applications/solvers/compressible/sonicFoam/createFields.H @@ -7,7 +7,7 @@ basicPsiThermo& thermo = pThermo(); volScalarField& p = thermo.p(); - volScalarField& h = thermo.h(); + volScalarField& e = thermo.e(); const volScalarField& psi = thermo.psi(); volScalarField rho @@ -35,7 +35,7 @@ mesh ); -# include "compressibleCreatePhi.H" + #include "compressibleCreatePhi.H" Info<< "Creating turbulence model\n" << endl; @@ -49,7 +49,3 @@ thermo ) ); - - Info<< "Creating field DpDt\n" << endl; - volScalarField DpDt = - fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); diff --git a/applications/solvers/compressible/sonicFoam/eEqn.H b/applications/solvers/compressible/sonicFoam/eEqn.H new file mode 100644 index 0000000000..1d1d6aef0b --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/eEqn.H @@ -0,0 +1,12 @@ +{ + solve + ( + fvm::ddt(rho, e) + + fvm::div(phi, e) + - fvm::laplacian(turbulence->alphaEff(), e) + == + - p*fvc::div(phi/fvc::interpolate(rho)) + ); + + thermo.correct(); +} diff --git a/applications/solvers/compressible/sonicFoam/hEqn.H b/applications/solvers/compressible/sonicFoam/hEqn.H deleted file mode 100644 index baa2dab343..0000000000 --- a/applications/solvers/compressible/sonicFoam/hEqn.H +++ /dev/null @@ -1,12 +0,0 @@ -{ - solve - ( - fvm::ddt(rho, h) - + fvm::div(phi, h) - - fvm::laplacian(turbulence->alphaEff(), h) - == - DpDt - ); - - thermo.correct(); -} diff --git a/applications/solvers/compressible/sonicFoam/pEqn.H b/applications/solvers/compressible/sonicFoam/pEqn.H new file mode 100644 index 0000000000..96a500d4c2 --- /dev/null +++ b/applications/solvers/compressible/sonicFoam/pEqn.H @@ -0,0 +1,37 @@ +rho = thermo.rho(); + +volScalarField rUA = 1.0/UEqn.A(); +U = rUA*UEqn.H(); + +surfaceScalarField phid +( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ) +); + +for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) +{ + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho*rUA, p) + ); + + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi = pEqn.flux(); + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +U -= rUA*fvc::grad(p); +U.correctBoundaryConditions(); diff --git a/applications/solvers/compressible/sonicFoam/readThermodynamicProperties.H b/applications/solvers/compressible/sonicFoam/readThermodynamicProperties.H deleted file mode 100644 index 1fc57fc5fd..0000000000 --- a/applications/solvers/compressible/sonicFoam/readThermodynamicProperties.H +++ /dev/null @@ -1,23 +0,0 @@ - Info<< "Reading thermodynamicProperties\n" << endl; - - IOdictionary thermodynamicProperties - ( - IOobject - ( - "thermodynamicProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - - dimensionedScalar R - ( - thermodynamicProperties.lookup("R") - ); - - dimensionedScalar Cv - ( - thermodynamicProperties.lookup("Cv") - ); diff --git a/applications/solvers/compressible/sonicFoam/readTransportProperties.H b/applications/solvers/compressible/sonicFoam/readTransportProperties.H deleted file mode 100644 index 1502e2033a..0000000000 --- a/applications/solvers/compressible/sonicFoam/readTransportProperties.H +++ /dev/null @@ -1,18 +0,0 @@ - Info<< "Reading transportProperties\n" << endl; - - IOdictionary transportProperties - ( - IOobject - ( - "transportProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - - dimensionedScalar mu - ( - transportProperties.lookup("mu") - ); diff --git a/applications/solvers/compressible/sonicFoam/sonicFoam.C b/applications/solvers/compressible/sonicFoam/sonicFoam.C index c9eda00fa3..df100262ef 100644 --- a/applications/solvers/compressible/sonicFoam/sonicFoam.C +++ b/applications/solvers/compressible/sonicFoam/sonicFoam.C @@ -58,64 +58,21 @@ int main(int argc, char *argv[]) #include "rhoEqn.H" - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - + turbulence->divDevRhoReff(U) - ); + #include "UEqn.H" - solve(UEqn == -fvc::grad(p)); - - #include "hEqn.H" + #include "eEqn.H" // --- PISO loop for (int corr=0; corrcorrect(); - rho = psi*p; + rho = thermo.rho(); runTime.write(); diff --git a/applications/solvers/incompressible/channelFoam/readTransportProperties.H b/applications/solvers/incompressible/channelFoam/readTransportProperties.H index c99d0ae4b5..87c93c9195 100644 --- a/applications/solvers/incompressible/channelFoam/readTransportProperties.H +++ b/applications/solvers/incompressible/channelFoam/readTransportProperties.H @@ -7,7 +7,8 @@ runTime.constant(), mesh, IOobject::MUST_READ, - IOobject::NO_WRITE + IOobject::NO_WRITE, + false ) ); diff --git a/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H b/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H index 3347db180b..aeb3137adf 100644 --- a/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H +++ b/applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H @@ -1,7 +1,7 @@ - Info<< "Creating merge patch pairs" << nl << endl; - if (mergePatchPairs.size()) { + Info<< "Creating merge patch pairs" << nl << endl; + // Create and add point and face zones and mesh modifiers List pz(mergePatchPairs.size()); List fz(3*mergePatchPairs.size()); diff --git a/applications/utilities/mesh/manipulation/cellSet/cellSetDict b/applications/utilities/mesh/manipulation/cellSet/cellSetDict index 95465d8f54..fd11ab871d 100644 --- a/applications/utilities/mesh/manipulation/cellSet/cellSetDict +++ b/applications/utilities/mesh/manipulation/cellSet/cellSetDict @@ -26,12 +26,31 @@ action new; topoSetSources ( + // Select by explicitly providing cell labels + labelToCell + { + value (12 13 56); // labels of cells + } + // Copy elements from cellSet cellToCell { set c1; } + // Cells in cell zone + zoneToCell + { + name ".*Zone"; // Name of cellZone, regular expressions allowed + } + + // Cells on master or slave side of faceZone + faceZoneToCell + { + name ".*Zone"; // Name of faceZone, regular expressions allowed + option master; // master/slave + } + // Select based on faceSet faceToCell { @@ -51,12 +70,6 @@ topoSetSources //option all; // cell with all points in pointSet } - // Select by explicitly providing cell labels - labelToCell - { - value (12 13 56); // labels of cells - } - // Select based on cellShape shapeToCell { @@ -87,7 +100,6 @@ topoSetSources radius 5.0; } - // Cells with centre within sphere sphereToCell { @@ -95,20 +107,6 @@ topoSetSources radius 5.0; } - // Cells in cell zone - zoneToCell - { - name ".*Zone"; // Name of cellZone, regular expressions allowed - } - - // values of field within certain range - fieldToCell - { - fieldName U; // Note: uses mag(U) since volVectorField - min 0.1; - max 0.5; - } - // Cells with cellCentre nearest to coordinates nearestToCell { @@ -129,6 +127,22 @@ topoSetSources // and near surf curvature // (set to -100 if not used) } + + // values of field within certain range + fieldToCell + { + fieldName U; // Note: uses mag(U) since volVectorField + min 0.1; + max 0.5; + } + + // Mesh region (non-face connected part of (subset of)mesh) + regionToCell + { + set c0; // name of cellSet giving mesh subset + insidePoint (1 2 3); // point inside region to select + } + ); diff --git a/applications/utilities/mesh/manipulation/pointSet/pointSetDict b/applications/utilities/mesh/manipulation/pointSet/pointSetDict index 1c704a1868..895dcd7f0a 100644 --- a/applications/utilities/mesh/manipulation/pointSet/pointSetDict +++ b/applications/utilities/mesh/manipulation/pointSet/pointSetDict @@ -51,18 +51,24 @@ topoSetSources value (12 13 56); // labels of points } - // Points with coordinate within box - boxToPoint - { - box (0 0 0) (1 1 1); - } - // All points in pointzone zoneToPoint { name ".*Zone"; // name of pointZone, regular expressions allowed } + // Points nearest to coordinates + nearestToPoint + { + points ((0 0 0) (1 1 1)); + } + + // Points with coordinate within box + boxToPoint + { + box (0 0 0) (1 1 1); + } + // Select based on surface surfaceToPoint { diff --git a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C index dfd7262fd1..dcb744e801 100644 --- a/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C +++ b/applications/utilities/mesh/manipulation/setsToZones/setsToZones.C @@ -28,8 +28,8 @@ Description There is one catch: for faceZones you also need to specify a flip condition which basically denotes the side of the face. In this app - it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) and - any face whose neighbour is in the cellSet gets a flip=true. + it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) which + is the masterCells of the zone. There are lots of situations in which this will go wrong but it is the best I can think of for now. @@ -201,13 +201,10 @@ int main(int argc, char *argv[]) if (!noFlipMap) { - word setName(set.name() + "Cells"); + word setName(set.name() + "SlaveCells"); Pout<< "Trying to load cellSet " << setName - << " to find out the flipMap." << nl - << "If the neighbour side of the face is in the cellSet" - << " the flipMap becomes true," << nl - << "in all other cases it stays false." + << " to find out the slave side of the zone." << nl << " If you do not care about the flipMap" << " (i.e. do not use the sideness)" << nl << "use the -noFlipMap command line option." @@ -230,7 +227,7 @@ int main(int argc, char *argv[]) && !cells.found(mesh.faceNeighbour()[faceI]) ) { - flip = false; + flip = true; } else if ( @@ -238,7 +235,7 @@ int main(int argc, char *argv[]) && cells.found(mesh.faceNeighbour()[faceI]) ) { - flip = true; + flip = false; } else { @@ -260,11 +257,11 @@ int main(int argc, char *argv[]) { if (cells.found(mesh.faceOwner()[faceI])) { - flip = false; + flip = true; } else { - flip = true; + flip = false; } } diff --git a/src/OpenFOAM/containers/Lists/List/List.C b/src/OpenFOAM/containers/Lists/List/List.C index 3fa3dcc874..08422fd9ba 100644 --- a/src/OpenFOAM/containers/Lists/List/List.C +++ b/src/OpenFOAM/containers/Lists/List/List.C @@ -35,8 +35,6 @@ License #include "BiIndirectList.H" #include "contiguous.H" -#include - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * // @@ -442,34 +440,6 @@ void Foam::List::transfer(SortableList& a) } -template -void Foam::sort(List& a) -{ - std::sort(a.begin(), a.end()); -} - - -template -void Foam::sort(List& a, const Cmp& cmp) -{ - std::sort(a.begin(), a.end(), cmp); -} - - -template -void Foam::stableSort(List& a) -{ - std::stable_sort(a.begin(), a.end()); -} - - -template -void Foam::stableSort(List& a, const Cmp& cmp) -{ - std::stable_sort(a.begin(), a.end(), cmp); -} - - // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // // Assignment to UList operator. Takes linear time. diff --git a/src/OpenFOAM/containers/Lists/List/List.H b/src/OpenFOAM/containers/Lists/List/List.H index a6bf29686a..b89f6c95a9 100644 --- a/src/OpenFOAM/containers/Lists/List/List.H +++ b/src/OpenFOAM/containers/Lists/List/List.H @@ -248,18 +248,6 @@ public: template List readList(Istream&); -template -void sort(List&); - -template -void sort(List&, const Cmp&); - -template -void stableSort(List&); - -template -void stableSort(List&, const Cmp&); - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/containers/Lists/UList/UList.C b/src/OpenFOAM/containers/Lists/UList/UList.C index 255a91ecfe..42d00bc41e 100644 --- a/src/OpenFOAM/containers/Lists/UList/UList.C +++ b/src/OpenFOAM/containers/Lists/UList/UList.C @@ -30,6 +30,8 @@ License #include "ListLoopM.H" #include "contiguous.H" +#include + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template @@ -116,6 +118,34 @@ Foam::label Foam::UList::byteSize() const } +template +void Foam::sort(UList& a) +{ + std::sort(a.begin(), a.end()); +} + + +template +void Foam::sort(UList& a, const Cmp& cmp) +{ + std::sort(a.begin(), a.end(), cmp); +} + + +template +void Foam::stableSort(UList& a) +{ + std::stable_sort(a.begin(), a.end()); +} + + +template +void Foam::stableSort(UList& a, const Cmp& cmp) +{ + std::stable_sort(a.begin(), a.end(), cmp); +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/containers/Lists/UList/UList.H b/src/OpenFOAM/containers/Lists/UList/UList.H index f7db8226ff..087f685ff6 100644 --- a/src/OpenFOAM/containers/Lists/UList/UList.H +++ b/src/OpenFOAM/containers/Lists/UList/UList.H @@ -320,6 +320,18 @@ public: ); }; +template +void sort(UList&); + +template +void sort(UList&, const Cmp&); + +template +void stableSort(UList&); + +template +void stableSort(UList&, const Cmp&); + // Reverse the first n elements of the list template inline void reverse(UList&, const label n); diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H index 0ed207057a..e394f679d0 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H @@ -106,6 +106,15 @@ public: //- Global sum of localSizes inline label size() const; + //- Size of procI data + inline label localSize(const label procI) const; + + //- From local to global on procI + inline label toGlobal(const label procI, const label i) const; + + //- Is on processor procI + inline bool isLocal(const label procI, const label i) const; + //- From global to local on procI inline label toLocal(const label procI, const label i) const; @@ -115,9 +124,6 @@ public: //- Start of procI data inline label offset(const label procI) const; - //- Size of procI data - inline label localSize(const label procI) const; - // IOstream Operators diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H index f9444fe54b..d7c44107c6 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H @@ -63,27 +63,39 @@ inline Foam::label Foam::globalIndex::size() const } -inline Foam::label Foam::globalIndex::toGlobal(const label i) const +inline Foam::label Foam::globalIndex::toGlobal +( + const label procI, + const label i +) const { - return - ( - Pstream::myProcNo() == 0 - ? i - : i + offsets_[Pstream::myProcNo()-1] - ); + return(procI == 0 ? i : i + offsets_[procI-1]); } + +inline Foam::label Foam::globalIndex::toGlobal(const label i) const +{ + return toGlobal(Pstream::myProcNo(), i); +} + + //- Is on local processor -inline bool Foam::globalIndex::isLocal(const label i) const +inline bool Foam::globalIndex::isLocal(const label procI, const label i) const { return - (i < offsets_[Pstream::myProcNo()]) - && (i >= (Pstream::myProcNo() == 0 ? 0 : offsets_[Pstream::myProcNo()-1])); + (i < offsets_[procI]) + && (i >= (procI == 0 ? 0 : offsets_[procI-1])); +} + + +inline bool Foam::globalIndex::isLocal(const label i) const +{ + return isLocal(Pstream::myProcNo(), i); } inline Foam::label Foam::globalIndex::toLocal(const label procI, const label i) - const +const { label localI = (procI == 0 ? i : i - offsets_[procI-1]); diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H index 690fb1083b..de1d9cc9f9 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H @@ -32,7 +32,6 @@ Description SourceFiles coupledPolyPatch.C - coupledPolyPatchMorph.C \*---------------------------------------------------------------------------*/ diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 33b3cab15e..e43d4b9d29 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -607,6 +607,18 @@ inline bool triangle::classify // system E0, E1 // + //Pout<< "alpha:" << alpha << endl; + //Pout<< "beta:" << beta << endl; + //Pout<< "hit:" << hit << endl; + //Pout<< "tol:" << tol << endl; + + if (hit) + { + // alpha,beta might get negative due to precision errors + alpha = max(0.0, min(1.0, alpha)); + beta = max(0.0, min(1.0, beta)); + } + nearType = NONE; nearLabel = -1; diff --git a/src/fvMotionSolver/Make/files b/src/fvMotionSolver/Make/files index 41f31c0ecc..f229ae322e 100644 --- a/src/fvMotionSolver/Make/files +++ b/src/fvMotionSolver/Make/files @@ -30,5 +30,6 @@ pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPo pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C +pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C LIB = $(FOAM_LIBBIN)/libfvMotionSolvers diff --git a/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000..bbafb189fa --- /dev/null +++ b/src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C @@ -0,0 +1,514 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "surfaceDisplacementPointPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "transformField.H" +#include "fvMesh.H" +#include "displacementLaplacianFvMotionSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +template<> +const char* +NamedEnum:: +names[] = +{ + "nearest", + "pointNormal", + "fixedNormal" +}; + +const NamedEnum + surfaceDisplacementPointPatchVectorField::projectModeNames_; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void surfaceDisplacementPointPatchVectorField::calcProjection +( + vectorField& displacement +) const +{ + const polyMesh& mesh = patch().boundaryMesh().mesh()(); + const pointField& localPoints = patch().localPoints(); + const labelList& meshPoints = patch().meshPoints(); + + //const scalar deltaT = mesh.time().deltaT().value(); + + // Construct large enough vector in direction of projectDir so + // we're guaranteed to hit something. + + //- Per point projection vector: + const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min()); + + // For case of fixed projection vector: + vector projectVec; + if (projectMode_ == FIXEDNORMAL) + { + vector n = projectDir_/mag(projectDir_); + projectVec = projectLen*n; + } + + + // Get fixed points (bit of a hack) + const pointZone* zonePtr = NULL; + + if (frozenPointsZone_.size() > 0) + { + const pointZoneMesh& pZones = mesh.pointZones(); + + zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)]; + + Pout<< "surfaceDisplacementPointPatchVectorField : Fixing all " + << zonePtr->size() << " points in pointZone " << zonePtr->name() + << endl; + } + + // Get the starting locations from the motionSolver + const pointField& points0 = mesh.lookupObject + ( + "dynamicMeshDict" + ).points0(); + + + pointField start(meshPoints.size()); + forAll(start, i) + { + start[i] = points0[meshPoints[i]] + displacement[i]; + } + + label nNotProjected = 0; + + if (projectMode_ == NEAREST) + { + List nearest; + labelList hitSurfaces; + surfaces().findNearest + ( + start, + scalarField(start.size(), sqr(projectLen)), + hitSurfaces, + nearest + ); + + forAll(nearest, i) + { + if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0)) + { + // Fixed point. Reset to point0 location. + displacement[i] = points0[meshPoints[i]] - localPoints[i]; + } + else if (nearest[i].hit()) + { + displacement[i] = + nearest[i].hitPoint() + - points0[meshPoints[i]]; + } + else + { + nNotProjected++; + + if (debug) + { + Pout<< " point:" << meshPoints[i] + << " coord:" << localPoints[i] + << " did not find any surface within " << projectLen + << endl; + } + } + } + } + else + { + // Do tests on all points. Combine later on. + + // 1. Check if already on surface + List nearest; + { + labelList nearestSurface; + surfaces().findNearest + ( + start, + scalarField(start.size(), sqr(SMALL)), + nearestSurface, + nearest + ); + } + + // 2. intersection. (combined later on with information from nearest + // above) + vectorField projectVecs(start.size(), projectVec); + + if (projectMode_ == POINTNORMAL) + { + projectVecs = projectLen*patch().pointNormals(); + } + + // Knock out any wedge component + scalarField offset(start.size(), 0.0); + if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents) + { + forAll(offset, i) + { + offset[i] = start[i][wedgePlane_]; + start[i][wedgePlane_] = 0; + projectVecs[i][wedgePlane_] = 0; + } + } + + List rightHit; + { + labelList rightSurf; + surfaces().findAnyIntersection + ( + start, + start+projectVecs, + rightSurf, + rightHit + ); + } + + List leftHit; + { + labelList leftSurf; + surfaces().findAnyIntersection + ( + start, + start-projectVecs, + leftSurf, + leftHit + ); + } + + // 3. Choose either -fixed, nearest, right, left. + forAll(displacement, i) + { + if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0)) + { + // Fixed point. Reset to point0 location. + displacement[i] = points0[meshPoints[i]] - localPoints[i]; + } + else if (nearest[i].hit()) + { + // Found nearest. + displacement[i] = + nearest[i].hitPoint() + - points0[meshPoints[i]]; + } + else + { + pointIndexHit interPt; + + if (rightHit[i].hit()) + { + if (leftHit[i].hit()) + { + if + ( + magSqr(rightHit[i].hitPoint()-start[i]) + < magSqr(leftHit[i].hitPoint()-start[i]) + ) + { + interPt = rightHit[i]; + } + else + { + interPt = leftHit[i]; + } + } + else + { + interPt = rightHit[i]; + } + } + else + { + if (leftHit[i].hit()) + { + interPt = leftHit[i]; + } + } + + + if (interPt.hit()) + { + if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents) + { + interPt.rawPoint()[wedgePlane_] += offset[i]; + } + displacement[i] = interPt.rawPoint()-points0[meshPoints[i]]; + } + else + { + nNotProjected++; + + if (debug) + { + Pout<< " point:" << meshPoints[i] + << " coord:" << localPoints[i] + << " did not find any intersection between" + << " ray from " << start[i]-projectVecs[i] + << " to " << start[i]+projectVecs[i] << endl; + } + } + } + } + } + + reduce(nNotProjected, sumOp