From 6eb1f6fefc338fdc04609811fa7a14e1e4f90628 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 8 Oct 2018 17:11:37 +0100 Subject: [PATCH 1/2] ENH: searchableSurface: handle multiple surfaces. Fixes #1034 --- .../searchableSurfacesQueries.C | 61 +++++++++++-------- 1 file changed, 36 insertions(+), 25 deletions(-) diff --git a/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C b/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C index f5c48914a2..df5234dc1a 100644 --- a/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C +++ b/src/meshTools/searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C @@ -516,37 +516,48 @@ void Foam::searchableSurfacesQueries::findNearest // - current surface : info+normal1 forAll(near, i) { - if (info[i].hit() && normal[i] != vector::zero) + if (info[i].hit()) { - if (mag(normal[i]&normal1[i]) < 1.0-1e-6) + if (normal[i] != vector::zero) { - plane pl0(near[i], normal[i], false); - plane pl1(info[i].hitPoint(), normal1[i], false); - - plane::ray r(pl0.planeIntersect(pl1)); - vector n = r.dir() / mag(r.dir()); - - // Calculate vector to move onto intersection line - vector d(r.refPoint()-near[i]); - d -= (d&n)*n; - - // Trim the max distance - scalar magD = mag(d); - if (magD > SMALL) + // Have previous hit. Find intersection + if (mag(normal[i]&normal1[i]) < 1.0-1e-6) { - scalar maxDist = Foam::sqrt(distSqr[i]); - if (magD > maxDist) - { - // Clip - d /= magD; - d *= maxDist; - } + plane pl0(near[i], normal[i], false); + plane pl1(info[i].hitPoint(), normal1[i], false); - near[i] += d; - normal[i] = normal1[i]; - constraint[i].applyConstraint(normal1[i]); + plane::ray r(pl0.planeIntersect(pl1)); + vector n = r.dir() / mag(r.dir()); + + // Calculate vector to move onto intersection line + vector d(r.refPoint()-near[i]); + d -= (d&n)*n; + + // Trim the max distance + scalar magD = mag(d); + if (magD > SMALL) + { + scalar maxDist = Foam::sqrt(distSqr[i]); + if (magD > maxDist) + { + // Clip + d /= magD; + d *= maxDist; + } + + near[i] += d; + normal[i] = normal1[i]; + constraint[i].applyConstraint(normal1[i]); + } } } + else + { + // First hit + near[i] = info[i].hitPoint(); + normal[i] = normal1[i]; + constraint[i].applyConstraint(normal1[i]); + } } } From 15747fde6b96946db63ec9457daf05d5b30c67e0 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 10 Oct 2018 17:03:04 +0100 Subject: [PATCH 2/2] ENH: overset: on-demand trigger of stencil-update. Fixes #1028. --- .../cellCellStencil/cellCellStencil.H | 15 ++- .../cellCellStencilTemplates.C | 62 +++++++++ .../cellVolumeWeightCellCellStencil.C | 8 +- .../cellVolumeWeightCellCellStencil.H | 8 +- .../inverseDistanceCellCellStencil.C | 127 ++++++++++-------- .../inverseDistanceCellCellStencil.H | 10 +- .../trackingInverseDistanceCellCellStencil.C | 16 ++- .../dynamicOversetFvMesh.C | 10 -- .../dynamicOversetFvMesh.H | 2 - 9 files changed, 178 insertions(+), 80 deletions(-) create mode 100644 src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H index 250afd7ba9..9bb1760188 100644 --- a/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H +++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -43,8 +43,8 @@ SourceFiles #include "scalarList.H" #include "mapDistribute.H" -#include "fvMesh.H" #include "pointList.H" +#include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -95,6 +95,11 @@ protected: //- Count ocurrences (in parallel) static labelList count(const label size, const labelUList& lst); + //- Helper: create volScalarField for postprocessing. + template + tmp createField(const word& name, UList&) const; + + private: // Private Member Functions @@ -220,6 +225,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "cellCellStencilTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C new file mode 100644 index 0000000000..d431890976 --- /dev/null +++ b/src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C @@ -0,0 +1,62 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 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 "zeroGradientFvPatchFields.H" + +// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * // + +template +Foam::tmp +Foam::cellCellStencil::createField(const word& name, UList& psi) const +{ + tmp tfld + ( + new volScalarField + ( + IOobject + ( + name, + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh_, + dimensionedScalar(dimless, Zero), + zeroGradientFvPatchScalarField::typeName + ) + ); + volScalarField& fld = tfld.ref(); + + forAll(psi, cellI) + { + fld[cellI] = psi[cellI]; + } + return tfld; +} + + +// ************************************************************************* // diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C index 7168900d6c..eca71ca548 100644 --- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C +++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C @@ -1113,8 +1113,10 @@ bool Foam::cellCellStencils::cellVolumeWeight::update() List> compactMap; - mapDistribute map(globalCells, cellStencil_, compactMap); - cellInterpolationMap_.transfer(map); + cellInterpolationMap_.reset + ( + new mapDistribute(globalCells, cellStencil_, compactMap) + ); // Dump interpolation stencil if (debug) @@ -1128,7 +1130,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update() OBJstream str(mesh_.time().timePath()/"stencil2.obj"); Info<< typeName << " : dumping to " << str.name() << endl; pointField cc(mesh_.cellCentres()); - cellInterpolationMap_.distribute(cc); + cellInterpolationMap().distribute(cc); forAll(interpolationCells_, compactI) { diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H index 321922a8a9..6de2e2fb15 100644 --- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H +++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2016-2017 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2016-2018 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -77,7 +77,7 @@ protected: labelList interpolationCells_; //- Fetch interpolated cells - mapDistribute cellInterpolationMap_; + autoPtr cellInterpolationMap_; //- Interpolation stencil labelListList cellStencil_; @@ -205,6 +205,10 @@ public: //- Return a communication schedule virtual const mapDistribute& cellInterpolationMap() const { + if (!cellInterpolationMap_.valid()) + { + const_cast(*this).update(); + } return cellInterpolationMap_; } diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C index 9e26175a19..e626bacede 100644 --- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C +++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C @@ -1500,7 +1500,7 @@ void Foam::cellCellStencils::inverseDistance::createStencil while (true) { - pointField samples(cellInterpolationMap_.constructSize(), greatPoint); + pointField samples(cellInterpolationMap().constructSize(), greatPoint); // Fill remote slots (override old content). We'll find out later // on which one has won and mark this one in doneAcceptor. @@ -1544,9 +1544,9 @@ void Foam::cellCellStencils::inverseDistance::createStencil Pstream::commsTypes::nonBlocking, List(), mesh_.nCells(), - cellInterpolationMap_.constructMap(), + cellInterpolationMap().constructMap(), false, - cellInterpolationMap_.subMap(), + cellInterpolationMap().subMap(), false, samples, minMagSqrEqOp(), @@ -1597,9 +1597,9 @@ void Foam::cellCellStencils::inverseDistance::createStencil // Transfer the information back to the acceptor: // - donorCellCells : stencil (with first element the original donor) // - donorWeights : weights for donorCellCells - cellInterpolationMap_.distribute(donorCellCells); - cellInterpolationMap_.distribute(donorWeights); - cellInterpolationMap_.distribute(samples); + cellInterpolationMap().distribute(donorCellCells); + cellInterpolationMap().distribute(donorWeights); + cellInterpolationMap().distribute(samples); // Check which acceptor has won and transfer forAll(interpolationCells_, i) @@ -1635,8 +1635,15 @@ void Foam::cellCellStencils::inverseDistance::createStencil // Re-do the mapDistribute List> compactMap; - mapDistribute map(globalCells, cellStencil_, compactMap); - cellInterpolationMap_.transfer(map); + cellInterpolationMap_.reset + ( + new mapDistribute + ( + globalCells, + cellStencil_, + compactMap + ) + ); } @@ -1982,6 +1989,12 @@ bool Foam::cellCellStencils::inverseDistance::update() } } + if (debug) + { + tmp tfld(createField("allCellTypes", allCellTypes)); + //tfld.ref().correctBoundaryConditions(); + tfld().write(); + } // Use the patch types and weights to decide what to do forAll(allPatchTypes, cellI) @@ -2017,15 +2030,55 @@ bool Foam::cellCellStencils::inverseDistance::update() } } + if (debug) + { + tmp tfld + ( + createField("allCellTypes_patch", allCellTypes) + ); + //tfld.ref().correctBoundaryConditions(); + tfld().write(); + } // Mark unreachable bits findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes); + if (debug) + { + tmp tfld + ( + createField("allCellTypes_hole", allCellTypes) + ); + //tfld.ref().correctBoundaryConditions(); + tfld().write(); + } + if (debug) + { + labelList stencilSize(mesh_.nCells()); + forAll(allStencil, celli) + { + stencilSize[celli] = allStencil[celli].size(); + } + tmp tfld(createField("allStencil_hole", stencilSize)); + //tfld.ref().correctBoundaryConditions(); + tfld().write(); + } + // Add buffer interpolation layer(s) around holes scalarField allWeight(mesh_.nCells(), 0.0); walkFront(layerRelax, allStencil, allCellTypes, allWeight); + if (debug) + { + tmp tfld + ( + createField("allCellTypes_front", allCellTypes) + ); + //tfld.ref().correctBoundaryConditions(); + tfld().write(); + } + // Check previous iteration cellTypes_ for any hole->calculated changes { @@ -2092,9 +2145,10 @@ bool Foam::cellCellStencils::inverseDistance::update() interpolationCells_.transfer(interpolationCells); List> compactMap; - mapDistribute map(globalCells, cellStencil_, compactMap); - cellInterpolationMap_.transfer(map); - + cellInterpolationMap_.reset + ( + new mapDistribute(globalCells, cellStencil_, compactMap) + ); cellInterpolationWeight_.transfer(allWeight); cellInterpolationWeight_.correctBoundaryConditions(); @@ -2110,7 +2164,7 @@ bool Foam::cellCellStencils::inverseDistance::update() Pout<< type() << " : dumping injectionStencil to " << str.name() << endl; pointField cc(mesh_.cellCentres()); - cellInterpolationMap_.distribute(cc); + cellInterpolationMap().distribute(cc); forAll(cellStencil_, celli) { @@ -2140,21 +2194,7 @@ bool Foam::cellCellStencils::inverseDistance::update() // Dump max weight { - volScalarField maxMagWeight - ( - IOobject - ( - "maxMagWeight", - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar(dimless, Zero), - zeroGradientFvPatchScalarField::typeName - ); + scalarField maxMagWeight(mesh_.nCells(), 0.0); forAll(cellStencil_, celli) { const scalarList& wghts = cellInterpolationWeights_[celli]; @@ -2180,44 +2220,25 @@ bool Foam::cellCellStencils::inverseDistance::update() << endl; } } - maxMagWeight.correctBoundaryConditions(); - maxMagWeight.write(); + tmp tfld(createField("maxMagWeight", maxMagWeight)); + tfld.ref().correctBoundaryConditions(); + tfld().write(); } // Dump cell types { - volScalarField volTypes - ( - IOobject - ( - "cellTypes", - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar(dimless, Zero), - zeroGradientFvPatchScalarField::typeName - ); - - forAll(volTypes.internalField(), cellI) - { - volTypes[cellI] = cellTypes_[cellI]; - } - volTypes.correctBoundaryConditions(); - volTypes.write(); + tmp tfld(createField("cellTypes", cellTypes_)); + tfld.ref().correctBoundaryConditions(); + tfld().write(); } - // Dump stencil mkDir(mesh_.time().timePath()); OBJstream str(mesh_.time().timePath()/"stencil.obj"); Pout<< type() << " : dumping to " << str.name() << endl; pointField cc(mesh_.cellCentres()); - cellInterpolationMap_.distribute(cc); + cellInterpolationMap().distribute(cc); forAll(cellStencil_, celli) { diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H index 22bf2531c8..ae7da72beb 100644 --- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H +++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -84,7 +84,7 @@ protected: labelList interpolationCells_; //- Fetch interpolated cells - mapDistribute cellInterpolationMap_; + autoPtr cellInterpolationMap_; //- Interpolation stencil labelListList cellStencil_; @@ -302,7 +302,11 @@ public: //- Return a communication schedule virtual const mapDistribute& cellInterpolationMap() const { - return cellInterpolationMap_; + if (!cellInterpolationMap_.valid()) + { + const_cast(*this).update(); + } + return cellInterpolationMap_(); } //- Per interpolated cell the neighbour cells (in terms of slots as diff --git a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C index e595d07c30..98d5df0525 100644 --- a/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C +++ b/src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C @@ -787,9 +787,15 @@ bool Foam::cellCellStencils::trackingInverseDistance::update() interpolationCells_.transfer(interpolationCells); List> compactMap; - mapDistribute map(globalCells, cellStencil_, compactMap); - cellInterpolationMap_.transfer(map); - + cellInterpolationMap_.reset + ( + new mapDistribute + ( + globalCells, + cellStencil_, + compactMap + ) + ); cellInterpolationWeight_.transfer(allWeight); cellInterpolationWeight_.correctBoundaryConditions(); @@ -802,7 +808,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update() Pout<< typeName << " : dumping injectionStencil to " << str.name() << endl; pointField cc(mesh_.cellCentres()); - cellInterpolationMap_.distribute(cc); + cellInterpolationMap().distribute(cc); forAll(cellStencil_, celli) { @@ -862,7 +868,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update() OBJstream str(mesh_.time().timePath()/"stencil.obj"); Pout<< typeName << " : dumping to " << str.name() << endl; pointField cc(mesh_.cellCentres()); - cellInterpolationMap_.distribute(cc); + cellInterpolationMap().distribute(cc); forAll(cellStencil_, celli) { diff --git a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C index 05969e3245..378516e334 100644 --- a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C +++ b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C @@ -44,15 +44,6 @@ namespace Foam bool Foam::dynamicOversetFvMesh::updateAddressing() const { const cellCellStencilObject& overlap = Stencil::New(*this); - // Make sure that stencil is not still in the initial, non-uptodate - // state. This gets triggered by e.g. Poisson wall distance that - // triggers the stencil even though time has not been updated (and hence - // mesh.update() has not been called. - if (!stencilIsUpToDate_) - { - stencilIsUpToDate_ = true; - const_cast(overlap).update(); - } // The (processor-local part of the) stencil determines the local // faces to add to the matrix. tbd: parallel @@ -236,7 +227,6 @@ Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io) { // Load stencil (but do not update) (void)Stencil::New(*this, false); - stencilIsUpToDate_ = false; } diff --git a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H index 70a064516c..21d875dea1 100644 --- a/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H +++ b/src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H @@ -65,8 +65,6 @@ protected: // lduAddressing (true) mutable bool active_; - mutable bool stencilIsUpToDate_; - //- Extended addressing (extended with local interpolation stencils) mutable autoPtr lduPtr_;