zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: mappedPatchBase: make robust w.r.t. warped faces

Browse Source
This commit is contained in:
mattijs
2011-12-13 12:43:55 +00:00
parent 2d7b3363ed
commit 2499ad10a7
9 changed files with 537 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/OpenFOAM/meshes/polyMesh/polyMesh.C
View File

@ -1172,6 +1172,8 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
)
).movePoints(points_);
}
meshSearchMeshObject::Delete(*this);
meshSearchFACECENTRETETSMeshObject::Delete(*this);
return sweptVols;
}

2
src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C
View File

@ -95,6 +95,8 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
)
).updateMesh(mpm);
}
meshSearchMeshObject::Delete(*this);
meshSearchFACECENTRETETSMeshObject::Delete(*this);
}

2
src/meshTools/Make/files
View File

@ -27,6 +27,8 @@ polyMeshZipUpCells/polyMeshZipUpCells.C
primitiveMeshGeometry/primitiveMeshGeometry.C
meshSearch/meshSearch.C
meshSearch/meshSearchFACECENTRETETSMeshObject.C
meshSearch/meshSearchMeshObject.C
meshTools/meshTools.C

317
src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.C
View File

@ -26,7 +26,7 @@ License
#include "mappedPatchBase.H"
#include "addToRunTimeSelectionTable.H"
#include "ListListOps.H"
#include "meshSearch.H"
#include "meshSearchMeshObject.H"
#include "meshTools.H"
#include "OFstream.H"
#include "Random.H"
@ -37,6 +37,7 @@ License
#include "Time.H"
#include "mapDistribute.H"
#include "SubField.H"
#include "triPointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,64 +81,144 @@ const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::facePoints
(
const polyPatch& pp
) const
{
const polyMesh& mesh = pp.boundaryMesh().mesh();
const labelList& tetBasePts = mesh.tetBasePtIs();
const pointField& p = pp.points();
const pointField& cellCentres = mesh.cellCentres();
const labelList& faceCells = pp.faceCells();
// Initialise to face-centre
tmp<pointField> tfacePoints(new pointField(patch_.faceCentres()));
pointField& facePoints = tfacePoints();
forAll(pp, faceI)
{
const face& f = pp[faceI];
if (f.size() <= 3)
{
// Point already on tets of cell.
}
else
{
// Find intersection of (face-centre-decomposition) centre to
// cell-centre with face-diagonal-decomposition triangles.
const point& cc = cellCentres[faceCells[faceI]];
vector d = facePoints[faceI]-cc;
const label fp0 = tetBasePts[faceI];
const point& basePoint = p[f[fp0]];
//bool foundPoint = false;
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); i++)
{
const point& thisPoint = p[f[fp]];
label nextFp = f.fcIndex(fp);
const point& nextPoint = p[f[nextFp]];
const triPointRef tri(basePoint, thisPoint, nextPoint);
pointHit hitInfo = tri.intersection
(
cc,
d,
intersection::HALF_RAY
);
if (hitInfo.hit() && hitInfo.distance() > 0)
{
facePoints[faceI] = hitInfo.hitPoint();
//foundPoint = true;
break;
}
fp = nextFp;
}
//if (!foundPoint)
//{
// Pout<< "cell:" << faceCells[faceI] << " at:" << cc
// << " face-centre:" << patch_.faceCentres()[faceI]
// << " did not find any intersection." << endl;
//}
}
}
return tfacePoints;
}
void Foam::mappedPatchBase::collectSamples
(
const pointField& facePoints,
pointField& samples,
labelList& patchFaceProcs,
labelList& patchFaces,
pointField& patchFc
) const
{
// Collect all sample points and the faces they come from.
List<pointField> globalFc(Pstream::nProcs());
List<pointField> globalSamples(Pstream::nProcs());
labelListList globalFaces(Pstream::nProcs());
{
List<pointField> globalFc(Pstream::nProcs());
globalFc[Pstream::myProcNo()] = facePoints;
Pstream::gatherList(globalFc);
Pstream::scatterList(globalFc);
// Rework into straight list
patchFc = ListListOps::combine<pointField>
(
globalFc,
accessOp<pointField>()
);
}
globalFc[Pstream::myProcNo()] = patch_.faceCentres();
globalSamples[Pstream::myProcNo()] = samplePoints();
globalFaces[Pstream::myProcNo()] = identity(patch_.size());
{
List<pointField> globalSamples(Pstream::nProcs());
globalSamples[Pstream::myProcNo()] = samplePoints(facePoints);
Pstream::gatherList(globalSamples);
Pstream::scatterList(globalSamples);
// Rework into straight list
samples = ListListOps::combine<pointField>
(
globalSamples,
accessOp<pointField>()
);
}
// Distribute to all processors
Pstream::gatherList(globalSamples);
Pstream::scatterList(globalSamples);
Pstream::gatherList(globalFaces);
Pstream::scatterList(globalFaces);
Pstream::gatherList(globalFc);
Pstream::scatterList(globalFc);
{
labelListList globalFaces(Pstream::nProcs());
globalFaces[Pstream::myProcNo()] = identity(patch_.size());
// Distribute to all processors
Pstream::gatherList(globalFaces);
Pstream::scatterList(globalFaces);
// Rework into straight list
samples = ListListOps::combine<pointField>
(
globalSamples,
accessOp<pointField>()
);
patchFaces = ListListOps::combine<labelList>
(
globalFaces,
accessOp<labelList>()
);
patchFc = ListListOps::combine<pointField>
(
globalFc,
accessOp<pointField>()
);
patchFaceProcs.setSize(patchFaces.size());
labelList nPerProc
(
ListListOps::subSizes
patchFaces = ListListOps::combine<labelList>
(
globalFaces,
accessOp<labelList>()
)
);
label sampleI = 0;
forAll(nPerProc, procI)
);
}
{
for (label i = 0; i < nPerProc[procI]; i++)
labelList nPerProc(Pstream::nProcs());
nPerProc[Pstream::myProcNo()] = patch_.size();
Pstream::gatherList(nPerProc);
Pstream::scatterList(nPerProc);
patchFaceProcs.setSize(patchFaces.size());
label sampleI = 0;
forAll(nPerProc, procI)
{
patchFaceProcs[sampleI++] = procI;
for (label i = 0; i < nPerProc[procI]; i++)
{
patchFaceProcs[sampleI++] = procI;
}
}
}
}
@ -173,8 +254,9 @@ void Foam::mappedPatchBase::findSamples
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
// Octree based search engine
meshSearch meshSearchEngine(mesh);
//- Note: face-diagonal decomposition
const meshSearchMeshObject& meshSearchEngine =
meshSearchMeshObject::New(mesh);
forAll(samples, sampleI)
{
@ -290,8 +372,9 @@ void Foam::mappedPatchBase::findSamples
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
// Octree based search engine
meshSearch meshSearchEngine(mesh);
//- Note: face-diagonal decomposition
const meshSearchMeshObject& meshSearchEngine =
meshSearchMeshObject::New(mesh);
forAll(samples, sampleI)
{
@ -355,23 +438,6 @@ void Foam::mappedPatchBase::findSamples
}
}
// Check for samples not being found
forAll(nearest, sampleI)
{
if (!nearest[sampleI].first().hit())
{
FatalErrorIn
(
"mappedPatchBase::findSamples"
"(const pointField&, labelList&"
", labelList&, pointField&)"
) << "Did not find sample " << samples[sampleI]
<< " on any processor of region " << sampleRegion_
<< exit(FatalError);
}
}
// Convert back into proc+local index
sampleProcs.setSize(samples.size());
sampleIndices.setSize(samples.size());
@ -379,21 +445,40 @@ void Foam::mappedPatchBase::findSamples
forAll(nearest, sampleI)
{
sampleProcs[sampleI] = nearest[sampleI].second().second();
sampleIndices[sampleI] = nearest[sampleI].first().index();
sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
if (!nearest[sampleI].first().hit())
{
sampleProcs[sampleI] = -1;
sampleIndices[sampleI] = -1;
sampleLocations[sampleI] = vector::max;
}
else
{
sampleProcs[sampleI] = nearest[sampleI].second().second();
sampleIndices[sampleI] = nearest[sampleI].first().index();
sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
}
}
}
void Foam::mappedPatchBase::calcMapping() const
{
static bool hasWarned = false;
if (mapPtr_.valid())
{
FatalErrorIn("mappedPatchBase::calcMapping() const")
<< "Mapping already calculated" << exit(FatalError);
}
// Get points on face (since cannot use face-centres - might be off
// face-diagonal decomposed tets.
tmp<pointField> patchPoints(facePoints(patch_));
// Get offsetted points
const pointField offsettedPoints = samplePoints(patchPoints());
// Do a sanity check
// Am I sampling my own patch? This only makes sense for a non-zero
// offset.
@ -405,8 +490,10 @@ void Foam::mappedPatchBase::calcMapping() const
);
// Check offset
vectorField d(samplePoints()-patch_.faceCentres());
if (sampleMyself && gAverage(mag(d)) <= ROOTVSMALL)
vectorField d(offsettedPoints-patchPoints());
bool coincident = (gAverage(mag(d)) <= ROOTVSMALL);
if (sampleMyself && coincident)
{
WarningIn
(
@ -438,7 +525,15 @@ void Foam::mappedPatchBase::calcMapping() const
labelList patchFaceProcs;
labelList patchFaces;
pointField patchFc;
collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
collectSamples
(
patchPoints,
samples,
patchFaceProcs,
patchFaces,
patchFc
);
// Find processor and cell/face samples are in and actual location.
labelList sampleProcs;
@ -446,6 +541,75 @@ void Foam::mappedPatchBase::calcMapping() const
pointField sampleLocations;
findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
// Check for samples that were not found. This will only happen for
// NEARESTCELL since finds cell containing a location
if (mode_ == NEARESTCELL)
{
label nNotFound = 0;
forAll(sampleProcs, sampleI)
{
if (sampleProcs[sampleI] == -1)
{
nNotFound++;
}
}
reduce(nNotFound, sumOp<label>());
if (nNotFound > 0)
{
if (!hasWarned)
{
WarningIn
(
"mappedPatchBase::mappedPatchBase\n"
"(\n"
" const polyPatch& pp,\n"
" const word& sampleRegion,\n"
" const sampleMode mode,\n"
" const word& samplePatch,\n"
" const vector& offset\n"
")\n"
) << "Did not find " << nNotFound
<< " out of " << sampleProcs.size() << " total samples."
<< " Sampling these on owner cell centre instead." << endl
<< "patch_:" << patch_.name() << endl
<< "sampleRegion_:" << sampleRegion_ << endl
<< "mode_:" << sampleModeNames_[mode_] << endl
<< "samplePatch_:" << samplePatch_ << endl
<< "offsetMode_:" << offsetModeNames_[offsetMode_] << endl;
hasWarned = true;
}
// Reset the samples that cannot be found to the cell centres.
pointField patchCc;
{
List<pointField> globalCc(Pstream::nProcs());
globalCc[Pstream::myProcNo()] = patch_.faceCellCentres();
Pstream::gatherList(globalCc);
Pstream::scatterList(globalCc);
patchCc = ListListOps::combine<pointField>
(
globalCc,
accessOp<pointField>()
);
}
forAll(sampleProcs, sampleI)
{
if (sampleProcs[sampleI] == -1)
{
// Reset to cell centres
samples[sampleI] = patchCc[sampleI];
}
}
// And re-search. Note: could be optimised to only search missing
// points.
findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
}
}
// Now we have all the data we need:
// - where sample originates from (so destination when mapping):
@ -613,7 +777,7 @@ void Foam::mappedPatchBase::calcAMI() const
/*
const polyPatch& nbr = samplePolyPatch();
// pointField nbrPoints(samplePoints());
// pointField nbrPoints(offsettedPoints());
pointField nbrPoints(nbr.localPoints());
if (debug)
@ -932,9 +1096,12 @@ const Foam::polyPatch& Foam::mappedPatchBase::samplePolyPatch() const
}
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints
(
const pointField& fc
) const
{
tmp<pointField> tfld(new pointField(patch_.faceCentres()));
tmp<pointField> tfld(new pointField(fc));
pointField& fld = tfld();
switch (offsetMode_)
@ -964,6 +1131,12 @@ Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
}
Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
{
return samplePoints(facePoints(patch_));
}
void Foam::mappedPatchBase::write(Ostream& os) const
{
os.writeKeyword("sampleMode") << sampleModeNames_[mode_]

10
src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBase.H
View File

@ -33,7 +33,7 @@ Description
sampleRegion region0;
// What to sample:
// - nearestCell : sample nearest cell
// - nearestCell : sample cell containing point
// - nearestPatchFace : nearest face on selected patch
// - nearestPatchFaceAMI : nearest face on selected patch
- patches need not conform
@ -204,9 +204,14 @@ protected:
// Protected Member Functions
//- Get the points from face-centre-decomposition face centres
// and project them onto the face-diagonal-decomposition triangles.
tmp<pointField> facePoints(const polyPatch&) const;
//- Collect single list of samples and originating processor+face.
void collectSamples
(
const pointField& facePoints,
pointField&,
labelList& patchFaceProcs,
labelList& patchFaces,
@ -222,6 +227,9 @@ protected:
pointField& sampleLocations // actual representative location
) const;
//- Get the sample points given the face points
tmp<pointField> samplePoints(const pointField&) const;
//- Calculate mapping
void calcMapping() const;

48
src/meshTools/meshSearch/meshSearchFACECENTRETETSMeshObject.C Normal file
View File

@ -0,0 +1,48 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "meshSearchFACECENTRETETSMeshObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(meshSearchFACECENTRETETSMeshObject, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::meshSearchFACECENTRETETSMeshObject::meshSearchFACECENTRETETSMeshObject
(
const polyMesh& mesh
)
:
MeshObject<polyMesh, meshSearchFACECENTRETETSMeshObject>(mesh),
meshSearch(mesh, polyMesh::FACECENTRETETS)
{}
// ************************************************************************* //

92
src/meshTools/meshSearch/meshSearchFACECENTRETETSMeshObject.H Normal file
View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
Foam::meshSearchFACECENTRETETSMeshObject
Description
MeshObject wrapper around meshSearch(mesh, polyMesh::FACECENTRETETS).
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef meshSearchFACECENTRETETSMeshObject_H
#define meshSearchFACECENTRETETSMeshObject_H
#include "MeshObject.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class meshSearchFACECENTRETETSMeshObject Declaration
\*---------------------------------------------------------------------------*/
class meshSearchFACECENTRETETSMeshObject
:
public MeshObject<polyMesh, meshSearchFACECENTRETETSMeshObject>,
public meshSearch
{
public:
// Declare name of the class and its debug switch
TypeName("meshSearchFACECENTRETETSMeshObject");
// Constructors
//- Constructor given polyMesh
explicit meshSearchFACECENTRETETSMeshObject(const polyMesh& mesh);
//- Destructor
virtual ~meshSearchFACECENTRETETSMeshObject()
{}
//
// // Member functions
//
// // Edit
//
// //- Update mesh topology using the morph engine
// void updateMesh();
//
// //- Correct weighting factors for moving mesh.
// bool movePoints();
//
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

45
src/meshTools/meshSearch/meshSearchMeshObject.C Normal file
View File

@ -0,0 +1,45 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "meshSearchMeshObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(meshSearchMeshObject, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::meshSearchMeshObject::meshSearchMeshObject(const polyMesh& mesh)
:
MeshObject<polyMesh, meshSearchMeshObject>(mesh),
meshSearch(mesh)
{}
// ************************************************************************* //

92
src/meshTools/meshSearch/meshSearchMeshObject.H Normal file
View File

@ -0,0 +1,92 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
Foam::meshSearchMeshObject
Description
MeshObject wrapper around meshSearch(mesh).
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef meshSearchMeshObject_H
#define meshSearchMeshObject_H
#include "MeshObject.H"
#include "meshSearch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class meshSearchMeshObject Declaration
\*---------------------------------------------------------------------------*/
class meshSearchMeshObject
:
public MeshObject<polyMesh, meshSearchMeshObject>,
public meshSearch
{
public:
// Declare name of the class and its debug switch
TypeName("meshSearchMeshObject");
// Constructors
//- Constructor given polyMesh
explicit meshSearchMeshObject(const polyMesh& mesh);
//- Destructor
virtual ~meshSearchMeshObject()
{}
//
// // Member functions
//
// // Edit
//
// //- Update mesh topology using the morph engine
// void updateMesh();
//
// //- Correct weighting factors for moving mesh.
// bool movePoints();
//
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //