songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: AMI - multiple updates

Browse Source 
- start of work to create a 1-to-1 face mapping across AMI patches
- faces are inserted according to the AMI addressing based on Horacio's method
- removed 'updated' flag and reworked some demand driven updates
- updated to handle 'walking' through baffles
- use bitSet instead of boolList
- moved update of meshPhi to movePoints() functions at fvPatch level
- moved scaling of areas to movePoints() functions at fvPatch level
- rehomed topology change code to own file
- added warning re: geometry construction

ACMI
- split srcMask into srcMask and srcAreaMask
  - former in range 0-1, and latter has bounding or tol to (1-tol) to avoid
    sigFpe's
This commit is contained in:
Andrew Heather 2019-03-07 20:25:23 +00:00 committed by Andrew Heather
parent b61dd6fd51
commit a13e00b5c4
43 changed files with 2435 additions and 654 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -45,6 +45,7 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
<< "Updating addressing and (optional) pointMesh/pointFields" << endl;
// Update boundaryMesh (note that patches themselves already ok)
// boundary_.updateMesh(mpm);
boundary_.updateMesh();
// Update zones

1
src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
View File

@ -63,6 +63,7 @@ void Foam::polyPatch::movePoints(PstreamBuffers&, const pointField& p)
primitivePatch::movePoints(p);
}
void Foam::polyPatch::updateMesh(PstreamBuffers&)
{
primitivePatch::clearGeom();

14
src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
View File

@ -55,6 +55,7 @@ namespace Foam
// Forward declarations
class polyBoundaryMesh;
class polyPatch;
class polyTopoChange;
class PstreamBuffers;
Ostream& operator<<(Ostream&, const polyPatch&);
@ -419,6 +420,19 @@ public:
labelList& rotation
) const;
//- For dynamic mesh cases - return true if this patch will change the
//- topology
virtual bool changeTopology() const
{
return false;
}
//- Collect topology changes in a polyTopoChange object
virtual bool setTopology(polyTopoChange&)
{
return false;
}
// Member operators

5
src/dynamicFvMesh/Make/files
View File

@ -10,4 +10,9 @@ dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C
simplifiedDynamicFvMesh/simplifiedDynamicFvMesh.C
dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C
LIB = $(FOAM_LIBBIN)/libdynamicFvMesh

3
src/dynamicFvMesh/dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
View File

@ -92,6 +92,9 @@ const Foam::motionSolver& Foam::dynamicMotionSolverFvMesh::motion() const
bool Foam::dynamicMotionSolverFvMesh::update()
{
// Scan through AMI patches and update
fvMesh::movePoints(motionPtr_->newPoints());
volVectorField* Uptr = getObjectPtr<volVectorField>("U");

216
src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.C Normal file
View File

@ -0,0 +1,216 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "dynamicMotionSolverFvMeshAMI.H"
#include "addToRunTimeSelectionTable.H"
#include "motionSolver.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "cyclicAMIPolyPatch.H"
#include "polyTopoChange.H"
#include "MeshObject.H"
#include "lduMesh.H"
#include "processorFvPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicMotionSolverFvMeshAMI, 0);
addToRunTimeSelectionTable
(
dynamicFvMesh,
dynamicMotionSolverFvMeshAMI,
IOobject
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicMotionSolverFvMeshAMI::dynamicMotionSolverFvMeshAMI
(
const IOobject& io
)
:
dynamicFvMesh(io),
motionPtr_(motionSolver::New(*this))
{}
Foam::dynamicMotionSolverFvMeshAMI::dynamicMotionSolverFvMeshAMI
(
const IOobject& io,
pointField&& points,
faceList&& faces,
labelList&& allOwner,
labelList&& allNeighbour,
const bool syncPar
)
:
dynamicFvMesh
(
io,
std::move(points),
std::move(faces),
std::move(allOwner),
std::move(allNeighbour),
syncPar
),
motionPtr_(motionSolver::New(*this))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::motionSolver& Foam::dynamicMotionSolverFvMeshAMI::motion() const
{
return *motionPtr_;
}
bool Foam::dynamicMotionSolverFvMeshAMI::update()
{
// Mesh not moved/changed yet
moving(false);
topoChanging(false);
if (debug)
{
for (const fvPatch& fvp : boundary())
{
if (!isA<processorFvPatch>(fvp))
{
Info<< "1 --- patch:" << fvp.patch().name()
<< " area:" << gSum(fvp.magSf()) << endl;
}
}
}
pointField newPoints(motionPtr_->curPoints());
polyBoundaryMesh& pbm = const_cast<polyBoundaryMesh&>(boundaryMesh());
// Scan all patches and see if we want to apply a mesh topology update
bool changeRequired = false;
for (polyPatch& pp : pbm)
{
DebugInfo
<< "pre-topology change: patch " << pp.name()
<< " size:" << returnReduce(pp.size(), sumOp<label>())
<< " mag(faceAreas):" << gSum(mag(pp.faceAreas())) << endl;
//changeRequired = pp.changeTopology(newPoints) || changeRequired;
changeRequired = pp.changeTopology() || changeRequired;
}
reduce(changeRequired, orOp<bool>());
if (changeRequired)
{
polyTopoChange polyTopo(*this);
// Set new point positions in polyTopo object
polyTopo.movePoints(newPoints);
// Accumulate the patch-based mesh changes on the current mesh
// Note:
// - updates the AMIs using the new points
// - creates a topo change object that removes old added faces and
// adds the new faces
for (polyPatch& pp : pbm)
{
pp.setTopology(polyTopo);
}
// Update geometry
// Note
// - changeMesh leads to polyMesh::resetPrimitives which will also
// trigger polyBoundaryMesh::updateMesh (init and update) and
// ::calcGeometry (with topoChanging = false)
// - BUT: mesh still corresponds to original (non-extended mesh) so
// we want to bypass these calls...
// - after changes topoChanging = true
autoPtr<mapPolyMesh> map =
polyTopo.changeMesh
(
*this,
true // We will be calling movePoints after this update
);
// Apply topology change - update fv geometry and map fields
// - polyMesh::updateMesh
// - fires initUpdateMesh and updateMesh in AMI BCs - called before
// mapFields
// - AMI addressing must be up-to-date - used by, e.g. FaceCellWave
// - will trigger (again) polyBoundaryMesh::updateMesh (init and update)
updateMesh(map());
// Move points and update derived properties
// Note:
// - resets face areas based on raw point locations!
// - polyBoundaryMesh::updateMesh (init and update)
// Note:
// - processorPolyPatches will trigger calculation of faceCentres
// (and therefore cell volumes), so need to update faceAreas in
// initMovePoints since proc patches will be evaluated later than
// AMI patches
if (map().hasMotionPoints())
{
movePoints(map().preMotionPoints());
}
}
else
{
fvMesh::movePoints(newPoints);
}
volVectorField* Uptr = getObjectPtr<volVectorField>("U");
if (Uptr)
{
Uptr->correctBoundaryConditions();
}
if (debug)
{
for (const fvPatch& fvp : boundary())
{
if (!isA<processorFvPatch>(fvp))
{
Info<< "2 --- patch:" << fvp.patch().name()
<< " area:" << gSum(fvp.magSf()) << endl;
}
}
}
return true;
}
// ************************************************************************* //

121
src/dynamicFvMesh/dynamicMotionSolverFvMeshAMI/dynamicMotionSolverFvMeshAMI.H Normal file
View File

@ -0,0 +1,121 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::dynamicMotionSolverFvMeshAMI
Description
The dynamicMotionSolverFvMeshAMI
SourceFiles
dynamicMotionSolverFvMeshAMI.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicMotionSolverFvMeshAMI_H
#define dynamicMotionSolverFvMeshAMI_H
#include "dynamicFvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class motionSolver;
/*---------------------------------------------------------------------------*\
Class dynamicMotionSolverFvMeshAMI Declaration
\*---------------------------------------------------------------------------*/
class dynamicMotionSolverFvMeshAMI
:
public dynamicFvMesh
{
// Private data
autoPtr<motionSolver> motionPtr_;
// Private Member Functions
//- No copy construct
dynamicMotionSolverFvMeshAMI
(
const dynamicMotionSolverFvMeshAMI&
) = delete;
//- No copy assignment
void operator=(const dynamicMotionSolverFvMeshAMI&) = delete;
public:
//- Runtime type information
TypeName("dynamicMotionSolverFvMeshAMI");
// Constructors
//- Construct from IOobject
dynamicMotionSolverFvMeshAMI(const IOobject& io);
//- Construct from components without boundary.
// Boundary is added using addFvPatches() member function
dynamicMotionSolverFvMeshAMI
(
const IOobject& io,
pointField&& points,
faceList&& faces,
labelList&& allOwner,
labelList&& allNeighbour,
const bool syncPar = true
);
//- Destructor
virtual ~dynamicMotionSolverFvMeshAMI() = default;
// Member Functions
//- Return the motionSolver
const motionSolver& motion() const;
//- Update the mesh for both mesh motion and topology change
virtual bool update();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

3
src/finiteVolume/Make/options
View File

@ -1,7 +1,8 @@
EXE_INC = \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
LIB_LIBS = \
-lOpenFOAM \

2
src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
View File

@ -173,7 +173,7 @@ public:
//- Return true if this patch field fixes a value
// Needed to check if a level has to be specified while solving
// Poissons equations
// Poisson equations
virtual bool fixesValue() const
{
const scalarField& mask =

21
src/finiteVolume/fvMesh/fvMesh.C
View File

@ -726,6 +726,8 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
{
DebugInFunction << endl;
// Grab old time volumes if the time has been incremented
// This will update V0, V00
if (curTimeIndex_ < time().timeIndex())
@ -733,6 +735,14 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
storeOldVol(V());
}
// Move the polyMesh and set the mesh motion fluxes to the swept-volumes
scalar rDeltaT = 1.0/time().deltaTValue();
tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
scalarField& sweptVols = tsweptVols.ref();
if (!phiPtr_)
{
// Create mesh motion flux
@ -761,14 +771,6 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
}
surfaceScalarField& phi = *phiPtr_;
// Move the polyMesh and set the mesh motion fluxes to the swept-volumes
scalar rDeltaT = 1.0/time().deltaTValue();
tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
scalarField& sweptVols = tsweptVols.ref();
phi.primitiveFieldRef() =
scalarField::subField(sweptVols, nInternalFaces());
phi.primitiveFieldRef() *= rDeltaT;
@ -776,7 +778,6 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
const fvPatchList& patches = boundary();
surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
forAll(patches, patchi)
{
phibf[patchi] = patches[patchi].patchSlice(sweptVols);
@ -805,6 +806,8 @@ Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
{
DebugInFunction << endl;
// Update polyMesh. This needs to keep volume existent!
polyMesh::updateMesh(mpm);

2
src/finiteVolume/fvMesh/fvMesh.H
View File

@ -91,6 +91,8 @@ class fvMesh
public fvSolution,
public data
{
protected:
// Private data
//- Boundary mesh

143
src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -26,9 +27,10 @@ License
\*---------------------------------------------------------------------------*/
#include "cyclicACMIFvPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "transform.H"
#include "surfaceFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -41,45 +43,13 @@ namespace Foam
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
void Foam::cyclicACMIFvPatch::updateAreas() const
void Foam::cyclicACMIFvPatch::resetPatchAreas(const fvPatch& fvp) const
{
if (cyclicACMIPolyPatch_.updated())
{
if (debug)
{
Pout<< "cyclicACMIFvPatch::updateAreas() : updating fv areas for "
<< name() << " and " << this->nonOverlapPatch().name()
<< endl;
}
const_cast<vectorField&>(fvp.Sf()) = fvp.patch().faceAreas();
const_cast<scalarField&>(fvp.magSf()) = mag(fvp.patch().faceAreas());
// owner couple
const_cast<vectorField&>(Sf()) = patch().faceAreas();
const_cast<scalarField&>(magSf()) = mag(patch().faceAreas());
// owner non-overlapping
const fvPatch& nonOverlapPatch = this->nonOverlapPatch();
const_cast<vectorField&>(nonOverlapPatch.Sf()) =
nonOverlapPatch.patch().faceAreas();
const_cast<scalarField&>(nonOverlapPatch.magSf()) =
mag(nonOverlapPatch.patch().faceAreas());
// neighbour couple
const cyclicACMIFvPatch& nbrACMI = neighbPatch();
const_cast<vectorField&>(nbrACMI.Sf()) =
nbrACMI.patch().faceAreas();
const_cast<scalarField&>(nbrACMI.magSf()) =
mag(nbrACMI.patch().faceAreas());
// neighbour non-overlapping
const fvPatch& nbrNonOverlapPatch = nbrACMI.nonOverlapPatch();
const_cast<vectorField&>(nbrNonOverlapPatch.Sf()) =
nbrNonOverlapPatch.patch().faceAreas();
const_cast<scalarField&>(nbrNonOverlapPatch.magSf()) =
mag(nbrNonOverlapPatch.patch().faceAreas());
// set the updated flag
cyclicACMIPolyPatch_.setUpdated(false);
}
DebugPout
<< fvp.patch().name() << " area:" << sum(fvp.magSf()) << endl;
}
@ -100,13 +70,13 @@ void Foam::cyclicACMIFvPatch::makeWeights(scalarField& w) const
)
);
scalar tol = cyclicACMIPolyPatch::tolerance();
const scalar tol = cyclicACMIPolyPatch::tolerance();
forAll(deltas, facei)
{
scalar di = deltas[facei];
scalar dni = nbrDeltas[facei];
scalar di = mag(deltas[facei]);
scalar dni = mag(nbrDeltas[facei]);
if (dni < tol)
{
@ -147,7 +117,6 @@ Foam::tmp<Foam::vectorField> Foam::cyclicACMIFvPatch::delta() const
vectorField nbrPatchD(interpolate(nbrPatch.coupledFvPatch::delta()));
tmp<vectorField> tpdv(new vectorField(patchD.size()));
vectorField& pdv = tpdv.ref();
@ -201,4 +170,94 @@ Foam::tmp<Foam::labelField> Foam::cyclicACMIFvPatch::internalFieldTransfer
}
void Foam::cyclicACMIFvPatch::movePoints()
{
if (!cyclicACMIPolyPatch_.owner())
{
return;
}
// Set the patch face areas to be consistent with the changes made at the
// polyPatch level
const fvPatch& nonOverlapPatch = this->nonOverlapPatch();
const cyclicACMIFvPatch& nbrACMI = neighbPatch();
const fvPatch& nbrNonOverlapPatch = nbrACMI.nonOverlapPatch();
resetPatchAreas(*this);
resetPatchAreas(nonOverlapPatch);
resetPatchAreas(nbrACMI);
resetPatchAreas(nbrNonOverlapPatch);
// Scale the mesh flux
const labelListList& newSrcAddr = AMI().srcAddress();
const labelListList& newTgtAddr = AMI().tgtAddress();
const fvMesh& mesh = boundaryMesh().mesh();
surfaceScalarField& meshPhi = const_cast<fvMesh&>(mesh).setPhi();
surfaceScalarField::Boundary& meshPhiBf = meshPhi.boundaryFieldRef();
// Note: phip and phiNonOverlapp will be different sizes if new faces
// have been added
scalarField& phip = meshPhiBf[cyclicACMIPolyPatch_.index()];
scalarField& phiNonOverlapp =
meshPhiBf[nonOverlapPatch.patch().index()];
forAll(phip, facei)
{
if (newSrcAddr[facei].empty())
{
// AMI patch with no connection to other coupled faces
phip[facei] = 0.0;
}
else
{
// Scale the mesh flux according to the area fraction
const face& fAMI = cyclicACMIPolyPatch_.localFaces()[facei];
// Note: using raw point locations to calculate the geometric
// area - faces areas are currently scaled (decoupled from
// mesh points)
const scalar geomArea = fAMI.mag(cyclicACMIPolyPatch_.localPoints());
phip[facei] *= magSf()[facei]/geomArea;
}
}
forAll(phiNonOverlapp, facei)
{
const scalar w = 1.0 - cyclicACMIPolyPatch_.srcMask()[facei];
phiNonOverlapp[facei] *= w;
}
scalarField& nbrPhip = meshPhiBf[nbrACMI.patch().index()];
scalarField& nbrPhiNonOverlapp =
meshPhiBf[nbrNonOverlapPatch.patch().index()];
forAll(nbrPhip, facei)
{
if (newTgtAddr[facei].empty())
{
nbrPhip[facei] = 0.0;
}
else
{
const face& fAMI = nbrACMI.patch().localFaces()[facei];
// Note: using raw point locations to calculate the geometric
// area - faces areas are currently scaled (decoupled from
// mesh points)
const scalar geomArea = fAMI.mag(nbrACMI.patch().localPoints());
nbrPhip[facei] *= nbrACMI.magSf()[facei]/geomArea;
}
}
forAll(nbrPhiNonOverlapp, facei)
{
const scalar w = 1.0 - cyclicACMIPolyPatch_.tgtMask()[facei];
nbrPhiNonOverlapp[facei] *= w;
}
}
// ************************************************************************* //

24
src/finiteVolume/fvMesh/fvPatches/constraint/cyclicACMI/cyclicACMIFvPatch.H
View File

@ -65,12 +65,15 @@ protected:
// Protected Member functions
//- Update the patch areas after AMI update
void updateAreas() const;
//- Helper function to reset the FV patch areas from the primitive patch
void resetPatchAreas(const fvPatch& fvp) const;
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Correct patches after moving points
virtual void movePoints();
public:
@ -134,12 +137,7 @@ public:
//- Return a reference to the AMI interpolator
virtual const AMIPatchToPatchInterpolation& AMI() const
{
const AMIPatchToPatchInterpolation& AMI =
cyclicACMIPolyPatch_.AMI();
updateAreas();
return AMI;
return cyclicACMIPolyPatch_.AMI();
}
//- Are the cyclic planes parallel
@ -169,8 +167,8 @@ public:
}
//- Return true if this patch is coupled. This is equivalent
// to the coupledPolyPatch::coupled() if parallel running or
// both sides present, false otherwise
//- to the coupledPolyPatch::coupled() if parallel running or
//- both sides present, false otherwise
virtual bool coupled() const;
//- Return delta (P to N) vectors across coupled patch
@ -183,8 +181,6 @@ public:
const Field<Type>& fld
) const
{
updateAreas();
return
cyclicACMIPolyPatch_.cyclicAMIPolyPatch::interpolate
(
@ -192,7 +188,7 @@ public:
);
}
//- Interpolate (make sure to have uptodate areas)
//- Interpolate (make sure to have up-to-date areas)
template<class Type>
tmp<Field<Type>> interpolate
(
@ -206,7 +202,7 @@ public:
// Interface transfer functions
//- Return the values of the given internal data adjacent to
// the interface as a field
//- the interface as a field
virtual tmp<labelField> interfaceInternalField
(
const labelUList& internalData

77
src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.C
View File

@ -29,6 +29,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "transform.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -82,8 +83,9 @@ void Foam::cyclicAMIFvPatch::makeWeights(scalarField& w) const
forAll(deltas, facei)
{
scalar di = deltas[facei];
scalar dni = nbrDeltas[facei];
// Note use of mag
scalar di = mag(deltas[facei]);
scalar dni = mag(nbrDeltas[facei]);
w[facei] = dni/(di + dni);
}
@ -174,4 +176,75 @@ Foam::tmp<Foam::labelField> Foam::cyclicAMIFvPatch::internalFieldTransfer
}
void Foam::cyclicAMIFvPatch::movePoints()
{
if (!owner() || !cyclicAMIPolyPatch_.createAMIFaces())
{
// Only manipulating patch face areas and mesh motion flux if the AMI
// creates additional faces
return;
}
// Update face data based on values set by the AMI manipulations
const_cast<vectorField&>(Sf()) = cyclicAMIPolyPatch_.faceAreas();
const_cast<vectorField&>(Cf()) = cyclicAMIPolyPatch_.faceCentres();
const_cast<scalarField&>(magSf()) = mag(Sf());
const cyclicAMIFvPatch& nbr = neighbPatch();
const_cast<vectorField&>(nbr.Sf()) = nbr.cyclicAMIPatch().faceAreas();
const_cast<vectorField&>(nbr.Cf()) = nbr.cyclicAMIPatch().faceCentres();
const_cast<scalarField&>(nbr.magSf()) = mag(nbr.Sf());
// Set consitent mesh motion flux
// TODO: currently maps src mesh flux to tgt - update to
// src = src + mapped(tgt) and tgt = tgt + mapped(src)?
const fvMesh& mesh = boundaryMesh().mesh();
surfaceScalarField& meshPhi = const_cast<fvMesh&>(mesh).setPhi();
surfaceScalarField::Boundary& meshPhiBf = meshPhi.boundaryFieldRef();
if (cyclicAMIPolyPatch_.owner())
{
scalarField& phip = meshPhiBf[patch().index()];
forAll(phip, facei)
{
const face& f = cyclicAMIPolyPatch_.localFaces()[facei];
// Note: using raw point locations to calculate the geometric
// area - faces areas are currently scaled by the AMI weights
// (decoupled from mesh points)
const scalar geomArea = f.mag(cyclicAMIPolyPatch_.localPoints());
const scalar scaledArea = magSf()[facei];
phip[facei] *= scaledArea/geomArea;
}
scalarField srcMeshPhi(phip);
if (Pstream::parRun())
{
AMI().srcMap().distribute(srcMeshPhi);
}
const labelListList& tgtToSrcAddr = AMI().tgtAddress();
scalarField& nbrPhip = meshPhiBf[nbr.index()];
forAll(tgtToSrcAddr, tgtFacei)
{
// Note: now have 1-to-1 mapping so tgtToSrcAddr[tgtFacei] is size 1
const label srcFacei = tgtToSrcAddr[tgtFacei][0];
nbrPhip[tgtFacei] = -srcMeshPhi[srcFacei];
}
DebugInfo
<< "patch:" << patch().name()
<< " sum(area):" << gSum(magSf())
<< " min(mag(faceAreas):" << gMin(magSf())
<< " sum(meshPhi):" << gSum(phip) << nl
<< " sum(nbrMeshPhi):" << gSum(nbrPhip) << nl
<< endl;
}
}
// ************************************************************************* //

3
src/finiteVolume/fvMesh/fvPatches/constraint/cyclicAMI/cyclicAMIFvPatch.H
View File

@ -68,6 +68,9 @@ protected:
//- Make patch weighting factors
void makeWeights(scalarField&) const;
//- Correct patches after moving points
virtual void movePoints();
public:

1
src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
View File

@ -172,7 +172,6 @@ void Foam::surfaceInterpolation::makeWeights() const
// ... and reference to the internal field of the weighting factors
scalarField& w = weights.primitiveFieldRef();
forAll(owner, facei)
{
// Note: mag in the dot-product.

142
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
View File

@ -33,6 +33,8 @@ License
#include "flipOp.H"
#include "profiling.H"
#define DEBUGAMI(msg){Pout<< "[" << __FILE__ << ":" << __LINE__ << "]: " << #msg << "=" << msg << endl;}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -49,10 +51,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolationMethodNames_
{ interpolationMethod::imPartialFaceAreaWeight, "partialFaceAreaWeightAMI" }
});
template<class SourcePatch, class TargetPatch>
bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_ =
false;
template<class SourcePatch, class TargetPatch>
template<class Patch>
Foam::tmp<Foam::scalarField>
@ -264,7 +268,6 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
}
}
// Agglomerate weights and indices
if (targetMapPtr.valid())
{
@ -298,7 +301,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
// the slots are equal to face indices.
// A mapDistribute has:
// - a subMap : these are face indices
// - a constructMap : these are from 'transferred-date' to slots
// - a constructMap : these are from 'transferred-data' to slots
labelListList tgtSubMap(Pstream::nProcs());
@ -616,10 +619,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
srcAddress_(),
srcWeights_(),
srcWeightsSum_(),
srcCentroids_(),
tgtMagSf_(),
tgtAddress_(),
tgtWeights_(),
tgtWeightsSum_(),
tgtCentroids_(),
triMode_(triMode),
srcMapPtr_(nullptr),
tgtMapPtr_(nullptr)
@ -649,10 +654,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
srcAddress_(),
srcWeights_(),
srcWeightsSum_(),
srcCentroids_(),
tgtMagSf_(),
tgtAddress_(),
tgtWeights_(),
tgtWeightsSum_(),
tgtCentroids_(),
triMode_(triMode),
srcMapPtr_(nullptr),
tgtMapPtr_(nullptr)
@ -683,10 +690,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
srcAddress_(),
srcWeights_(),
srcWeightsSum_(),
srcCentroids_(),
tgtMagSf_(),
tgtAddress_(),
tgtWeights_(),
tgtWeightsSum_(),
tgtCentroids_(),
triMode_(triMode),
srcMapPtr_(nullptr),
tgtMapPtr_(nullptr)
@ -717,10 +726,12 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
srcAddress_(),
srcWeights_(),
srcWeightsSum_(),
srcCentroids_(),
tgtMagSf_(),
tgtAddress_(),
tgtWeights_(),
tgtWeightsSum_(),
tgtCentroids_(),
triMode_(triMode),
srcMapPtr_(nullptr),
tgtMapPtr_(nullptr)
@ -834,13 +845,6 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::AMIInterpolation<SourcePatch, TargetPatch>::~AMIInterpolation()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -932,6 +936,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
(
srcAddress_,
srcWeights_,
srcCentroids_,
tgtAddress_,
tgtWeights_
);
@ -1009,9 +1014,16 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
AMIPtr->normaliseWeights(true, *this);
// Cache maps and reset addresses
List<Map<label>> cMap;
srcMapPtr_.reset(new mapDistribute(globalSrcFaces, tgtAddress_, cMap));
tgtMapPtr_.reset(new mapDistribute(globalTgtFaces, srcAddress_, cMap));
List<Map<label>> cMapSrc;
srcMapPtr_.reset
(
new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
);
List<Map<label>> cMapTgt;
tgtMapPtr_.reset
(
new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
);
}
else
{
@ -1033,6 +1045,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
(
srcAddress_,
srcWeights_,
srcCentroids_,
tgtAddress_,
tgtWeights_
);
@ -1056,6 +1069,42 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
(
autoPtr<mapDistribute>&& srcToTgtMap,
autoPtr<mapDistribute>&& tgtToSrcMap,
labelListList&& srcAddress,
scalarListList&& srcWeights,
labelListList&& tgtAddress,
scalarListList&& tgtWeights
)
{
DebugInFunction<< endl;
srcAddress_.transfer(srcAddress);
srcWeights_.transfer(srcWeights);
tgtAddress_.transfer(tgtAddress);
tgtWeights_.transfer(tgtWeights);
// Reset the sums of the weights
srcWeightsSum_.setSize(srcWeights_.size());
forAll(srcWeights_, facei)
{
srcWeightsSum_[facei] = sum(srcWeights_[facei]);
}
tgtWeightsSum_.setSize(tgtWeights_.size());
forAll(tgtWeights_, facei)
{
tgtWeightsSum_[facei] = sum(tgtWeights_[facei]);
}
srcMapPtr_ = srcToTgtMap;
tgtMapPtr_ = tgtToSrcMap;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::append
(
@ -1397,7 +1446,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::interpolateToSource
FatalErrorInFunction
<< "Employing default values when sum of weights falls below "
<< lowWeightCorrection_
<< " but supplied default field size is not equal to target "
<< " but supplied default field size is not equal to source "
<< "patch size" << nl
<< " default values = " << defaultValues.size() << nl
<< " source patch = " << srcAddress_.size() << nl
@ -1693,6 +1742,73 @@ const
}
template<class SourcePatch, class TargetPatch>
bool Foam::AMIInterpolation<SourcePatch, TargetPatch>::checkSymmetricWeights
(
const bool log
) const
{
if (Pstream::parRun() && (singlePatchProc_ == -1))
{
Log << "Checks only valid for serial running (currently)" << endl;
return true;
}
bool symmetricSrc = true;
Log << " Checking for missing src face in tgt lists" << nl;
forAll(srcAddress_, srcFacei)
{
const labelList& tgtIds = srcAddress_[srcFacei];
for (const label tgtFacei : tgtIds)
{
if (!tgtAddress_[tgtFacei].found(srcFacei))
{
symmetricSrc = false;
Log << " srcFacei:" << srcFacei
<< " not found in tgtToSrc list for tgtFacei:"
<< tgtFacei << nl;
}
}
}
if (symmetricSrc)
{
Log << " - symmetric" << endl;
}
bool symmetricTgt = true;
Log << " Checking for missing tgt face in src lists" << nl;
forAll(tgtAddress_, tgtFacei)
{
const labelList& srcIds = tgtAddress_[tgtFacei];
for (const label srcFacei : srcIds)
{
if (!srcAddress_[srcFacei].found(tgtFacei))
{
symmetricTgt = false;
Log << " tgtFacei:" << tgtFacei
<< " not found in srcToTgt list for srcFacei:"
<< srcFacei << nl;
}
}
}
if (symmetricTgt)
{
Log << " - symmetric" << endl;
}
return symmetricSrc && symmetricTgt;
}
template<class SourcePatch, class TargetPatch>
void Foam::AMIInterpolation<SourcePatch, TargetPatch>::writeFaceConnectivity
(

35
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
View File

@ -61,6 +61,7 @@ SourceFiles
#include "globalIndex.H"
#include "ops.H"
#include "Enum.H"
#include "pointList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -146,6 +147,9 @@ private:
//- Sum of weights of target faces per source face
scalarField srcWeightsSum_;
//- Centroid of target faces per source face
pointListList srcCentroids_;
// Target patch
@ -161,6 +165,9 @@ private:
//- Sum of weights of source faces per target face
scalarField tgtWeightsSum_;
//- Centroid of source faces per target face
pointListList tgtCentroids_;
//- Face triangulation mode
const faceAreaIntersect::triangulationMode triMode_;
@ -347,7 +354,7 @@ public:
//- Destructor
~AMIInterpolation();
~AMIInterpolation() = default;
// Typedef to SourcePatch type this AMIInterpolation is instantiated on
typedef SourcePatch sourcePatchType;
@ -399,6 +406,12 @@ public:
//- patch weights (i.e. the sum before normalisation)
inline scalarField& srcWeightsSum();
//- Return const access to source patch face centroids
inline const pointListList& srcCentroids() const;
//- Return access to source patch face centroids
inline pointListList& srcCentroids();
//- Source map pointer - valid only if singlePatchProc = -1
//- This gets source data into a form to be consumed by
//- tgtAddress, tgtWeights
@ -448,6 +461,22 @@ public:
const TargetPatch& tgtPatch
);
void update
(
autoPtr<mapDistribute>&& srcToTgtMap,
autoPtr<mapDistribute>&& tgtToSrcMap,
labelListList&& srcAddress,
scalarListList&& srcWeights,
labelListList&& tgtAddress,
scalarListList&& tgtWeights
);
void setAreas(const scalarList& srcMagSf, const scalarList& tgtMagSf)
{
srcMagSf_ = srcMagSf;
tgtMagSf_ = tgtMagSf;
}
//- Append additional addressing and weights
void append
(
@ -582,6 +611,10 @@ public:
// Checks
//- Check if src addresses are present in tgt addresses and
//- viceversa
bool checkSymmetricWeights(const bool log) const;
//- Write face connectivity as OBJ file
void writeFaceConnectivity
(

16
src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
View File

@ -115,6 +115,22 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcWeightsSum()
}
template<class SourcePatch, class TargetPatch>
inline const Foam::pointListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcCentroids() const
{
return srcCentroids_;
}
template<class SourcePatch, class TargetPatch>
inline Foam::pointListList&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcCentroids()
{
return srcCentroids_;
}
template<class SourcePatch, class TargetPatch>
inline const Foam::mapDistribute&
Foam::AMIInterpolation<SourcePatch, TargetPatch>::srcMap() const

18
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.C
View File

@ -31,6 +31,8 @@ License
#include "mapDistribute.H"
#include "unitConversion.H"
#include "findNearestMaskedOp.H"
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -231,17 +233,25 @@ void Foam::AMIMethod<SourcePatch, TargetPatch>::resetTree()
template<class SourcePatch, class TargetPatch>
Foam::label Foam::AMIMethod<SourcePatch, TargetPatch>::findTargetFace
(
const label srcFacei
const label srcFacei,
const UList<label>& excludeFaces,
const label srcFacePti
) const
{
label targetFacei = -1;
const pointField& srcPts = srcPatch_.points();
const face& srcFace = srcPatch_[srcFacei];
const point srcPt = srcFace.centre(srcPts);
const scalar srcFaceArea = srcMagSf_[srcFacei];
pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
findNearestMaskedOp<TargetPatch> fnOp(*treePtr_, excludeFaces);
boundBox bb(srcPts, srcFace, false);
const point srcPt =
srcFacePti == -1 ? bb.centre() : srcPts[srcFace[srcFacePti]];
const pointIndexHit sample =
treePtr_->findNearest(srcPt, magSqr(bb.max() - bb.centre()), fnOp);
if (sample.hit())
{

12
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/AMIMethod/AMIMethod.H
View File

@ -45,6 +45,7 @@ SourceFiles
#include "treeDataPrimitivePatch.H"
#include "treeBoundBoxList.H"
#include "runTimeSelectionTables.H"
#include "pointList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -101,7 +102,7 @@ protected:
List<scalar> tgtMagSf_;
//- Labels of faces that are not overlapped by any target faces
//- (should be empty for correct functioning)
//- (should be empty for correct functioning for fully covered AMIs)
labelList srcNonOverlap_;
//- Octree used to find face seeds
@ -145,8 +146,12 @@ protected:
//- Reset the octree for the target patch face search
void resetTree();
//- Find face on target patch that overlaps source face
label findTargetFace(const label srcFacei) const;
label findTargetFace
(
const label srcFacei,
const UList<label>& excludeFaces = UList<label>::null(),
const label srcFacePti = -1
) const;
//- Add faces neighbouring facei to the ID list
void appendNbrFaces
@ -251,6 +256,7 @@ public:
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei = -1,

8
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
View File

@ -205,13 +205,6 @@ Foam::directAMI<SourcePatch, TargetPatch>::directAMI
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::directAMI<SourcePatch, TargetPatch>::~directAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -219,6 +212,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei,

3
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.H
View File

@ -118,7 +118,7 @@ public:
//- Destructor
virtual ~directAMI();
virtual ~directAMI() = default;
// Member Functions
@ -130,6 +130,7 @@ public:
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei = -1,

300
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.C
View File

@ -28,6 +28,7 @@ License
#include "faceAreaWeightAMI.H"
#include "profiling.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
@ -36,6 +37,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
(
List<DynamicList<label>>& srcAddr,
List<DynamicList<scalar>>& srcWght,
List<DynamicList<point>>& srcCtr,
List<DynamicList<label>>& tgtAddr,
List<DynamicList<scalar>>& tgtWght,
label srcFacei,
@ -60,7 +62,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
seedFaces[srcFacei] = tgtFacei;
// list to keep track of whether src face can be mapped
boolList mapFlag(nFacesRemaining, true);
bitSet mapFlag(nFacesRemaining, true);
// reset starting seed
label startSeedi = 0;
@ -68,6 +70,9 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
DynamicList<label> nonOverlapFaces;
do
{
nbrFaces.clear();
visitedFaces.clear();
// Do advancing front starting from srcFacei,tgtFacei
bool faceProcessed = processSourceFace
(
@ -79,11 +84,12 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcAddressing
srcAddr,
srcWght,
srcCtr,
tgtAddr,
tgtWght
);
mapFlag[srcFacei] = false;
mapFlag.unset(srcFacei);
nFacesRemaining--;
@ -122,9 +128,10 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
// list of faces currently visited for srcFacei to avoid multiple hits
DynamicList<label>& visitedFaces,
// temporary storage for addressing and weights
// temporary storage for addressing, weights and centroid
List<DynamicList<label>>& srcAddr,
List<DynamicList<scalar>>& srcWght,
List<DynamicList<point>>& srcCtr,
List<DynamicList<label>>& tgtAddr,
List<DynamicList<scalar>>& tgtWght
)
@ -136,9 +143,6 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
return false;
}
nbrFaces.clear();
visitedFaces.clear();
// append initial target face and neighbours
nbrFaces.append(tgtStartFacei);
this->appendNbrFaces
@ -158,16 +162,24 @@ bool Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::processSourceFace
// process new target face
label tgtFacei = nbrFaces.remove();
visitedFaces.append(tgtFacei);
scalar area = interArea(srcFacei, tgtFacei);
scalar interArea = 0;
vector interCentroid(Zero);
calcInterArea(srcFacei, tgtFacei, interArea, interCentroid);
// store when intersection fractional area > tolerance
if (area/this->srcMagSf_[srcFacei] > faceAreaIntersect::tolerance())
if
(
interArea/this->srcMagSf_[srcFacei]
> faceAreaIntersect::tolerance()
)
{
srcAddr[srcFacei].append(tgtFacei);
srcWght[srcFacei].append(area);
srcWght[srcFacei].append(interArea);
srcCtr[srcFacei].append(interCentroid);
tgtAddr[tgtFacei].append(srcFacei);
tgtWght[tgtFacei].append(area);
tgtWght[tgtFacei].append(interArea);
this->appendNbrFaces
(
@ -199,10 +211,10 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
label& startSeedi,
label& srcFacei,
label& tgtFacei,
const boolList& mapFlag,
const bitSet& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
bool errorOnNotFound
const bool errorOnNotFound
) const
{
addProfiling(ami, "faceAreaWeightAMI::setNextFaces");
@ -216,15 +228,14 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
bool valuesSet = false;
for (label faceS: srcNbrFaces)
{
if (mapFlag[faceS] && seedFaces[faceS] == -1)
if (mapFlag.test(faceS) && seedFaces[faceS] == -1)
{
for (label faceT : visitedFaces)
{
const scalar threshold =
this->srcMagSf_[faceS]*faceAreaIntersect::tolerance();
// store when intersection fractional area > tolerance
// Check that faces have enough overlap for robust walking
// Store when intersection fractional area > threshold
if (overlaps(faceS, faceT, threshold))
{
seedFaces[faceS] = faceT;
@ -248,25 +259,31 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
else
{
// try to use existing seed
bool foundNextSeed = false;
for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
label facei = startSeedi;
if (!mapFlag.test(startSeedi))
{
if (mapFlag[facei])
facei = mapFlag.find_next(facei);
}
const label startSeedi0 = facei;
bool foundNextSeed = false;
while (facei != -1)
{
if (!foundNextSeed)
{
if (!foundNextSeed)
{
startSeedi = facei;
foundNextSeed = true;
}
if (seedFaces[facei] != -1)
{
srcFacei = facei;
tgtFacei = seedFaces[facei];
return;
}
startSeedi = facei;
foundNextSeed = true;
}
if (seedFaces[facei] != -1)
{
srcFacei = facei;
tgtFacei = seedFaces[facei];
return;
}
facei = mapFlag.find_next(facei);
}
// perform new search to find match
@ -276,25 +293,18 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
<< "target face" << endl;
}
foundNextSeed = false;
for (label facei = startSeedi; facei < mapFlag.size(); ++facei)
facei = startSeedi0;
while (facei != -1)
{
if (mapFlag[facei])
srcFacei = facei;
tgtFacei = this->findTargetFace(srcFacei, visitedFaces);
if (tgtFacei >= 0)
{
if (!foundNextSeed)
{
startSeedi = facei + 1;
foundNextSeed = true;
}
srcFacei = facei;
tgtFacei = this->findTargetFace(srcFacei);
if (tgtFacei >= 0)
{
return;
}
return;
}
facei = mapFlag.find_next(facei);
}
if (errorOnNotFound)
@ -307,16 +317,16 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
template<class SourcePatch, class TargetPatch>
Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calcInterArea
(
const label srcFacei,
const label tgtFacei
const label tgtFacei,
scalar& area,
vector& centroid
) const
{
addProfiling(ami, "faceAreaWeightAMI::interArea");
scalar area = 0;
const pointField& srcPoints = this->srcPatch_.points();
const pointField& tgtPoints = this->tgtPatch_.points();
@ -331,7 +341,7 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
|| (this->tgtMagSf_[tgtFacei] < ROOTVSMALL)
)
{
return area;
return;
}
// create intersection object
@ -359,7 +369,7 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
if (magN > ROOTVSMALL)
{
area = inter.calc(src, tgt, n/magN);
inter.calc(src, tgt, n/magN, area, centroid);
}
else
{
@ -371,13 +381,24 @@ Foam::scalar Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::interArea
<< endl;
}
if
(
AMIInterpolation<SourcePatch, TargetPatch>::cacheIntersections_
&& debug
)
{
static OBJstream tris("intersectionTris.obj");
const auto& triPts = inter.triangles();
for (const auto& tp : triPts)
{
tris.write(triPointRef(tp[0], tp[1], tp[2]), false);
}
}
if ((debug > 1) && (area > 0))
{
this->writeIntersectionOBJ(area, src, tgt, srcPoints, tgtPoints);
}
return area;
}
@ -452,95 +473,70 @@ restartUncoveredSourceFace
(
List<DynamicList<label>>& srcAddr,
List<DynamicList<scalar>>& srcWght,
List<DynamicList<point>>& srcCtr,
List<DynamicList<label>>& tgtAddr,
List<DynamicList<scalar>>& tgtWght
)
{
addProfiling(ami, "faceAreaWeightAMI::restartUncoveredSourceFace");
// Collect all src faces with a low weight
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: exclude faces in srcNonOverlap_ for ACMI?
label nBelowMinWeight = 0;
const scalar minWeight = 0.95;
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFacei to avoid multiple hits
DynamicList<label> visitedFaces(10);
labelHashSet lowWeightFaces(100);
forAll(srcWght, srcFacei)
{
scalar s = sum(srcWght[srcFacei]);
scalar t = s/this->srcMagSf_[srcFacei];
const scalar s = sum(srcWght[srcFacei]);
const scalar t = s/this->srcMagSf_[srcFacei];
if (t < 0.5)
if (t < minWeight)
{
lowWeightFaces.insert(srcFacei);
}
}
++nBelowMinWeight;
if (debug)
{
Pout<< "faceAreaWeightAMI: restarting search on "
<< lowWeightFaces.size() << " faces since sum of weights < 0.5"
<< endl;
}
const face& f = this->srcPatch_[srcFacei];
if (lowWeightFaces.size() > 0)
{
// Erase all the lowWeight source faces from the target
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DynamicList<label> okSrcFaces(10);
DynamicList<scalar> okSrcWeights(10);
forAll(tgtAddr, tgtFacei)
{
okSrcFaces.clear();
okSrcWeights.clear();
DynamicList<label>& srcFaces = tgtAddr[tgtFacei];
DynamicList<scalar>& srcWeights = tgtWght[tgtFacei];
forAll(srcFaces, i)
forAll(f, fpi)
{
if (!lowWeightFaces.found(srcFaces[i]))
const label tgtFacei =
this->findTargetFace(srcFacei, srcAddr[srcFacei], fpi);
if (tgtFacei != -1)
{
okSrcFaces.append(srcFaces[i]);
okSrcWeights.append(srcWeights[i]);
nbrFaces.clear();
visitedFaces = srcAddr[srcFacei];
(void)processSourceFace
(
srcFacei,
tgtFacei,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
srcCtr,
tgtAddr,
tgtWght
);
}
}
if (okSrcFaces.size() < srcFaces.size())
{
srcFaces.transfer(okSrcFaces);
srcWeights.transfer(okSrcWeights);
}
}
}
// Restart search from best hit
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// list of tgt face neighbour faces
DynamicList<label> nbrFaces(10);
// list of faces currently visited for srcFacei to avoid multiple hits
DynamicList<label> visitedFaces(10);
for (const label srcFacei : lowWeightFaces)
{
const label tgtFacei = this->findTargetFace(srcFacei);
if (tgtFacei != -1)
{
//bool faceProcessed =
processSourceFace
(
srcFacei,
tgtFacei,
nbrFaces,
visitedFaces,
srcAddr,
srcWght,
tgtAddr,
tgtWght
);
// ? Check faceProcessed to see if restarting has worked.
}
}
if (debug && nBelowMinWeight)
{
WarningInFunction
<< "Restarted search on " << nBelowMinWeight
<< " faces since sum of weights < " << minWeight
<< endl;
}
}
@ -572,16 +568,45 @@ Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::faceAreaWeightAMI
{
this->triangulatePatch(srcPatch, srcTris_, this->srcMagSf_);
this->triangulatePatch(tgtPatch, tgtTris_, this->tgtMagSf_);
if (debug)
{
static label nAMI = 0;
// Write out triangulated surfaces as OBJ files
OBJstream srcTriObj("srcTris_" + Foam::name(nAMI) + ".obj");
const pointField& srcPts = srcPatch.points();
forAll(srcTris_, facei)
{
const DynamicList<face>& faces = srcTris_[facei];
for (const face& f : faces)
{
srcTriObj.write
(
triPointRef(srcPts[f[0]], srcPts[f[1]], srcPts[f[2]])
);
}
}
OBJstream tgtTriObj("tgtTris_" + Foam::name(nAMI) + ".obj");
const pointField& tgtPts = tgtPatch.points();
forAll(tgtTris_, facei)
{
const DynamicList<face>& faces = tgtTris_[facei];
for (const face& f : faces)
{
tgtTriObj.write
(
triPointRef(tgtPts[f[0]], tgtPts[f[1]], tgtPts[f[2]])
);
}
}
++nAMI;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::~faceAreaWeightAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -589,6 +614,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei,
@ -608,6 +634,9 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
tgtFacei
);
srcCentroids.setSize(srcAddress.size());
if (!ok)
{
return;
@ -616,6 +645,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
// temporary storage for addressing and weights
List<DynamicList<label>> srcAddr(this->srcPatch_.size());
List<DynamicList<scalar>> srcWght(srcAddr.size());
List<DynamicList<point>> srcCtr(srcAddr.size());
List<DynamicList<label>> tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar>> tgtWght(tgtAddr.size());
@ -623,6 +653,7 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
srcAddr,
srcWght,
srcCtr,
tgtAddr,
tgtWght,
srcFacei,
@ -644,17 +675,18 @@ void Foam::faceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
srcAddr,
srcWght,
srcCtr,
tgtAddr,
tgtWght
);
}
// transfer data to persistent storage
// Transfer data to persistent storage
forAll(srcAddr, i)
{
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i].transfer(srcWght[i]);
srcCentroids[i].transfer(srcCtr[i]);
}
forAll(tgtAddr, i)
{

21
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/faceAreaWeightAMI/faceAreaWeightAMI.H
View File

@ -80,11 +80,13 @@ protected:
// Marching front
//- Calculate addressing and weights using temporary storage
//- Calculate addressing, weights and centroids using temporary
//- storage
virtual void calcAddressing
(
List<DynamicList<label>>& srcAddress,
List<DynamicList<scalar>>& srcWeights,
List<DynamicList<point>>& srcCentroids,
List<DynamicList<label>>& tgtAddress,
List<DynamicList<scalar>>& tgtWeights,
label srcFacei,
@ -100,6 +102,7 @@ protected:
DynamicList<label>& visitedFaces,
List<DynamicList<label>>& srcAddr,
List<DynamicList<scalar>>& srcWght,
List<DynamicList<point>>& srcCtr,
List<DynamicList<label>>& tgtAddr,
List<DynamicList<scalar>>& tgtWght
);
@ -109,6 +112,7 @@ protected:
(
List<DynamicList<label>>& srcAddr,
List<DynamicList<scalar>>& srcWght,
List<DynamicList<point>>& srcCtr,
List<DynamicList<label>>& tgtAddr,
List<DynamicList<scalar>>& tgtWght
);
@ -119,20 +123,22 @@ protected:
label& startSeedi,
label& srcFacei,
label& tgtFacei,
const boolList& mapFlag,
const bitSet& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
bool errorOnNotFound = true
const bool errorOnNotFound = true
) const;
// Evaluation
//- Area of intersection between source and target faces
virtual scalar interArea
virtual void calcInterArea
(
const label srcFacei,
const label tgtFacei
const label tgtFacei,
scalar& area,
vector& centroid
) const;
//- Return true if faces overlap
@ -165,18 +171,19 @@ public:
//- Destructor
virtual ~faceAreaWeightAMI();
virtual ~faceAreaWeightAMI() = default;
// Member Functions
// Manipulation
//- Update addressing and weights
//- Update addressing, weights and centroids
virtual void calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei = -1,

8
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.C
View File

@ -190,13 +190,6 @@ Foam::mapNearestAMI<SourcePatch, TargetPatch>::mapNearestAMI
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::mapNearestAMI<SourcePatch, TargetPatch>::~mapNearestAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -204,6 +197,7 @@ void Foam::mapNearestAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei,

3
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/mapNearestAMI/mapNearestAMI.H
View File

@ -123,7 +123,7 @@ public:
//- Destructor
virtual ~mapNearestAMI();
virtual ~mapNearestAMI() = default;
// Member Functions
@ -135,6 +135,7 @@ public:
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei = -1,

16
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
View File

@ -35,7 +35,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::setNextFaces
label& startSeedi,
label& srcFacei,
label& tgtFacei,
const boolList& mapFlag,
const bitSet& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
const bool errorOnNotFound
@ -78,14 +78,6 @@ partialFaceAreaWeightAMI
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::
~partialFaceAreaWeightAMI()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -100,6 +92,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei,
@ -122,9 +115,12 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
return;
}
srcCentroids.setSize(srcAddress.size());
// temporary storage for addressing and weights
List<DynamicList<label>> srcAddr(this->srcPatch_.size());
List<DynamicList<scalar>> srcWght(srcAddr.size());
List<DynamicList<point>> srcCtr(srcAddr.size());
List<DynamicList<label>> tgtAddr(this->tgtPatch_.size());
List<DynamicList<scalar>> tgtWght(tgtAddr.size());
@ -132,6 +128,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
(
srcAddr,
srcWght,
srcCtr,
tgtAddr,
tgtWght,
srcFacei,
@ -143,6 +140,7 @@ void Foam::partialFaceAreaWeightAMI<SourcePatch, TargetPatch>::calculate
{
srcAddress[i].transfer(srcAddr[i]);
srcWeights[i].transfer(srcWght[i]);
srcCentroids[i].transfer(srcCtr[i]);
}
forAll(tgtAddr, i)
{

30
src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
View File

@ -65,19 +65,22 @@ private:
//- No copy assignment
void operator=(const partialFaceAreaWeightAMI&) = delete;
// Marching front
//- Set the source and target seed faces
virtual void setNextFaces
(
label& startSeedi,
label& srcFacei,
label& tgtFacei,
const boolList& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
bool errorOnNotFound = true
) const;
protected:
// Protected Member Functions
//- Set the source and target seed faces
virtual void setNextFaces
(
label& startSeedi,
label& srcFacei,
label& tgtFacei,
const bitSet& mapFlag,
labelList& seedFaces,
const DynamicList<label>& visitedFaces,
const bool errorOnNotFound = true
) const;
public:
@ -100,7 +103,7 @@ public:
//- Destructor
virtual ~partialFaceAreaWeightAMI();
virtual ~partialFaceAreaWeightAMI() = default;
// Member Functions
@ -118,6 +121,7 @@ public:
(
labelListList& srcAddress,
scalarListList& srcWeights,
pointListList& srcCentroids,
labelListList& tgtAddress,
scalarListList& tgtWeights,
label srcFacei = -1,

62
src/meshTools/AMIInterpolation/AMIInterpolation/findNearestMaskedOp.H Normal file
View File

@ -0,0 +1,62 @@
#include "indexedOctree.H"
#include "treeDataPrimitivePatch.H"
namespace Foam
{
template<class PatchType>
class findNearestMaskedOp
{
const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree_;
const labelUList& excludeIndices_;
public:
findNearestMaskedOp
(
const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
const labelUList& excludeIndices
)
:
tree_(tree),
excludeIndices_(excludeIndices)
{}
void operator()
(
const labelUList& indices,
const point& sample,
scalar& nearestDistSqr,
label& minIndex,
point& nearestPoint
) const
{
const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
const PatchType& patch = shape.patch();
const pointField& points = patch.points();
forAll(indices, i)
{
const label index = indices[i];
if (!excludeIndices_.found(index))
{
const typename PatchType::FaceType& f = patch[index];
pointHit nearHit = f.nearestPoint(sample, points);
scalar distSqr = sqr(nearHit.distance());
if (distSqr < nearestDistSqr)
{
nearestDistSqr = distSqr;
minIndex = index;
nearestPoint = nearHit.rawPoint();
}
}
}
}
};
} // End namespace Foam

135
src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
View File

@ -219,13 +219,15 @@ void Foam::faceAreaIntersect::triSliceWithPlane
}
Foam::scalar Foam::faceAreaIntersect::triangleIntersect
void Foam::faceAreaIntersect::triangleIntersect
(
const triPoints& src,
const point& tgt0,
const point& tgt1,
const point& tgt2,
const vector& n
const vector& n,
scalar& area,
vector& centroid
) const
{
// Work storage
@ -241,7 +243,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
const scalar srcArea(triArea(src));
if (srcArea < ROOTVSMALL)
{
return 0.0;
return;
}
// Typical length scale
@ -255,7 +257,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
scalar s = mag(tgt1 - tgt0);
if (s < ROOTVSMALL)
{
return 0.0;
return;
}
// Note: outer product with n pre-scaled with edge length. This is
@ -268,7 +270,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
// Triangle either zero edge length (s == 0) or
// perpendicular to face normal n. In either case zero
// overlap area
return 0.0;
return;
}
else
{
@ -279,7 +281,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
if (nWorkTris1 == 0)
{
return 0.0;
return;
}
// Edge 1
@ -290,7 +292,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
scalar s = mag(tgt2 - tgt1);
if (s < ROOTVSMALL)
{
return 0.0;
return;
}
const vector n1((tgt1 - tgt2)^(-s*n));
const scalar magSqrN1(magSqr(n1));
@ -300,7 +302,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
// Triangle either zero edge length (s == 0) or
// perpendicular to face normal n. In either case zero
// overlap area
return 0.0;
return;
}
else
{
@ -315,7 +317,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
if (nWorkTris2 == 0)
{
return 0.0;
return;
}
}
}
@ -328,7 +330,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
scalar s = mag(tgt2 - tgt0);
if (s < ROOTVSMALL)
{
return 0.0;
return;
}
const vector n2((tgt2 - tgt0)^(-s*n));
const scalar magSqrN2(magSqr(n2));
@ -338,7 +340,7 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
// Triangle either zero edge length (s == 0) or
// perpendicular to face normal n. In either case zero
// overlap area
return 0.0;
return;
}
else
{
@ -353,23 +355,25 @@ Foam::scalar Foam::faceAreaIntersect::triangleIntersect
if (nWorkTris1 == 0)
{
return 0.0;
return;
}
else
{
// Calculate area of sub-triangles
scalar area = 0.0;
for (label i = 0; i < nWorkTris1; ++i)
{
area += triArea(workTris1[i]);
// Area of intersection
const scalar currArea = triArea(workTris1[i]);
area += currArea;
// Area-weighted centroid of intersection
centroid += currArea*triCentroid(workTris1[i]);
if (cacheTriangulation_)
{
triangles_.append(workTris1[i]);
}
}
return area;
}
}
}
@ -443,12 +447,13 @@ void Foam::faceAreaIntersect::triangulate
}
Foam::scalar Foam::faceAreaIntersect::calc
void Foam::faceAreaIntersect::calc
(
const face& faceA,
const face& faceB,
const vector& n
const vector& n,
scalar& area,
vector& centroid
) const
{
if (cacheTriangulation_)
@ -456,8 +461,10 @@ Foam::scalar Foam::faceAreaIntersect::calc
triangles_.clear();
}
area = 0.0;
centroid = vector::zero;
// Intersect triangles
scalar totalArea = 0.0;
for (const face& triA : trisA_)
{
triPoints tpA = getTriPoints(pointsA_, triA, false);
@ -466,32 +473,38 @@ Foam::scalar Foam::faceAreaIntersect::calc
{
if (reverseB_)
{
totalArea +=
triangleIntersect
(
tpA,
pointsB_[triB[0]],
pointsB_[triB[1]],
pointsB_[triB[2]],
n
);
triangleIntersect
(
tpA,
pointsB_[triB[0]],
pointsB_[triB[1]],
pointsB_[triB[2]],
n,
area,
centroid
);
}
else
{
totalArea +=
triangleIntersect
(
tpA,
pointsB_[triB[2]],
pointsB_[triB[1]],
pointsB_[triB[0]],
n
);
triangleIntersect
(
tpA,
pointsB_[triB[2]],
pointsB_[triB[1]],
pointsB_[triB[0]],
n,
area,
centroid
);
}
}
}
return totalArea;
// Area weighed centroid
if (area > 0)
{
centroid /= area;
}
}
@ -503,8 +516,10 @@ bool Foam::faceAreaIntersect::overlaps
const scalar threshold
) const
{
scalar area = 0.0;
vector centroid(Zero);
// Intersect triangles
scalar totalArea = 0.0;
for (const face& triA : trisA_)
{
const triPoints tpA = getTriPoints(pointsA_, triA, false);
@ -513,30 +528,32 @@ bool Foam::faceAreaIntersect::overlaps
{
if (reverseB_)
{
totalArea +=
triangleIntersect
(
tpA,
pointsB_[triB[0]],
pointsB_[triB[1]],
pointsB_[triB[2]],
n
);
triangleIntersect
(
tpA,
pointsB_[triB[0]],
pointsB_[triB[1]],
pointsB_[triB[2]],
n,
area,
centroid
);
}
else
{
totalArea +=
triangleIntersect
(
tpA,
pointsB_[triB[2]],
pointsB_[triB[1]],
pointsB_[triB[0]],
n
);
triangleIntersect
(
tpA,
pointsB_[triB[2]],
pointsB_[triB[1]],
pointsB_[triB[0]],
n,
area,
centroid
);
}
if (totalArea > threshold)
if (area > threshold)
{
return true;
}

19
src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
View File

@ -97,6 +97,7 @@ private:
// Static data members
//- Tolerance
static scalar tol;
@ -132,6 +133,9 @@ private:
//- Return triangle area
inline scalar triArea(const triPoints& t) const;
//- Return triangle centre
inline vector triCentroid(const triPoints& t) const;
//- Slice triangle with plane and generate new cut sub-triangles
void triSliceWithPlane
(
@ -143,13 +147,15 @@ private:
) const;
//- Return area of intersection of triangles src and tgt
scalar triangleIntersect
void triangleIntersect
(
const triPoints& src,
const point& tgt0,
const point& tgt1,
const point& tgt2,
const vector& n
const vector& n,
scalar& area,
vector& centroid
) const;
@ -199,12 +205,15 @@ public:
DynamicList<face>& faces
);
//- Return area of intersection of faceA with faceB
scalar calc
//- Return area of intersection of faceA with faceB and effective
//- face centre
void calc
(
const face& faceA,
const face& faceB,
const vector& n
const vector& n,
scalar& area,
vector& centroid
) const;
//- Return area of intersection of faceA with faceB

9
src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
View File

@ -113,6 +113,15 @@ inline Foam::scalar Foam::faceAreaIntersect::triArea(const triPoints& t) const
}
inline Foam::vector Foam::faceAreaIntersect::triCentroid
(
const triPoints& t
) const
{
return (t[0] + t[1] + t[2])/3;
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar& Foam::faceAreaIntersect::tolerance()

379
src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
Copyright (C) 2017-2018 OpenCFD Ltd.
Copyright (C) 2017-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -41,195 +41,225 @@ namespace Foam
addToRunTimeSelectionTable(polyPatch, cyclicACMIPolyPatch, dictionary);
}
const Foam::scalar Foam::cyclicACMIPolyPatch::tolerance_ = 1e-10;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
void Foam::cyclicACMIPolyPatch::reportCoverage
(
const word& name,
const scalarField& weightSum
) const
{
label nUncovered = 0;
label nCovered = 0;
for (const scalar sum : weightSum)
{
if (sum < tolerance_)
{
++nUncovered;
}
else if (sum > scalar(1) - tolerance_)
{
++nCovered;
}
}
reduce(nUncovered, sumOp<label>());
reduce(nCovered, sumOp<label>());
label nTotal = returnReduce(weightSum.size(), sumOp<label>());
Info<< "ACMI: Patch " << name << " uncovered/blended/covered = "
<< nUncovered << ", " << nTotal-nUncovered-nCovered
<< ", " << nCovered << endl;
}
void Foam::cyclicACMIPolyPatch::resetAMI
(
const AMIPatchToPatchInterpolation::interpolationMethod&
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
) const
{
if (owner())
{
const polyPatch& nonOverlapPatch = this->nonOverlapPatch();
resetAMI(boundaryMesh().mesh().points(), AMIMethod);
}
void Foam::cyclicACMIPolyPatch::resetAMI
(
const UList<point>& points,
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
) const
{
if (!owner())
{
return;
}
const polyPatch& nonOverlapPatch = this->nonOverlapPatch();
if (debug)
{
Pout<< "cyclicACMIPolyPatch::resetAMI : recalculating weights"
<< " for " << name() << " and " << nonOverlapPatch.name()
<< endl;
}
WarningInFunction<< "DEACTIVATED clearGeom()" << endl;
// if (boundaryMesh().mesh().hasCellCentres())
if (0 && boundaryMesh().mesh().hasCellCentres())
{
if (debug)
{
Pout<< "cyclicACMIPolyPatch::resetAMI : recalculating weights"
Pout<< "cyclicACMIPolyPatch::resetAMI : clearing cellCentres"
<< " for " << name() << " and " << nonOverlapPatch.name()
<< endl;
}
if (boundaryMesh().mesh().hasCellCentres())
{
if (debug)
{
Pout<< "cyclicACMIPolyPatch::resetAMI : clearing cellCentres"
<< " for " << name() << " and " << nonOverlapPatch.name()
<< endl;
}
//WarningInFunction
// << "The mesh already has cellCentres calculated when"
// << " resetting ACMI " << name() << "." << endl
// << "This is a problem since ACMI adapts the face areas"
// << " (to close cells) so this has" << endl
// << "to be done before cell centre calculation." << endl
// << "This can happen if e.g. the cyclicACMI is after"
// << " any processor patches in the boundary." << endl;
const_cast<polyMesh&>
(
boundaryMesh().mesh()
).primitiveMesh::clearGeom();
}
// Trigger re-building of faceAreas
(void)boundaryMesh().mesh().faceAreas();
// Calculate the AMI using partial face-area-weighted. This leaves
// the weights as fractions of local areas (sum(weights) = 1 means
// face is fully covered)
cyclicAMIPolyPatch::resetAMI
//WarningInFunction
// << "The mesh already has cellCentres calculated when"
// << " resetting ACMI " << name() << "." << endl
// << "This is a problem since ACMI adapts the face areas"
// << " (to close cells) so this has" << endl
// << "to be done before cell centre calculation." << endl
// << "This can happen if e.g. the cyclicACMI is after"
// << " any processor patches in the boundary." << endl;
const_cast<polyMesh&>
(
AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
);
boundaryMesh().mesh()
).primitiveMesh::clearGeom();
}
AMIPatchToPatchInterpolation& AMI =
const_cast<AMIPatchToPatchInterpolation&>(this->AMI());
// Output some stats. AMIInterpolation will have already output the
// average weights ("sum(weights) min:1 max:1 average:1")
// Trigger re-building of faceAreas
(void)boundaryMesh().mesh().faceAreas();
// Calculate the AMI using partial face-area-weighted. This leaves
// the weights as fractions of local areas (sum(weights) = 1 means
// face is fully covered)
cyclicAMIPolyPatch::resetAMI
(
points,
AMIPatchToPatchInterpolation::imPartialFaceAreaWeight
);
const AMIPatchToPatchInterpolation& AMI = this->AMI();
// Output some statistics
reportCoverage("source", AMI.srcWeightsSum());
reportCoverage("target", AMI.tgtWeightsSum());
// Set the mask fields
// Note:
// - assumes that the non-overlap patches are decomposed using the same
// decomposition as the coupled patches (per side)
srcMask_ = min(scalar(1), max(scalar(0), AMI.srcWeightsSum()));
tgtMask_ = min(scalar(1), max(scalar(0), AMI.tgtWeightsSum()));
if (debug)
{
Pout<< "resetAMI" << endl;
{
const scalarField& wghtsSum = AMI.srcWeightsSum();
label nUncovered = 0;
label nCovered = 0;
forAll(wghtsSum, facei)
{
scalar sum = wghtsSum[facei];
if (sum < tolerance_)
{
nUncovered++;
}
else if (sum > scalar(1)-tolerance_)
{
nCovered++;
}
}
reduce(nUncovered, sumOp<label>());
reduce(nCovered, sumOp<label>());
label nTotal = returnReduce(wghtsSum.size(), sumOp<label>());
Info<< "ACMI: Patch source uncovered/blended/covered = "
<< nUncovered << ", " << nTotal-nUncovered-nCovered
<< ", " << nCovered << endl;
const cyclicACMIPolyPatch& patch = *this;
Pout<< "patch:" << patch.name() << " size:" << patch.size()
<< " non-overlap patch: " << patch.nonOverlapPatch().name()
<< " size:" << patch.nonOverlapPatch().size()
<< " mask size:" << patch.srcMask().size() << endl;
}
{
const scalarField& wghtsSum = AMI.tgtWeightsSum();
label nUncovered = 0;
label nCovered = 0;
forAll(wghtsSum, facei)
{
scalar sum = wghtsSum[facei];
if (sum < tolerance_)
{
nUncovered++;
}
else if (sum > scalar(1)-tolerance_)
{
nCovered++;
}
}
reduce(nUncovered, sumOp<label>());
reduce(nCovered, sumOp<label>());
label nTotal = returnReduce(wghtsSum.size(), sumOp<label>());
Info<< "ACMI: Patch target uncovered/blended/covered = "
<< nUncovered << ", " << nTotal-nUncovered-nCovered
<< ", " << nCovered << endl;
const cyclicACMIPolyPatch& patch = this->neighbPatch();
Pout<< "patch:" << patch.name() << " size:" << patch.size()
<< " non-overlap patch: " << patch.nonOverlapPatch().name()
<< " size:" << patch.nonOverlapPatch().size()
<< " mask size:" << patch.neighbPatch().tgtMask().size()
<< endl;
}
srcMask_ =
min(scalar(1) - tolerance_, max(tolerance_, AMI.srcWeightsSum()));
tgtMask_ =
min(scalar(1) - tolerance_, max(tolerance_, AMI.tgtWeightsSum()));
}
}
// Adapt owner side areas. Note that in uncoupled situations (e.g.
// decomposePar) srcMask, tgtMask can be zero size.
if (srcMask_.size())
void Foam::cyclicACMIPolyPatch::scalePatchFaceAreas()
{
if (!owner() || !canResetAMI())
{
return;
}
scalePatchFaceAreas(*this);
scalePatchFaceAreas(this->neighbPatch());
}
void Foam::cyclicACMIPolyPatch::scalePatchFaceAreas
(
const cyclicACMIPolyPatch& acmipp
)
{
// Primitive patch face areas have been cleared/reset based on the raw
// points - need to reset to avoid double-accounting of face areas
DebugPout
<< "rescaling non-overlap patch areas" << endl;
const scalar maxTol = scalar(1) - tolerance_;
const scalarField& mask = acmipp.mask();
const polyPatch& nonOverlapPatch = acmipp.nonOverlapPatch();
vectorField::subField noSf = nonOverlapPatch.faceAreas();
if (mask.size() != noSf.size())
{
WarningInFunction
<< "Inconsistent sizes for patch: " << acmipp.name()
<< " - not manipulating patches" << nl
<< " - size: " << size() << nl
<< " - non-overlap patch size: " << noSf.size() << nl
<< " - mask size: " << mask.size() << nl
<< "This is OK for decomposition but should be considered fatal "
<< "at run-time" << endl;
return;
}
forAll(noSf, facei)
{
const scalar w = min(maxTol, max(tolerance_, mask[facei]));
noSf[facei] *= scalar(1) - w;
}
if (!createAMIFaces_)
{
// Note: for topological update (createAMIFaces_ = true)
// AMI coupled patch face areas are updated as part of the topological
// updates, e.g. by the calls to cyclicAMIPolyPatch's setTopology and
// initMovePoints
DebugPout
<< "scaling coupled patch areas" << endl;
// Scale the coupled patch face areas
vectorField::subField Sf = acmipp.faceAreas();
forAll(Sf, facei)
{
vectorField::subField Sf = faceAreas();
vectorField::subField noSf = nonOverlapPatch.faceAreas();
forAll(Sf, facei)
{
Sf[facei] *= srcMask_[facei];
noSf[facei] *= 1.0 - srcMask_[facei];
}
}
// Adapt slave side areas
if (tgtMask_.size())
{
const cyclicACMIPolyPatch& cp =
refCast<const cyclicACMIPolyPatch>(this->neighbPatch());
const polyPatch& pp = cp.nonOverlapPatch();
vectorField::subField Sf = cp.faceAreas();
vectorField::subField noSf = pp.faceAreas();
forAll(Sf, facei)
{
Sf[facei] *= tgtMask_[facei];
noSf[facei] *= 1.0 - tgtMask_[facei];
}
Sf[facei] *= max(tolerance_, mask[facei]);
}
// Re-normalise the weights since the effect of overlap is already
// accounted for in the area.
// accounted for in the area
auto& weights = const_cast<scalarListList&>(acmipp.weights());
auto& weightsSum = const_cast<scalarField&>(acmipp.weightsSum());
forAll(weights, i)
{
scalarListList& srcWeights = AMI.srcWeights();
scalarField& srcWeightsSum = AMI.srcWeightsSum();
forAll(srcWeights, i)
scalarList& wghts = weights[i];
if (wghts.size())
{
scalarList& wghts = srcWeights[i];
if (wghts.size())
{
scalar& sum = srcWeightsSum[i];
scalar& sum = weightsSum[i];
forAll(wghts, j)
{
wghts[j] /= sum;
}
sum = 1.0;
forAll(wghts, j)
{
wghts[j] /= sum;
}
sum = 1.0;
}
}
{
scalarListList& tgtWeights = AMI.tgtWeights();
scalarField& tgtWeightsSum = AMI.tgtWeightsSum();
forAll(tgtWeights, i)
{
scalarList& wghts = tgtWeights[i];
if (wghts.size())
{
scalar& sum = tgtWeightsSum[i];
forAll(wghts, j)
{
wghts[j] /= sum;
}
sum = 1.0;
}
}
}
// Set the updated flag
updated_ = true;
}
}
@ -244,9 +274,9 @@ void Foam::cyclicACMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
// Note: calculates transformation and triggers face centre calculation
cyclicAMIPolyPatch::initGeometry(pBufs);
// Initialise the AMI early to make sure we adapt the face areas before the
// cell centre calculation gets triggered.
resetAMI();
// On start-up there are no topological updates so scale the face areas
// - Note: resetAMI called by cyclicAMIPolyPatch::initGeometry
scalePatchFaceAreas();
}
@ -272,10 +302,10 @@ void Foam::cyclicACMIPolyPatch::initMovePoints
}
// Note: calculates transformation and triggers face centre calculation
// - Note: resetAMI called by cyclicAMIPolyPatch::initMovePoints
cyclicAMIPolyPatch::initMovePoints(pBufs, p);
// Initialise the AMI early. See initGeometry.
resetAMI();
scalePatchFaceAreas();
}
@ -289,6 +319,8 @@ void Foam::cyclicACMIPolyPatch::movePoints
{
Pout<< "cyclicACMIPolyPatch::movePoints : " << name() << endl;
}
// When topology is changing, this will scale the duplicate AMI faces
cyclicAMIPolyPatch::movePoints(pBufs, p);
}
@ -352,8 +384,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
nonOverlapPatchName_(word::null),
nonOverlapPatchID_(-1),
srcMask_(),
tgtMask_(),
updated_(false)
tgtMask_()
{
AMIRequireMatch_ = false;
@ -375,8 +406,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
nonOverlapPatchName_(dict.lookup("nonOverlapPatch")),
nonOverlapPatchID_(-1),
srcMask_(),
tgtMask_(),
updated_(false)
tgtMask_()
{
AMIRequireMatch_ = false;
@ -403,8 +433,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
nonOverlapPatchName_(pp.nonOverlapPatchName_),
nonOverlapPatchID_(-1),
srcMask_(),
tgtMask_(),
updated_(false)
tgtMask_()
{
AMIRequireMatch_ = false;
@ -428,8 +457,7 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
nonOverlapPatchName_(nonOverlapPatchName),
nonOverlapPatchID_(-1),
srcMask_(),
tgtMask_(),
updated_(false)
tgtMask_()
{
AMIRequireMatch_ = false;
@ -459,19 +487,12 @@ Foam::cyclicACMIPolyPatch::cyclicACMIPolyPatch
nonOverlapPatchName_(pp.nonOverlapPatchName_),
nonOverlapPatchID_(-1),
srcMask_(),
tgtMask_(),
updated_(false)
tgtMask_()
{
AMIRequireMatch_ = false;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cyclicACMIPolyPatch::~cyclicACMIPolyPatch()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::cyclicACMIPolyPatch& Foam::cyclicACMIPolyPatch::neighbPatch() const

54
src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -60,9 +60,6 @@ private:
// Private data
//- Fraction of face area below which face is considered disconnected
static const scalar tolerance_;
//- Name of non-overlapping patch
const word nonOverlapPatchName_;
@ -75,21 +72,40 @@ private:
//- Mask/weighting for target patch
mutable scalarField tgtMask_;
//- Flag to indicate that AMI has been updated
mutable bool updated_;
protected:
// Protected Member Functions
//- Reset the AMI interpolator
//- Report coverage statics, e.g. number of uncovered/blended/covered
//- faces
void reportCoverage
(
const word& name,
const scalarField& weightSum
) const;
//- Reset the AMI interpolator, supply patch points
virtual void resetAMI
(
const UList<point>& points,
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
AMIPatchToPatchInterpolation::imFaceAreaWeight
) const;
//- Reset the AMI interpolator, use current patch points
virtual void resetAMI
(
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
AMIPatchToPatchInterpolation::imFaceAreaWeight
) const;
//- Scale patch face areas to maintain physical area
virtual void scalePatchFaceAreas();
//- Scale patch face areas to maintain physical area
virtual void scalePatchFaceAreas(const cyclicACMIPolyPatch& acmipp);
//- Initialise the calculation of the patch geometry
virtual void initGeometry(PstreamBuffers&);
@ -111,12 +127,6 @@ protected:
//- Clear geometry
virtual void clearGeom();
//- Return the mask/weighting for the source patch
virtual const scalarField& srcMask() const;
//- Return the mask/weighting for the target patch
virtual const scalarField& tgtMask() const;
public:
@ -235,19 +245,13 @@ public:
//- Destructor
virtual ~cyclicACMIPolyPatch();
virtual ~cyclicACMIPolyPatch() = default;
// Member Functions
// Access
//- Reset the updated flag
inline void setUpdated(bool flag) const;
//- Return access to the updated flag
inline bool updated() const;
//- Return a reference to the neighbour patch
virtual const cyclicACMIPolyPatch& neighbPatch() const;
@ -263,9 +267,15 @@ public:
//- Return a reference to the non-overlapping patch
inline polyPatch& nonOverlapPatch();
//- Mask field where 1 = overlap, 0 = no-overlap
//- Mask field where 1 = overlap(coupled), 0 = no-overlap
inline const scalarField& mask() const;
//- Return the mask/weighting for the source patch
virtual const scalarField& srcMask() const;
//- Return the mask/weighting for the target patch
virtual const scalarField& tgtMask() const;
//- Overlap tolerance
inline static scalar tolerance();

18
src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatchI.H
View File

@ -27,18 +27,6 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline void Foam::cyclicACMIPolyPatch::setUpdated(const bool flag) const
{
updated_ = flag;
}
inline bool Foam::cyclicACMIPolyPatch::updated() const
{
return updated_;
}
inline const Foam::word& Foam::cyclicACMIPolyPatch::nonOverlapPatchName() const
{
return nonOverlapPatchName_;
@ -69,10 +57,8 @@ inline const Foam::scalarField& Foam::cyclicACMIPolyPatch::mask() const
{
return srcMask_;
}
else
{
return neighbPatch().tgtMask();
}
return neighbPatch().tgtMask();
}

299
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
View File

@ -27,13 +27,10 @@ License
\*---------------------------------------------------------------------------*/
#include "cyclicAMIPolyPatch.H"
#include "transformField.H"
#include "SubField.H"
#include "polyMesh.H"
#include "Time.H"
#include "addToRunTimeSelectionTable.H"
#include "faceAreaIntersect.H"
#include "ops.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -45,6 +42,7 @@ namespace Foam
addToRunTimeSelectionTable(polyPatch, cyclicAMIPolyPatch, dictionary);
}
const Foam::scalar Foam::cyclicAMIPolyPatch::tolerance_ = 1e-10;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@ -295,77 +293,114 @@ void Foam::cyclicAMIPolyPatch::resetAMI
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
) const
{
if (owner())
resetAMI(boundaryMesh().mesh().points(), AMIMethod);
}
void Foam::cyclicAMIPolyPatch::resetAMI
(
const UList<point>& points,
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod
) const
{
DebugInFunction << endl;
if (!owner())
{
AMIPtr_.clear();
return;
}
const polyPatch& nbr = neighbPatch();
pointField nbrPoints
(
neighbPatch().boundaryMesh().mesh().points(),
neighbPatch().meshPoints()
);
AMIPtr_.clear();
if (debug)
const cyclicAMIPolyPatch& nbr = neighbPatch();
pointField srcPoints(points, meshPoints());
pointField nbrPoints(points, neighbPatch().meshPoints());
if (debug)
{
const Time& t = boundaryMesh().mesh().time();
OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
}
label patchSize0 = size();
label nbrPatchSize0 = nbr.size();
if (createAMIFaces_)
{
// AMI is created based on the original patch faces (non-extended patch)
if (srcFaceIDs_.size())
{
const Time& t = boundaryMesh().mesh().time();
OFstream os(t.path()/name() + "_neighbourPatch-org.obj");
meshTools::writeOBJ(os, neighbPatch().localFaces(), nbrPoints);
patchSize0 = srcFaceIDs_.size();
}
// Transform neighbour patch to local system
transformPosition(nbrPoints);
primitivePatch nbrPatch0
(
SubList<face>
(
nbr.localFaces(),
nbr.size()
),
nbrPoints
);
if (debug)
if (tgtFaceIDs_.size())
{
const Time& t = boundaryMesh().mesh().time();
OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
OFstream osO(t.path()/name() + "_ownerPatch.obj");
meshTools::writeOBJ(osO, this->localFaces(), localPoints());
nbrPatchSize0 = tgtFaceIDs_.size();
}
}
// Construct/apply AMI interpolation to determine addressing and weights
AMIPtr_.reset
// Transform neighbour patch to local system
transformPosition(nbrPoints);
primitivePatch nbrPatch0
(
SubList<face>
(
new AMIPatchToPatchInterpolation
(
*this,
nbrPatch0,
surfPtr(),
faceAreaIntersect::tmMesh,
AMIRequireMatch_,
AMIMethod,
AMILowWeightCorrection_,
AMIReverse_
)
);
nbr.localFaces(),
nbrPatchSize0
),
nbrPoints
);
primitivePatch patch0
(
SubList<face>
(
localFaces(),
patchSize0
),
srcPoints
);
if (debug)
{
Pout<< "cyclicAMIPolyPatch : " << name()
<< " constructed AMI with " << nl
<< " " << "srcAddress:" << AMIPtr_().srcAddress().size()
<< nl
<< " " << "tgAddress :" << AMIPtr_().tgtAddress().size()
<< nl << endl;
}
if (debug)
{
const Time& t = boundaryMesh().mesh().time();
OFstream osN(t.path()/name() + "_neighbourPatch-trans.obj");
meshTools::writeOBJ(osN, nbrPatch0.localFaces(), nbrPoints);
OFstream osO(t.path()/name() + "_ownerPatch.obj");
meshTools::writeOBJ(osO, this->localFaces(), localPoints());
}
// Construct/apply AMI interpolation to determine addressing and weights
AMIPtr_.reset
(
new AMIPatchToPatchInterpolation
(
patch0, // *this,
nbrPatch0,
surfPtr(),
faceAreaIntersect::tmMesh,
AMIRequireMatch_,
AMIMethod,
AMILowWeightCorrection_,
AMIReverse_
)
);
// Set the updated flag
updated_ = true;
if (debug)
{
AMIPtr_->checkSymmetricWeights(true);
}
}
void Foam::cyclicAMIPolyPatch::calcTransforms()
{
DebugInFunction << endl;
const cyclicAMIPolyPatch& half0 = *this;
vectorField half0Areas(half0.size());
forAll(half0, facei)
@ -402,22 +437,29 @@ void Foam::cyclicAMIPolyPatch::calcTransforms()
void Foam::cyclicAMIPolyPatch::initGeometry(PstreamBuffers& pBufs)
{
// The AMI is no longer valid. Leave it up to demand-driven calculation
AMIPtr_.clear();
DebugInFunction << endl;
// Only resetting the AMI if not creating additional AMI faces
// - When creating additional faces the AMI is updated externally by the
// dynamic mesh via the setTopology function.
if (!createAMIFaces_ && canResetAMI())
{
resetAMI(AMIMethod_);
}
polyPatch::initGeometry(pBufs);
// Early calculation of transforms so e.g. cyclicACMI can use them.
// Note: also triggers primitiveMesh face centre. Note that cell
// centres should -not- be calculated
// since e.g. cyclicACMI override face areas
// since e.g. cyclicACMI overrides face areas
calcTransforms();
}
void Foam::cyclicAMIPolyPatch::calcGeometry(PstreamBuffers& pBufs)
{
// All geometry done inside initGeometry
DebugInFunction << endl;
}
@ -427,14 +469,31 @@ void Foam::cyclicAMIPolyPatch::initMovePoints
const pointField& p
)
{
// The AMI is no longer valid. Leave it up to demand-driven calculation
AMIPtr_.clear();
polyPatch::initMovePoints(pBufs, p);
DebugInFunction << endl;
// See below. Clear out any local geometry
primitivePatch::movePoints(p);
// Note: processorPolyPatch::initMovePoints calls
// processorPolyPatch::initGeometry which will trigger calculation of
// patch faceCentres() and cell volumes...
if (createAMIFaces_)
{
// Note: AMI should have been updated in setTopology
// faceAreas() and faceCentres() have been reset and will be
// recalculated on-demand using the mesh points and no longer
// correspond to the scaled areas!
restoreScaledGeometry();
// deltas need to be recalculated to use new face centres!
}
else
{
resetAMI(p, AMIMethod_);
}
// Early calculation of transforms. See above.
calcTransforms();
}
@ -446,30 +505,55 @@ void Foam::cyclicAMIPolyPatch::movePoints
const pointField& p
)
{
polyPatch::movePoints(pBufs, p);
DebugInFunction << endl;
// All transformation tensors already done in initMovePoints
// Note: not calling movePoints since this will undo our manipulations!
// polyPatch::movePoints(pBufs, p);
/*
polyPatch::movePoints
-> primitivePatch::movePoints
-> primitivePatch::clearGeom:
deleteDemandDrivenData(localPointsPtr_);
deleteDemandDrivenData(faceCentresPtr_);
deleteDemandDrivenData(faceAreasPtr_);
deleteDemandDrivenData(magFaceAreasPtr_);
deleteDemandDrivenData(faceNormalsPtr_);
deleteDemandDrivenData(pointNormalsPtr_);
*/
}
void Foam::cyclicAMIPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
{
// The AMI is no longer valid. Leave it up to demand-driven calculation
AMIPtr_.clear();
DebugInFunction << endl;
polyPatch::initUpdateMesh(pBufs);
if (createAMIFaces_ && boundaryMesh().mesh().topoChanging() && owner())
{
setAMIFaces();
}
}
void Foam::cyclicAMIPolyPatch::updateMesh(PstreamBuffers& pBufs)
{
DebugInFunction << endl;
// Note: this clears out cellCentres(), faceCentres() and faceAreas()
polyPatch::updateMesh(pBufs);
}
void Foam::cyclicAMIPolyPatch::clearGeom()
{
AMIPtr_.clear();
DebugInFunction << endl;
if (!updatingAMI_)
{
//DEBUG("*** CLEARING AMI ***");
AMIPtr_.clear();
}
polyPatch::clearGeom();
}
@ -488,6 +572,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
)
:
coupledPolyPatch(name, size, start, index, bm, patchType, transform),
updated_(false),
nbrPatchName_(word::null),
nbrPatchID_(-1),
rotationAxis_(Zero),
@ -501,7 +586,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
AMIRequireMatch_(true),
AMILowWeightCorrection_(-1.0),
surfPtr_(nullptr),
surfDict_(fileName("surface"))
surfDict_(fileName("surface")),
createAMIFaces_(false),
updatingAMI_(true),
srcFaceIDs_(),
tgtFaceIDs_(),
faceAreas0_(),
faceCentres0_()
{
// Neighbour patch might not be valid yet so no transformation
// calculation possible
@ -518,6 +609,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
)
:
coupledPolyPatch(name, dict, index, bm, patchType),
updated_(false),
nbrPatchName_(dict.getOrDefault<word>("neighbourPatch", "")),
coupleGroup_(dict),
nbrPatchID_(-1),
@ -545,7 +637,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
AMIRequireMatch_(true),
AMILowWeightCorrection_(dict.getOrDefault("lowWeightCorrection", -1.0)),
surfPtr_(nullptr),
surfDict_(dict.subOrEmptyDict("surface"))
surfDict_(dict.subOrEmptyDict("surface")),
createAMIFaces_(false),
updatingAMI_(true),
srcFaceIDs_(),
tgtFaceIDs_(),
faceAreas0_(),
faceCentres0_()
{
if (nbrPatchName_ == word::null && !coupleGroup_.valid())
{
@ -605,6 +703,20 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
// Neighbour patch might not be valid yet so no transformation
// calculation possible
// If topology change, recover the sizes of the original patches
label srcSize0 = 0;
if (dict.readIfPresent("srcSize", srcSize0))
{
srcFaceIDs_.setSize(srcSize0);
createAMIFaces_ = true;
}
label tgtSize0 = 0;
if (dict.readIfPresent("tgtSize", tgtSize0))
{
tgtFaceIDs_.setSize(tgtSize0);
createAMIFaces_ = true;
}
}
@ -615,6 +727,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
)
:
coupledPolyPatch(pp, bm),
updated_(false),
nbrPatchName_(pp.nbrPatchName_),
coupleGroup_(pp.coupleGroup_),
nbrPatchID_(-1),
@ -629,7 +742,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
AMIRequireMatch_(pp.AMIRequireMatch_),
AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
surfPtr_(nullptr),
surfDict_(pp.surfDict_)
surfDict_(pp.surfDict_),
createAMIFaces_(false),
updatingAMI_(true),
srcFaceIDs_(),
tgtFaceIDs_(),
faceAreas0_(),
faceCentres0_()
{
// Neighbour patch might not be valid yet so no transformation
// calculation possible
@ -647,6 +766,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
)
:
coupledPolyPatch(pp, bm, index, newSize, newStart),
updated_(false),
nbrPatchName_(nbrPatchName),
coupleGroup_(pp.coupleGroup_),
nbrPatchID_(-1),
@ -661,7 +781,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
AMIRequireMatch_(pp.AMIRequireMatch_),
AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
surfPtr_(nullptr),
surfDict_(pp.surfDict_)
surfDict_(pp.surfDict_),
createAMIFaces_(false),
updatingAMI_(true),
srcFaceIDs_(),
tgtFaceIDs_(),
faceAreas0_(),
faceCentres0_()
{
if (nbrPatchName_ == name())
{
@ -686,6 +812,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
)
:
coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
updated_(false),
nbrPatchName_(pp.nbrPatchName_),
coupleGroup_(pp.coupleGroup_),
nbrPatchID_(-1),
@ -700,13 +827,13 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
AMIRequireMatch_(pp.AMIRequireMatch_),
AMILowWeightCorrection_(pp.AMILowWeightCorrection_),
surfPtr_(nullptr),
surfDict_(pp.surfDict_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cyclicAMIPolyPatch::~cyclicAMIPolyPatch()
surfDict_(pp.surfDict_),
createAMIFaces_(false),
updatingAMI_(true),
srcFaceIDs_(),
tgtFaceIDs_(),
faceAreas0_(),
faceCentres0_()
{}
@ -1105,6 +1232,12 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
{
surfDict_.writeEntry(surfDict_.dictName(), os);
}
if (createAMIFaces_)
{
os.writeEntry("srcSize", srcFaceIDs_.size());
os.writeEntry("tgtSize", tgtFaceIDs_.size());
}
}

95
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H
View File

@ -74,6 +74,9 @@ protected:
// Protected data
//- Flag to indicate that AMI has been updated
mutable bool updated_;
//- Name of other half
mutable word nbrPatchName_;
@ -129,9 +132,47 @@ protected:
const dictionary surfDict_;
// Change of topology as AMI is updated
//- Flag to indicate that new AMI faces will created
// Set by the call to changeTopology
mutable bool createAMIFaces_;
mutable bool updatingAMI_;
labelListList srcFaceIDs_;
labelListList tgtFaceIDs_;
//- Temporary storage for AMI face areas
mutable vectorField faceAreas0_;
//- Temporary storage for AMI face centres
mutable vectorField faceCentres0_;
// Protected Member Functions
//- Reset the AMI interpolator
// Topology change
virtual bool removeAMIFaces(polyTopoChange& topoChange);
virtual bool addAMIFaces(polyTopoChange& topoChange);
virtual void setAMIFaces();
virtual void restoreScaledGeometry();
//- Reset the AMI interpolator, supply patch points
virtual void resetAMI
(
const UList<point>& points,
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
AMIPatchToPatchInterpolation::imFaceAreaWeight
) const;
//- Reset the AMI interpolator, use current patch points
virtual void resetAMI
(
const AMIPatchToPatchInterpolation::interpolationMethod& AMIMethod =
@ -274,15 +315,45 @@ public:
//- Destructor
virtual ~cyclicAMIPolyPatch();
virtual ~cyclicAMIPolyPatch() = default;
// Member Functions
// Access
//- Tolerance used e.g. for area calculations/limits
static const scalar tolerance_;
//- Flag to indicate whether the AMI can be reset
inline bool canResetAMI() const;
//- Reset the updated flag
inline void setUpdated(bool flag) const;
//- Return access to the updated flag
inline bool updated() const;
//- Return access to the createAMIFaces flag
inline bool createAMIFaces() const;
//- Return access to the updated flag
inline bool updatingAMI() const;
void clearAMI() const
{
AMIPtr_.clear();
}
//- Return true if this patch changes the mesh topology
// True when createAMIFaces is true
virtual bool changeTopology() const;
//- Set topology changes in the polyTopoChange object
virtual bool setTopology(polyTopoChange& topoChange);
//- Is patch 'coupled'. Note that on AMI the geometry is not
// coupled but the fields are!
//- coupled but the fields are!
virtual bool coupled() const
{
return false;
@ -306,9 +377,23 @@ public:
//- Return a reference to the AMI interpolator
const AMIPatchToPatchInterpolation& AMI() const;
//- Helper function to return the weights
inline const scalarListList& weights() const;
//- Helper function to return the weights sum
inline const scalarField& weightsSum() const;
//- Return true if applying the low weight correction
bool applyLowWeightCorrection() const;
//- Return access to the initial face areas
// Used for topology change
inline vectorField& faceAreas0() const;
//- Return access to the initial face centres
// Used for topology change
inline vectorField& faceCentres0() const;
// Transformations
@ -388,7 +473,7 @@ public:
);
//- Initialize ordering for primitivePatch. Does not
// refer to *this (except for name() and type() etc.)
//- refer to *this (except for name() and type() etc.)
virtual void initOrder
(
PstreamBuffers&,
@ -409,7 +494,7 @@ public:
) const;
//- Return face index on neighbour patch which shares point p
// following trajectory vector n
//- following trajectory vector n
label pointFace
(
const label facei,

63
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2013 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,6 +28,36 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline bool Foam::cyclicAMIPolyPatch::canResetAMI() const
{
return Pstream::parRun() || !boundaryMesh().mesh().time().processorCase();
}
inline void Foam::cyclicAMIPolyPatch::setUpdated(const bool flag) const
{
updated_ = flag;
}
inline bool Foam::cyclicAMIPolyPatch::updated() const
{
return updated_;
}
inline bool Foam::cyclicAMIPolyPatch::createAMIFaces() const
{
return createAMIFaces_;
}
inline bool Foam::cyclicAMIPolyPatch::updatingAMI() const
{
return updatingAMI_;
}
inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
{
if (nbrPatchName_.empty())
@ -39,6 +70,38 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
return nbrPatchName_;
}
inline const Foam::scalarListList& Foam::cyclicAMIPolyPatch::weights() const
{
if (owner())
{
return AMI().srcWeights();
}
return neighbPatch().AMI().tgtWeights();
}
inline const Foam::scalarField& Foam::cyclicAMIPolyPatch::weightsSum() const
{
if (owner())
{
return AMI().srcWeightsSum();
}
return neighbPatch().AMI().tgtWeightsSum();
}
inline Foam::vectorField& Foam::cyclicAMIPolyPatch::faceAreas0() const
{
return faceAreas0_;
}
inline Foam::vectorField& Foam::cyclicAMIPolyPatch::faceCentres0() const
{
return faceCentres0_;
}
inline const Foam::vector& Foam::cyclicAMIPolyPatch::rotationAxis() const
{

676
src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C Normal file
View File

@ -0,0 +1,676 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "cyclicAMIPolyPatch.H"
#include "SubField.H"
#include "vectorList.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
void Foam::cyclicAMIPolyPatch::restoreScaledGeometry()
{
DebugInFunction << endl;
// Note: only used for topology update (createAMIFaces_ flag)
if (!createAMIFaces_)
{
FatalErrorInFunction
<< "Attempted to perform topology update when createAMIFaces_ "
<< "flag is set to false"
<< abort(FatalError);
}
if (boundaryMesh().mesh().hasCellVolumes())
{
WarningInFunction
<< "Mesh already has volumes set!"
<< endl;
}
vectorField::subField faceAreas = this->faceAreas();
vectorField::subField faceCentres = this->faceCentres();
if (debug)
{
Info<< "Patch:" << name() << " before: sum(mag(faceAreas)):"
<< gSum(mag(faceAreas)) << endl;
Info<< "Patch:" << name() << " before: sum(mag(faceAreas0)):"
<< gSum(mag(faceAreas0_)) << endl;
}
faceAreas = faceAreas0_;
Info<< "WARNING: ==============================" << endl;
Info<< "WARNING: Not updating face cell centres" << endl;
Info<< "WARNING: ==============================" << endl;
//faceCentres = faceCentres0_;
faceAreas0_.clear();
faceCentres0_.clear();
if (debug)
{
Info<< "Patch:" << name() << " after: sum(mag(faceAreas)):"
<< gSum(mag(faceAreas)) << endl;
Info<< "Patch:" << name() << " after: sum(mag(faceAreas0)):"
<< gSum(mag(faceAreas0_)) << endl;
}
}
bool Foam::cyclicAMIPolyPatch::removeAMIFaces(polyTopoChange& topoChange)
{
DebugInFunction << endl;
// Note: only used for topology update (createAMIFaces_ flag)
if (!createAMIFaces_)
{
FatalErrorInFunction
<< "Attempted to perform topology update when createAMIFaces_ "
<< "flag is set to false"
<< abort(FatalError);
}
if (!owner())
{
return false;
}
bool changeRequired = false;
// Remove any faces that we inserted to create the 1-to-1 match...
const cyclicAMIPolyPatch& nbr = neighbPatch();
const label newSrcFaceStart = srcFaceIDs_.size();
if (newSrcFaceStart != 0)
{
for (label facei = newSrcFaceStart; facei < size(); ++facei)
{
changeRequired = true;
label meshFacei = start() + facei;
topoChange.removeFace(meshFacei, -1);
}
}
const label newTgtFaceStart = tgtFaceIDs_.size();
if (newTgtFaceStart != 0)
{
for (label facei = newTgtFaceStart; facei < nbr.size(); ++facei)
{
changeRequired = true;
label meshFacei = nbr.start() + facei;
topoChange.removeFace(meshFacei, -1);
}
}
srcFaceIDs_.clear();
tgtFaceIDs_.clear();
return changeRequired;
}
bool Foam::cyclicAMIPolyPatch::addAMIFaces(polyTopoChange& topoChange)
{
DebugInFunction << endl;
// Note: only used for topology update (createAMIFaces_ flag = true)
if (!createAMIFaces_)
{
FatalErrorInFunction
<< "Attempted to perform topology update when createAMIFaces_ "
<< "flag is set to false"
<< abort(FatalError);
}
bool changedFaces = false;
const cyclicAMIPolyPatch& nbr = neighbPatch();
polyMesh& mesh = const_cast<polyMesh&>(boundaryMesh().mesh());
const faceZoneMesh& faceZones = mesh.faceZones();
// First face address and weight are used to manipulate the
// original face - all other addresses and weights are used to
// create additional faces
const labelListList& srcToTgtAddr = AMI().srcAddress();
const labelListList& tgtToSrcAddr = AMI().tgtAddress();
const label nSrcFace = srcToTgtAddr.size();
const label nTgtFace = tgtToSrcAddr.size();
srcFaceIDs_.setSize(nSrcFace);
tgtFaceIDs_.setSize(nTgtFace);
label nNewSrcFaces = 0;
forAll(srcToTgtAddr, srcFacei)
{
const labelList& tgtAddr = srcToTgtAddr[srcFacei];
// No tgt faces linked to srcFacei (ACMI)
if (tgtAddr.empty()) continue;
srcFaceIDs_[srcFacei].setSize(tgtAddr.size());
srcFaceIDs_[srcFacei][0] = srcFacei;
const label meshFacei = start() + srcFacei;
for (label addri = 1; addri < tgtAddr.size(); ++addri)
{
changedFaces = true;
// Note: new faces reuse originating face points
// - but areas are scaled by the weights (later)
// New source face for each target face address
srcFaceIDs_[srcFacei][addri] = nNewSrcFaces + nSrcFace;
++nNewSrcFaces;
(void)topoChange.addFace
(
mesh.faces()[meshFacei], // modified face
mesh.faceOwner()[meshFacei], // owner
-1, // neighbour
-1, // master point
-1, // master edge
meshFacei, // master face
false, // face flip
index(), // patch for face
faceZones.whichZone(meshFacei), // zone for original face
false // face flip in zone
);
}
}
label nNewTgtFaces = 0;
forAll(tgtToSrcAddr, tgtFacei)
{
const labelList& srcAddr = tgtToSrcAddr[tgtFacei];
// No src faces linked to tgtFacei (ACMI)
if (srcAddr.empty()) continue;
tgtFaceIDs_[tgtFacei].setSize(srcAddr.size());
tgtFaceIDs_[tgtFacei][0] = tgtFacei;
const label meshFacei = nbr.start() + tgtFacei;
for (label addri = 1; addri < srcAddr.size(); ++addri)
{
changedFaces = true;
// Note: new faces reuse originating face points
// - but areas are scaled by the weights (later)
// New target face for each source face address
tgtFaceIDs_[tgtFacei][addri] = nNewTgtFaces + nTgtFace;
++nNewTgtFaces;
(void)topoChange.addFace
(
mesh.faces()[meshFacei], // modified face
mesh.faceOwner()[meshFacei], // owner
-1, // neighbour
-1, // master point
-1, // master edge
meshFacei, // master face
false, // face flip
nbr.index(), // patch for face
faceZones.whichZone(meshFacei), // zone for original face
false // face flip in zone
);
}
}
Info<< "AMI: Patch " << name() << " additional faces: "
<< returnReduce(nNewSrcFaces, sumOp<label>()) << nl
<< "AMI: Patch " << nbr.name() << " additional faces: "
<< returnReduce(nNewTgtFaces, sumOp<label>())
<< endl;
if (debug)
{
Pout<< "New faces - " << name() << ": " << nNewSrcFaces
<< " " << nbr.name() << ": " << nNewTgtFaces << endl;
}
return returnReduce(changedFaces, orOp<bool>());
}
void Foam::cyclicAMIPolyPatch::setAMIFaces()
{
// Create new mesh faces so that there is a 1-to-1 correspondence
// between faces on each side of the AMI
// Note: only used for topology update (createAMIFaces_ flag)
if (!createAMIFaces_)
{
FatalErrorInFunction
<< "Attempted to perform topology update when createAMIFaces_ "
<< "flag is set to false"
<< abort(FatalError);
}
DebugInFunction << endl;
if (!owner())
{
return;
}
const cyclicAMIPolyPatch& nbr = neighbPatch();
vectorField& nbrFaceAreas0 = nbr.faceAreas0();
vectorField& nbrFaceCentres0 = nbr.faceCentres0();
// Scale the new face areas and set the centroids
// Note:
// - storing local copies so that they can be re-applied after the call to
// movePoints that will reset any changes to the areas and centroids
//
// - For AMI, src and tgt patches should be the same
// - For ACMI they are likely to be different!
faceAreas0_ = faceAreas();
faceCentres0_ = faceCentres();
nbrFaceAreas0 = nbr.faceAreas();
nbrFaceCentres0 = nbr.faceCentres();
// Original AMI info (based on the mesh state when the AMI was evaluated)
const labelListList& srcToTgtAddr0 = AMIPtr_->srcAddress();
const labelListList& tgtToSrcAddr0 = AMIPtr_->tgtAddress();
const pointListList& srcCtr0 = AMIPtr_->srcCentroids();
const scalarListList& srcToTgtWght0 = AMIPtr_->srcWeights();
// New addressing on new mesh (extended by polyTopoChange)
labelListList srcToTgtAddr1(size(), labelList());
labelListList tgtToSrcAddr1(nbr.size(), labelList());
// Need to calc new parallel maps (mesh has changed since AMI was computed)
autoPtr<mapDistribute> srcToTgtMap1;
autoPtr<mapDistribute> tgtToSrcMap1;
if (AMIPtr_->singlePatchProc() == -1)
{
// Parallel running
// Global index based on old patch sizes (when AMI was computed)
globalIndex globalSrcFaces0(srcToTgtAddr0.size());
globalIndex globalTgtFaces0(tgtToSrcAddr0.size());
// Global index based on new patch sizes
globalIndex globalSrcFaces1(size());
globalIndex globalTgtFaces1(nbr.size());
// Gather source side info
// =======================
// Note: using new global index for addressing, and distributed using
// the old AMI map
labelListList newTgtGlobalFaces(tgtFaceIDs_);
forAll(newTgtGlobalFaces, tgtFacei)
{
globalTgtFaces1.inplaceToGlobal(newTgtGlobalFaces[tgtFacei]);
}
AMIPtr_->tgtMap().distribute(newTgtGlobalFaces);
// Now have new tgt face indices for each src face
labelList globalSrcFaceIDs(identity(srcToTgtAddr0.size()));
globalSrcFaces0.inplaceToGlobal(globalSrcFaceIDs);
AMIPtr_->srcMap().distribute(globalSrcFaceIDs);
// globalSrcFaceIDs now has remote data for each srcFacei0 known to the
// tgt patch
List<List<point>> globalSrcCtrs0(srcCtr0);
AMIPtr_->srcMap().distribute(globalSrcCtrs0);
labelList globalTgtFaceIDs(identity(tgtToSrcAddr0.size()));
globalTgtFaces0.inplaceToGlobal(globalTgtFaceIDs);
AMIPtr_->tgtMap().distribute(globalTgtFaceIDs);
// globalTgtFaceIDs now has remote data for each tgtFacei0 known to the
// src patch
// For debug - send tgt face centres and compare against mapped src
// face centres
//List<List<point>> globalTgtCtrs0(tgtCtr0);
//AMIPtr_->tgtMap().distribute(globalTgtCtrs0);
labelListList globalTgtToSrcAddr(tgtToSrcAddr0);
forAll(tgtToSrcAddr0, tgtFacei0)
{
forAll(tgtToSrcAddr0[tgtFacei0], addri)
{
const label globalSrcFacei =
globalSrcFaceIDs[tgtToSrcAddr0[tgtFacei0][addri]];
globalTgtToSrcAddr[tgtFacei0][addri] = globalSrcFacei;
}
}
AMIPtr_->tgtMap().distribute(globalTgtToSrcAddr);
labelListList globalSrcToTgtAddr(srcToTgtAddr0);
forAll(srcToTgtAddr0, srcFacei0)
{
forAll(srcToTgtAddr0[srcFacei0], addri)
{
const label globalTgtFacei =
globalTgtFaceIDs[srcToTgtAddr0[srcFacei0][addri]];
globalSrcToTgtAddr[srcFacei0][addri] = globalTgtFacei;
}
}
AMIPtr_->srcMap().distribute(globalSrcToTgtAddr);
label nError = 0;
forAll(srcToTgtAddr0, srcFacei0)
{
const labelList& newSrcFaces = srcFaceIDs_[srcFacei0];
forAll(newSrcFaces, i)
{
const label srcFacei1 = newSrcFaces[i];
// What index did srcFacei0 appear in tgtToSrc0 list?
// - if first index, all ok
// - else tgt face has been moved to according to tgtFaceIDs_
const label tgtFacei0 = srcToTgtAddr0[srcFacei0][i];
const label addri =
globalTgtToSrcAddr[tgtFacei0].find
(
globalSrcFaceIDs[srcFacei0]
);
if (addri == -1)
{
++nError;
continue;
if (debug)
{
Pout<< "Unable to find global source face "
<< globalSrcFaceIDs[srcFacei0]
<< " in globalTgtToSrcAddr[" << tgtFacei0 << "]: "
<< globalTgtToSrcAddr[tgtFacei0]
<< endl;
}
}
const label tgtFacei1 = newTgtGlobalFaces[tgtFacei0][addri];
// Sanity check to see that we've picked the correct face
// point tgtCtr0(globalTgtCtrs0[tgtFacei0][addri]);
// Pout<< "srcCtr:" << srcCtr0[srcFacei0][i]
// << " tgtCtr:" << tgtCtr0 << endl;
srcToTgtAddr1[srcFacei1] = labelList(1, tgtFacei1);
faceAreas0_[srcFacei1] *= srcToTgtWght0[srcFacei0][i];
faceCentres0_[srcFacei1] = srcCtr0[srcFacei0][i];
}
}
if (nError)
{
FatalErrorInFunction
<< "Unable to find " << nError << " global source faces"
<< abort(FatalError);
}
// Gather Target side info
// =======================
labelListList newSrcGlobalFaces(srcFaceIDs_);
forAll(newSrcGlobalFaces, srcFacei)
{
globalSrcFaces1.inplaceToGlobal(newSrcGlobalFaces[srcFacei]);
}
AMIPtr_->srcMap().distribute(newSrcGlobalFaces);
// Now have new src face indices for each tgt face
forAll(tgtToSrcAddr0, tgtFacei0)
{
const labelList& newTgtFaces = tgtFaceIDs_[tgtFacei0];
forAll(newTgtFaces, i)
{
const label srcFacei0 = tgtToSrcAddr0[tgtFacei0][i];
const label addri =
globalSrcToTgtAddr[srcFacei0].find
(
globalTgtFaceIDs[tgtFacei0]
);
if (addri == -1)
{
++nError;
continue;
if (debug)
{
Pout<< "Unable to find global target face "
<< globalTgtFaceIDs[tgtFacei0]
<< " in globalSrcToTgtAddr[" << srcFacei0 << "]: "
<< globalSrcToTgtAddr[srcFacei0]
<< endl;
}
}
const label srcFacei1 = newSrcGlobalFaces[srcFacei0][addri];
// Sanity check to see that we've picked the correct face
point srcCtr0(globalSrcCtrs0[srcFacei0][addri]);
reverseTransformPosition(srcCtr0, srcFacei0);
const label tgtFacei1 = newTgtFaces[i];
tgtToSrcAddr1[tgtFacei1] = labelList(1, srcFacei1);
nbrFaceCentres0[tgtFacei1] = srcCtr0;
}
}
if (nError)
{
FatalErrorInFunction
<< "Unable to find " << nError << " global target faces"
<< abort(FatalError);
}
// Update the maps
{
List<Map<label>> cMap;
srcToTgtMap1.reset
(
new mapDistribute(globalSrcFaces1, tgtToSrcAddr1, cMap)
);
}
{
List<Map<label>> cMap;
tgtToSrcMap1.reset
(
new mapDistribute(globalTgtFaces1, srcToTgtAddr1, cMap)
);
}
// Reset tgt patch areas using the new map
vectorList newSrcGlobalFaceAreas(faceAreas0_);
srcToTgtMap1->distribute(newSrcGlobalFaceAreas);
forAll(nbrFaceAreas0, tgtFacei)
{
if (!tgtToSrcAddr1[tgtFacei].empty())
{
const label srcFacei = tgtToSrcAddr1[tgtFacei][0];
nbrFaceAreas0[tgtFacei] = -newSrcGlobalFaceAreas[srcFacei];
}
}
}
else
{
label nError = 0;
forAll(srcToTgtAddr0, srcFacei0)
{
const labelList& srcFaceTgtAddr0 = srcToTgtAddr0[srcFacei0];
const scalarList& srcFaceTgtWght0 = srcToTgtWght0[srcFacei0];
const pointList& srcFaceTgtCtr0 = srcCtr0[srcFacei0];
forAll(srcFaceTgtAddr0, addri)
{
const label srcFacei1 = srcFaceIDs_[srcFacei0][addri];
// Find which slot srcFacei0 appears in tgt->src addressing
const label tgtFacei0 = srcFaceTgtAddr0[addri];
const label tgtAddri0 =
tgtToSrcAddr0[tgtFacei0].find(srcFacei0);
if (tgtAddri0 == -1)
{
++nError;
continue;
if (debug)
{
Pout<< "Unable to find source face " << srcFacei0
<< " in tgtToSrcAddr0[" << tgtFacei0 << "]: "
<< tgtToSrcAddr0[tgtFacei0]
<< endl;
}
}
const label tgtFacei1 = tgtFaceIDs_[tgtFacei0][tgtAddri0];
faceAreas0_[srcFacei1] *= srcFaceTgtWght0[addri];
nbrFaceAreas0[tgtFacei1] = -faceAreas0_[srcFacei1];
point pt(srcFaceTgtCtr0[addri]);
faceCentres0_[srcFacei1] = pt;
reverseTransformPosition(pt, srcFacei0);
nbrFaceCentres0[tgtFacei1] = pt;
// SANITY CHECK
// Info<< "srcPt:" << srcFaceCentres[srcFacei1]
// << " tgtPt:" << tgtFaceCentres[tgtFacei1] << endl;
srcToTgtAddr1[srcFacei1] = labelList(1, tgtFacei1);
tgtToSrcAddr1[tgtFacei1] = labelList(1, srcFacei1);
}
}
if (nError)
{
FatalErrorInFunction
<< "Unable to find " << nError
<< " source faces in tgtToSrcAddr0"
<< abort(FatalError);
}
}
scalarListList newSrcToTgtWeights(srcToTgtAddr1.size());
forAll(srcToTgtAddr1, facei)
{
if (srcToTgtAddr1[facei].size())
{
newSrcToTgtWeights[facei] = scalarList(1, scalar(1));
}
else
{
// No connection - effect of face removed by setting area to a
// a small value
faceAreas0_[facei] *= tolerance_;
}
}
scalarListList newTgtToSrcWeights(tgtToSrcAddr1.size());
forAll(tgtToSrcAddr1, facei)
{
if (tgtToSrcAddr1[facei].size())
{
newTgtToSrcWeights[facei] = scalarList(1, scalar(1));
}
else
{
// No connection - effect of face removed by setting area to a
// a small value
nbrFaceAreas0[facei] *= tolerance_;
}
}
// Update the AMI addressing and weights to reflect the new 1-to-1
// correspondence
AMIPtr_->update
(
std::move(srcToTgtMap1),
std::move(tgtToSrcMap1),
std::move(srcToTgtAddr1),
std::move(newSrcToTgtWeights),
std::move(tgtToSrcAddr1),
std::move(newTgtToSrcWeights)
);
AMIPtr_->setAreas(mag(faceAreas0_), mag(nbrFaceAreas0));
if (debug)
{
Pout<< "cyclicAMIPolyPatch : " << name()
<< " constructed AMI with " << nl
<< " " << "srcAddress:" << AMIPtr_().srcAddress().size()
<< nl
<< " " << "tgAddress :" << AMIPtr_().tgtAddress().size()
<< nl << endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::cyclicAMIPolyPatch::changeTopology() const
{
DebugInFunction << endl;
createAMIFaces_ = true;
return true;
}
bool Foam::cyclicAMIPolyPatch::setTopology(polyTopoChange& topoChange)
{
DebugInFunction << endl;
if (createAMIFaces_ && owner())
{
// Calculate the AMI using the new points
// Note: mesh still has old points
resetAMI(topoChange.points(), AMIMethod_);
removeAMIFaces(topoChange);
addAMIFaces(topoChange);
return true;
}
return false;
}
// ************************************************************************* //

1
src/meshTools/Make/files
View File

@ -271,6 +271,7 @@ AMICycPatches=$(AMI)/patches/cyclicAMI
$(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterface.C
$(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterfaceField.C
$(AMICycPatches)/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
$(AMICycPatches)/cyclicAMIPolyPatch/cyclicAMIPolyPatchTopologyChange.C
$(AMICycPatches)/cyclicAMIPointPatch/cyclicAMIPointPatch.C
$(AMICycPatches)/cyclicAMIPointPatchField/cyclicAMIPointPatchFields.C

3
src/meshTools/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
LIB_LIBS = \
-lfileFormats \

27
src/meshTools/polyTopoChange/polyTopoChange.H
View File

@ -94,7 +94,7 @@ class IOobject;
template<class T, class Container> class CompactListList;
/*---------------------------------------------------------------------------*\
Class polyTopoChange Declaration
Class polyTopoChange Declaration
\*---------------------------------------------------------------------------*/
class polyTopoChange
@ -121,7 +121,7 @@ class polyTopoChange
DynamicList<label> pointMap_;
//- For all original and added points contains new point label.
// (used to map return value of addPoint to new mesh point)
//- (used to map return value of addPoint to new mesh point)
DynamicList<label> reversePointMap_;
//- Zone of point
@ -324,10 +324,10 @@ class polyTopoChange
);
//- Remove all unused/removed points/faces/cells and update
// face ordering (always), cell ordering (bandcompression,
// orderCells=true),
// point ordering (sorted into internal and boundary points,
// orderPoints=true)
//- face ordering (always), cell ordering (bandcompression,
//- orderCells=true),
//- point ordering (sorted into internal and boundary points,
//- orderPoints=true)
void compact
(
const bool orderCells,
@ -433,7 +433,7 @@ public:
// Constructors
//- Construct without mesh. Either specify nPatches or use
// setNumPatches before trying to make a mesh (makeMesh, changeMesh)
//- setNumPatches before trying to make a mesh (makeMesh, changeMesh)
polyTopoChange(const label nPatches, const bool strict = true);
//- Construct from mesh. Adds all points/face/cells from mesh.
@ -471,15 +471,15 @@ public:
}
//- Is point removed?
// Considered removed if point is GREAT.
//- Considered removed if point is GREAT.
inline bool pointRemoved(const label pointi) const;
//- Is face removed?
// Considered removed if face is empty
//- Considered removed if face is empty
inline bool faceRemoved(const label facei) const;
//- Is cell removed?
// Considered removed if the cellMap is -2
//- Considered removed if the cellMap is -2
inline bool cellRemoved(const label celli) const;
@ -489,7 +489,7 @@ public:
void clear();
//- Add all points/faces/cells of mesh. Additional offset for patch
// or zone ids.
//- or zone ids.
void addMesh
(
const polyMesh& mesh,
@ -500,7 +500,7 @@ public:
);
//- Explicitly pre-size the dynamic storage for expected mesh
// size for if construct-without-mesh
//- size for if construct-without-mesh
void setCapacity
(
const label nPoints,
@ -590,7 +590,7 @@ public:
void removeCell(const label celli, const label mergeCelli);
//- Explicitly set the number of patches if construct-without-mesh
// used.
//- used.
inline void setNumPatches(const label nPatches);
// Other
@ -628,7 +628,6 @@ public:
const bool orderCells = false,
const bool orderPoints = false
);
};