Partial overlap GGI conservative interpolation feature. Author: Vuko Vukcevic. Merge: Hrvoje Jasak.

This commit is contained in:
Hrvoje Jasak 2018-03-20 13:47:44 +00:00
commit 6a41a3a022
34 changed files with 2572 additions and 459 deletions

View file

@ -14,5 +14,8 @@ multiMixerFvMesh/multiMixerFvMesh.C
dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
dynamicPolyRefinementFvMesh/refinementSelection/refinementSelection/refinementSelection.C
dynamicPolyRefinementFvMesh/refinementSelection/fieldBoundsRefinement/fieldBoundsRefinement.C
dynamicPolyRefinementFvMesh/refinementSelection/minCellVolumeRefinement/minCellVolumeRefinement.C
dynamicPolyRefinementFvMesh/refinementSelection/minPatchDistanceRefinement/minPatchDistanceRefinement.C
dynamicPolyRefinementFvMesh/refinementSelection/compositeRefinementSelection/compositeRefinementSelection.C
LIB = $(FOAM_LIBBIN)/libtopoChangerFvMesh

View file

@ -0,0 +1,209 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "compositeRefinementSelection.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(compositeRefinementSelection, 0);
addToRunTimeSelectionTable
(
refinementSelection,
compositeRefinementSelection,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::compositeRefinementSelection::compositeRefinementSelection
(
const fvMesh& mesh,
const dictionary& dict
)
:
refinementSelection(mesh, dict),
baseRefinementSelections_()
{
// Read basic refinement selections
PtrList<entry> baseRefSelectionEntries
(
coeffDict().lookup("baseRefinementSelections")
);
baseRefinementSelections_.setSize(baseRefSelectionEntries.size());
forAll (baseRefinementSelections_, brsI)
{
baseRefinementSelections_.set
(
brsI,
refinementSelection::New
(
mesh,
baseRefSelectionEntries[brsI].dict()
)
);
}
}
// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
Foam::compositeRefinementSelection::~compositeRefinementSelection()
{}
// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
Foam::Xfer<Foam::labelList>
Foam::compositeRefinementSelection::refinementCellCandidates() const
{
// Final refinement cell candidates are defined as the intersection of all
// sets (obtained with different refinement selection algorithms)
// In order to define the intersection in a straightforward and efficient
// way, we will create a labelField for all cells counting the number of
// times this cell has been selected for refinement. If all selection
// algorithms have selected the cell, this cell is marked as a final
// refinement candidate.
// Create field counting the number of selections
labelField nSelections(mesh().nCells(), 0);
// Loop through all base refinement selections
forAll (baseRefinementSelections_, brsI)
{
// Get refinement candidates from this base selection algorithm. Note:
// list is transferred
const labelList curRefCandidates
(
baseRefinementSelections_[brsI].refinementCellCandidates()
);
// Increment the number of selections for selected cells
forAll (curRefCandidates, i)
{
++nSelections[curRefCandidates[i]];
}
}
// Create storage for collection of final cell candidates. Assume that
// one fifth of the cells will be marked to prevent excessive resizing
dynamicLabelList refinementCandidates(mesh().nCells()/5);
// Get number of active selection algorithms
const label nBaseSelections = baseRefinementSelections_.size();
// Loop through all cells and collect final refinement candidates
forAll (nSelections, cellI)
{
if (nSelections[cellI] == nBaseSelections)
{
// Cell has been marked by all selection algorithms, append it
refinementCandidates.append(cellI);
}
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(refinementCandidates.size(), sumOp<label>())
<< " cells as refinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return refinementCandidates.xfer();
}
Foam::Xfer<Foam::labelList>
Foam::compositeRefinementSelection::unrefinementPointCandidates() const
{
// Final unrefinement point candidates are defined as the intersection of
// all sets (obtained with different refinement selection algorithms) In
// order to define the intersection in a straightforward and efficient way,
// we will create a labelField for all points counting the number of times
// this point has been selected for unrefinement. If all selection
// algorithms have selected the point, this point is marked as a final
// unrefinement split point candidate.
// Create a field counting the number of selections
labelField nSelections(mesh().nPoints(), 0);
// Loop through all base refinement selections
forAll (baseRefinementSelections_, brsI)
{
// Get unrefinement candidates from this base selection algorithm. Note:
// list is transferred
const labelList curUnrefCandidates
(
baseRefinementSelections_[brsI].unrefinementPointCandidates()
);
// Increment the number of selections for selected points
forAll (curUnrefCandidates, i)
{
++nSelections[curUnrefCandidates[i]];
}
}
// Create storage for collection of final point candidates. Assume that one
// tenth of the points will be marked to prevent excessive resizing
dynamicLabelList unrefinementCandidates(mesh().nPoints()/10);
// Get number of active selection algorithms
const label nBaseSelections = baseRefinementSelections_.size();
// Loop through all points and collect final unrefinement candidates
forAll (nSelections, pointI)
{
if (nSelections[pointI] == nBaseSelections)
{
// Point has been marked by all selection algorithms, append it
unrefinementCandidates.append(pointI);
}
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(unrefinementCandidates.size(), sumOp<label>())
<< " points as unrefinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return unrefinementCandidates.xfer();
}
// ************************************************************************* //

View file

@ -0,0 +1,113 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::compositeRefinementSelection
Description
Selection of refinement cells based on an arbitrary number of combined
"basic" selection algorithms.
SourceFiles
compositeRefinementSelection.C
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#ifndef compositeRefinementSelection_H
#define compositeRefinementSelection_H
#include "refinementSelection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class compositeRefinementSelection Declaration
\*---------------------------------------------------------------------------*/
class compositeRefinementSelection
:
public refinementSelection
{
// Private data
//- List of basic refinement selection algorithms
PtrList<refinementSelection> baseRefinementSelections_;
// Private Member Functions
//- Disallow default bitwise copy construct
compositeRefinementSelection(const compositeRefinementSelection&);
//- Disallow default bitwise assignment
void operator=(const compositeRefinementSelection&);
public:
//- Runtime type information
TypeName("compositeRefinementSelection");
// Constructors
//- Construct from components
compositeRefinementSelection
(
const fvMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~compositeRefinementSelection();
// Member Functions
// Selection of refinement/unrefinement candidates
//- Return transferable list of cells to refine
virtual Xfer<labelList> refinementCellCandidates() const;
//- Return transferable list of split points to unrefine
virtual Xfer<labelList> unrefinementPointCandidates() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,128 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "minCellVolumeRefinement.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(minCellVolumeRefinement, 0);
addToRunTimeSelectionTable
(
refinementSelection,
minCellVolumeRefinement,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::minCellVolumeRefinement::minCellVolumeRefinement
(
const fvMesh& mesh,
const dictionary& dict
)
:
refinementSelection(mesh, dict),
minCellV_(readScalar(coeffDict().lookup("minCellVolume")))
{}
// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
Foam::minCellVolumeRefinement::~minCellVolumeRefinement()
{}
// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
Foam::Xfer<Foam::labelList>
Foam::minCellVolumeRefinement::refinementCellCandidates() const
{
// Get mesh volumes
const scalarField& cellV = mesh().V().field();
// Create storage for collection of cells. Assume that almost all of the
// cells will be marked to prevent excessive resizing.
dynamicLabelList refinementCandidates(mesh().nCells());
// Loop through cells and collect refinement candidates
forAll (cellV, cellI)
{
if (cellV[cellI] > minCellV_)
{
// Cell is larger than the specified minimum, append cell for
// potential refinement
refinementCandidates.append(cellI);
}
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(refinementCandidates.size(), sumOp<label>())
<< " cells as refinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return refinementCandidates.xfer();
}
Foam::Xfer<Foam::labelList>
Foam::minCellVolumeRefinement::unrefinementPointCandidates() const
{
// Mark all points as unrefinement candidates since only split points may be
// considered for actual unrefinement and since this refinement criterion
// will be usually used in combination with others. VV, 15/Mar/2018.
// All points are unrefinement candidates
labelList unrefinementCandidates(mesh().nPoints());
forAll (unrefinementCandidates, pointI)
{
unrefinementCandidates[pointI] = pointI;
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(unrefinementCandidates.size(), sumOp<label>())
<< " points as unrefinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return unrefinementCandidates.xfer();
}
// ************************************************************************* //

View file

@ -0,0 +1,109 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::minCellVolumeRefinement
Description
Selection of refinement cells based on a minimum cell volume. All cells
with volume larger than specified minimum get selected.
SourceFiles
minCellVolumeRefinement.C
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#ifndef minCellVolumeRefinement_H
#define minCellVolumeRefinement_H
#include "refinementSelection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class minCellVolumeRefinement Declaration
\*---------------------------------------------------------------------------*/
class minCellVolumeRefinement
:
public refinementSelection
{
// Private data
//- Minimum cell volume
scalar minCellV_;
// Private Member Functions
//- Disallow default bitwise copy construct
minCellVolumeRefinement(const minCellVolumeRefinement&);
//- Disallow default bitwise assignment
void operator=(const minCellVolumeRefinement&);
public:
//- Runtime type information
TypeName("minCellVolumeRefinement");
// Constructors
//- Construct from components
minCellVolumeRefinement(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~minCellVolumeRefinement();
// Member Functions
// Selection of refinement/unrefinement candidates
//- Return transferable list of cells to refine
virtual Xfer<labelList> refinementCellCandidates() const;
//- Return transferable list of split points to unrefine
virtual Xfer<labelList> unrefinementPointCandidates() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,161 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "minPatchDistanceRefinement.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
#include "patchWave.H"
#include "polyPatchID.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(minPatchDistanceRefinement, 0);
addToRunTimeSelectionTable
(
refinementSelection,
minPatchDistanceRefinement,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::minPatchDistanceRefinement::minPatchDistanceRefinement
(
const fvMesh& mesh,
const dictionary& dict
)
:
refinementSelection(mesh, dict),
minPatchDistance_(readScalar(coeffDict().lookup("minPatchDistance")))
{}
// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
Foam::minPatchDistanceRefinement::~minPatchDistanceRefinement()
{}
// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
Foam::Xfer<Foam::labelList>
Foam::minPatchDistanceRefinement::refinementCellCandidates() const
{
// Note: recalculate distance every time due to probable topo changes
// Read patch names from dictionary
const wordList patchNames(coeffDict().lookup("distancePatches"));
// Collect patchIDs into a hash set
labelHashSet patchIDs(patchNames.size());
forAll (patchNames, patchI)
{
// Get polyPatchID
const polyPatchID pID(patchNames[patchI], mesh().boundaryMesh());
if (pID.active())
{
patchIDs.insert(pID.index());
}
else
{
FatalIOErrorIn
(
"Xfer<labelList> minPatchDistanceRefinement::"
"refinenementPatchCandidates() const",
coeffDict()
) << "Cannot find patch " << patchNames[patchI]
<< " in the mesh." << nl
<< "Available patches are: " << mesh().boundaryMesh().names()
<< abort(FatalIOError);
}
}
// Calculate distance from specified patches and do not correct for accurate
// near wall distance (false argument)
patchWave waveDistance(mesh(), patchIDs, false);
const scalarField& patchDistance = waveDistance.distance();
// Create storage for collection of cells. Assume that almost all of the
// cells will be marked to prevent excessive resizing.
dynamicLabelList refinementCandidates(mesh().nCells());
// Loop through cells and collect refinement candidates
forAll (patchDistance, cellI)
{
if (patchDistance[cellI] > minPatchDistance_)
{
// Cell is far from the patches, append it for potential refinement
refinementCandidates.append(cellI);
}
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(refinementCandidates.size(), sumOp<label>())
<< " cells as refinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return refinementCandidates.xfer();
}
Foam::Xfer<Foam::labelList>
Foam::minPatchDistanceRefinement::unrefinementPointCandidates() const
{
// Mark all points as unrefinement candidates since only split points may be
// considered for actual unrefinement and since this refinement criterion
// will be usually used in combination with others. VV, 15/Mar/2018.
// All points are unrefinement candidates
labelList unrefinementCandidates(mesh().nPoints());
forAll (unrefinementCandidates, pointI)
{
unrefinementCandidates[pointI] = pointI;
}
// Print out some information
Info<< "Selection algorithm " << type() << " selected "
<< returnReduce(unrefinementCandidates.size(), sumOp<label>())
<< " points as unrefinement candidates."
<< endl;
// Return the list in the Xfer container to prevent copying
return unrefinementCandidates.xfer();
}
// ************************************************************************* //

View file

@ -0,0 +1,109 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::minPatchDistanceRefinement
Description
Selection of refinement cells based on a minimum distance from a set of
patches. All cells farther away from the minimum distance get selected.
SourceFiles
minPatchDistanceRefinement.C
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#ifndef minPatchDistanceRefinement_H
#define minPatchDistanceRefinement_H
#include "refinementSelection.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class minPatchDistanceRefinement Declaration
\*---------------------------------------------------------------------------*/
class minPatchDistanceRefinement
:
public refinementSelection
{
// Private data
//- Minimum distance to patch
scalar minPatchDistance_;
// Private Member Functions
//- Disallow default bitwise copy construct
minPatchDistanceRefinement(const minPatchDistanceRefinement&);
//- Disallow default bitwise assignment
void operator=(const minPatchDistanceRefinement&);
public:
//- Runtime type information
TypeName("minPatchDistanceRefinement");
// Constructors
//- Construct from components
minPatchDistanceRefinement(const fvMesh& mesh, const dictionary& dict);
//- Destructor
virtual ~minPatchDistanceRefinement();
// Member Functions
// Selection of refinement/unrefinement candidates
//- Return transferable list of cells to refine
virtual Xfer<labelList> refinementCellCandidates() const;
//- Return transferable list of split points to unrefine
virtual Xfer<labelList> unrefinementPointCandidates() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -119,6 +119,7 @@ $(constraintFvPatchFields)/symmetry/symmetryFvPatchFields.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchFields.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C
$(constraintFvPatchFields)/wedge/wedgeFvPatchVectorNFields.C
$(constraintFvPatchFields)/ggi/ggiFvPatchScalarField.C
$(constraintFvPatchFields)/ggi/ggiFvPatchFields.C
$(constraintFvPatchFields)/ggi/ggiFvPatchVectorNFields.C
$(constraintFvPatchFields)/jumpGgi/jumpGgiFvPatchFields.C

View file

@ -156,15 +156,22 @@ tmp<Field<Type> > cyclicGgiFvPatchField<Type>::patchNeighbourField() const
if (cyclicGgiPatch_.bridgeOverlap())
{
// Symmetry treatment used for overlap
vectorField nHat = this->patch().nf();
// Use mirrored neighbour field for interpolation. Note: mirroring needs
// to take into account the weights, i.e. how "far" we are actually
// mirroring. VV, 19/Jan/2018.
const Field<Type> mirrorField =
transform
(
(I - sqr(this->patch().nf())/
(1.0 - cyclicGgiPatch_.fvPatch::weights())),
this->patchInternalField()
);
// Use mirrored internal field for neighbour
// HJ, 27/Jan/2009
Field<Type> bridgeField =
transform(I - 2.0*sqr(nHat), this->patchInternalField());
// Set mirror values to fully uncovered faces
cyclicGgiPatch_.setUncoveredFaces(mirrorField, pnf);
cyclicGgiPatch_.bridge(bridgeField, pnf);
// Add part of the mirror field to partially covered faces
cyclicGgiPatch_.addToPartialFaces(mirrorField, pnf);
}
return tpnf;
@ -188,19 +195,8 @@ void cyclicGgiFvPatchField<Type>::initEvaluate
+ (1.0 - this->patch().weights())*this->patchNeighbourField()
);
if (cyclicGgiPatch_.bridgeOverlap())
{
// Symmetry treatment used for overlap
vectorField nHat = this->patch().nf();
Field<Type> bridgeField =
(
this->patchInternalField()
+ transform(I - 2.0*sqr(nHat), this->patchInternalField())
)/2.0;
cyclicGgiPatch_.bridge(bridgeField, pf);
}
// Note: bridging and correction of partially overlapping faces taken into
// account in patchNeighbourField(). VV, 16/Oct/2017.
Field<Type>::operator=(pf);
}

View file

@ -38,6 +38,7 @@ Note on parallelisation
#include "ggiFvPatchField.H"
#include "symmTransformField.H"
#include "coeffFields.H"
#include "fvMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -181,25 +182,18 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
// Use mirrored neighbour field for interpolation. Note: mirroring needs
// to take into account the weights, i.e. how "far" we are actually
// mirroring. VV, 19/Jan/2018.
const Field<Type> bridgeField =
const Field<Type> mirrorField =
transform
(
(I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
this->patchInternalField()
);
// if (pTraits<Type>::rank == 0)
// {
// // Scale the field for scalars to ensure conservative and consistent
// // flux on both sides
// ggiPatch_.scaleForPartialCoverage(bridgeField, pnf);
// }
// else
{
// Bridge the field for higher order tensors to correctly take into
// account mirroring
ggiPatch_.bridge(bridgeField, pnf);
}
// Set mirror values to fully uncovered faces
ggiPatch_.setUncoveredFaces(mirrorField, pnf);
// Add part of the mirror field to partially covered faces
ggiPatch_.addToPartialFaces(mirrorField, pnf);
}
return tpnf;
@ -277,14 +271,16 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
// Use mirrored neighbour field for interpolation. Note: mirroring needs
// to take into account the weights, i.e. how "far" we are actually
// mirroring. VV, 19/Jan/2018.
const scalarField bridgeField =
const scalarField mirrorField =
transform
(
(I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
ggiPatch_.patchInternalField(psiInternal)
);
ggiPatch_.bridge(bridgeField, pnf);
// Only set fully uncovered faces. Partially covered faces taken into
// account by manipulating value and gradient matrix coefficients
ggiPatch_.setUncoveredFaces(mirrorField, pnf);
}
// Multiply the field by coefficients and add into the result
@ -353,14 +349,16 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
// Use mirrored neighbour field for interpolation. Note: mirroring needs
// to take into account the weights, i.e. how "far" we are actually
// mirroring. VV, 19/Jan/2018.
const Field<Type> bridgeField =
const Field<Type> mirrorField =
transform
(
(I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
ggiPatch_.patchInternalField(psiInternal)
);
ggiPatch_.bridge(bridgeField, pnf);
// Only set fully uncovered faces. Partially covered faces taken into
// account by manipulating value and gradient matrix coefficients
ggiPatch_.setUncoveredFaces(mirrorField, pnf);
}
// Multiply neighbour field with coeffs and re-use pnf for result
@ -396,7 +394,207 @@ void ggiFvPatchField<Type>::updateInterfaceMatrix
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
{
// Does nothing
}
template<class Type>
void Foam::ggiFvPatchField<Type>::patchFlux
(
GeometricField<Type, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<Type>& matrix
) const
{
// Since we have adjusted the internal/boundary coefficients in the
// manipulateMatrix member function below, we must not use
// patchNeighbourField for reconstructing the flux. We only need to use
// interpolated shadow field. VV, 6/Mar/2018.
// Get patch ID
const label patchI = this->patch().index();
// Get internal field
const Field<Type>& iField = this->internalField();
// Get shadow face-cells and assemble shadow field
const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
Field<Type> sField(sfc.size());
forAll (sField, i)
{
sField[i] = iField[sfc[i]];
}
// Interpolate shadow to this side. Note: must not bridge since internal
// coeffs and boundary coeffs take it into account
Field<Type> neighbourField(ggiPatch_.interpolate(sField));
if (ggiPatch_.bridgeOverlap())
{
const Field<Type> mirrorField =
transform
(
(I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
this->patchInternalField()
);
// Only set fully uncovered faces. Partially covered faces taken into
// account by manipulating value and gradient matrix coefficients
ggiPatch_.setUncoveredFaces(mirrorField, neighbourField);
}
// Calculate the flux with correct neighbour field (fully uncovered faces
// bridged, while partially uncovered faces taken into account by
// manipulating value and gradient matrix coefficients in order to ensure
// conservation for both convection and diffusion part across partially
// overlapping faces). VV, 14/Mar/2018.
flux.boundaryField()[patchI] =
cmptMultiply
(
matrix.internalCoeffs()[patchI],
this->patchInternalField()
)
- cmptMultiply
(
matrix.boundaryCoeffs()[patchI],
neighbourField
);
}
template<class Type>
void Foam::ggiFvPatchField<Type>::manipulateValueCoeffs
(
fvMatrix<Type>& matrix
) const
{
// Conservative treatment for convection across bridged overlap
// Note: master corrects both sets of coefficients (master and slave)
if (ggiPatch_.bridgeOverlap() && ggiPatch_.master())
{
// Get this patch and shadow patch index. In this case, this patch is
// always master and shadow patch is always slave (second condition in
// the if statement)
const label mPatchI = this->patch().index();
const label sPatchI = ggiPatch_.shadowIndex();
// Get all matrix coefficients
Field<Type>& masterIC = matrix.internalCoeffs()[mPatchI];
Field<Type>& masterBC = matrix.boundaryCoeffs()[mPatchI];
Field<Type>& slaveIC = matrix.internalCoeffs()[sPatchI];
Field<Type>& slaveBC = matrix.boundaryCoeffs()[sPatchI];
// Get surface area magnitudes
const scalarField& magSfMaster = ggiPatch_.magSf();
const scalarField& magSfSlave = ggiPatch_.shadow().magSf();
// Interpolate all coeffs from each side on the other side. Note: need
// to normalize by surface areas since it is "hidden" inside the
// convection coefficient (usually Sf & Uf), before interpolation. After
// interpolation, need to multiply with surface areas on this side.
// Interpolation from slave to this (master) side
const Field<Type> masterInterpolatedIC =
magSfSlave*ggiPatch_.shadow().interpolate(masterIC/magSfMaster);
const Field<Type> masterInterpolatedBC =
magSfSlave*ggiPatch_.shadow().interpolate(masterBC/magSfMaster);
// Interpolation from this (master) side to slave
const Field<Type> slaveInterpolatedIC =
magSfMaster*ggiPatch_.interpolate(slaveIC/magSfSlave);
const Field<Type> slaveInterpolatedBC =
magSfMaster*ggiPatch_.interpolate(slaveBC/magSfSlave);
// Set partially covered master coeffs using slave data
ggiPatch_.setPartialFaces(slaveInterpolatedBC, masterIC);
ggiPatch_.setPartialFaces(slaveInterpolatedIC, masterBC);
// Scale partially overlapping boundary coeffs to ensure conservation
ggiPatch_.scalePartialFaces(masterBC);
// Set partially covered slave coeffs using master data
ggiPatch_.shadow().setPartialFaces(masterInterpolatedBC, slaveIC);
ggiPatch_.shadow().setPartialFaces(masterInterpolatedIC, slaveBC);
// Scale partially overlapping boundary coeffs to ensure conservation
ggiPatch_.shadow().scalePartialFaces(slaveBC);
}
}
template<class Type>
void Foam::ggiFvPatchField<Type>::manipulateGradientCoeffs
(
fvMatrix<Type>& matrix
) const
{
// Conservative treatment for convection across bridged overlap
// Note: master corrects both sets of coefficients (master and slave)
if (ggiPatch_.bridgeOverlap() && ggiPatch_.master())
{
// Get this patch and shadow patch index. In this case, this patch is
// always master and shadow patch is always slave (second condition in
// the if statement)
const label mPatchI = this->patch().index();
const label sPatchI = ggiPatch_.shadowIndex();
// Get all matrix coefficients
Field<Type>& masterIC = matrix.internalCoeffs()[mPatchI];
Field<Type>& masterBC = matrix.boundaryCoeffs()[mPatchI];
Field<Type>& slaveIC = matrix.internalCoeffs()[sPatchI];
Field<Type>& slaveBC = matrix.boundaryCoeffs()[sPatchI];
// Get surface area magnitudes
const scalarField& magSfMaster = ggiPatch_.magSf();
const scalarField& magSfSlave = ggiPatch_.shadow().magSf();
// Interpolate all coeffs from each side on the other side. Note: need
// to normalize by surface areas since it is "hidden" inside the
// diffusion coefficient (usually |Sf|/|df|), before interpolation.
// After interpolation, need to multiply with surface areas on this
// side.
// Interpolation from slave to this (master) side
const Field<Type> masterInterpolatedIC =
magSfSlave*ggiPatch_.shadow().interpolate(masterIC/magSfMaster);
const Field<Type> masterInterpolatedBC =
magSfSlave*ggiPatch_.shadow().interpolate(masterBC/magSfMaster);
const Field<Type> slaveInterpolatedIC =
magSfMaster*ggiPatch_.interpolate(slaveIC/magSfSlave);
const Field<Type> slaveInterpolatedBC =
magSfMaster*ggiPatch_.interpolate(slaveBC/magSfSlave);
// Set partially covered master coeffs using slave data
ggiPatch_.setPartialFaces(slaveInterpolatedBC, masterIC);
ggiPatch_.setPartialFaces(slaveInterpolatedIC, masterBC);
// Scale partially overlapping master coeffs to ensure conservation
ggiPatch_.scalePartialFaces(masterIC);
ggiPatch_.scalePartialFaces(masterBC);
// Add to partially overlapping master internal coeffs to ensure
// conservation. Note: boundary coeffs must be negated
ggiPatch_.addToPartialFaces((-masterBC)(), masterIC);
// Set partially covered slave coeffs using master data
ggiPatch_.shadow().setPartialFaces(masterInterpolatedBC, slaveIC);
ggiPatch_.shadow().setPartialFaces(masterInterpolatedIC, slaveBC);
// Scale partially overlapping slave coeffs to ensure conservation
ggiPatch_.shadow().scalePartialFaces(slaveIC);
ggiPatch_.shadow().scalePartialFaces(slaveBC);
// Add to partially overlapping slave internal coeffs to ensure
// conservation. Note: boundary coeffs must be negated
ggiPatch_.shadow().addToPartialFaces((-slaveBC)(), slaveIC);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -47,6 +47,11 @@ SourceFiles
namespace Foam
{
// Forward declaration
template<class Type>
class fvMatrix;
/*---------------------------------------------------------------------------*\
Class ggiFvPatchField Declaration
\*---------------------------------------------------------------------------*/
@ -228,6 +233,29 @@ public:
) const;
// Patch interpolation and patch flux
//- Calculate patch flux
virtual void patchFlux
(
GeometricField<Type, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<Type>& matrix
) const;
// Matrix manipulation
//- Manipulate value and boundary coefficients for convection. Used
// to ensure conservation across partially overlapping GGI faces if
// the bridging is switched on
virtual void manipulateValueCoeffs(fvMatrix<Type>& matrix) const;
//- Manipulate value and boundary coefficients for diffusion. Used
// to ensure conservation across partially overlapping GGI faces if
// the bridging is switched on
virtual void manipulateGradientCoeffs(fvMatrix<Type>& matrix) const;
// GGI coupled interface functions
//- Does the patch field perform the transfromation
@ -262,6 +290,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "ggiFvPatchScalarField.H"
#ifdef NoRepository
# include "ggiFvPatchField.C"
#endif

View file

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved
Description
Specilalisation of patchFlux member function for scalars needed to
ensure conservation across GGI patches that are partially covered if the
bridge overlap switch is on.
\*---------------------------------------------------------------------------*/
#include "ggiFvPatchField.H"
#include "surfaceFields.H"
#include "fvScalarMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<>
void Foam::ggiFvPatchField<scalar>::patchFlux
(
GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<scalar>& matrix
) const
{
// Since we have adjusted the internal/boundary coefficients in the
// manipulateMatrix member function below, we must not use
// patchNeighbourField for reconstructing the flux. We only need to use
// interpolated shadow field. VV, 6/Mar/2018.
// Get patch ID
const label patchI = this->patch().index();
// Get internal field
const scalarField& iField = this->internalField();
// Get shadow face-cells and assemble shadow field
const unallocLabelList& sfc = ggiPatch_.shadow().faceCells();
scalarField sField(sfc.size());
forAll (sField, i)
{
sField[i] = iField[sfc[i]];
}
// Interpolate shadow to this side. Note: must not bridge since internal
// coeffs and boundary coeffs take it into account
scalarField neighbourField(ggiPatch_.interpolate(sField));
if (ggiPatch_.bridgeOverlap())
{
// Only set fully uncovered faces. Partially covered faces taken into
// account by manipulating value and gradient matrix coefficients. Note:
// mirror field is the same as patch internal field
ggiPatch_.setUncoveredFaces
(
this->patchInternalField()(),
neighbourField
);
}
// Calculate the flux with correct neighbour field (fully uncovered faces
// bridged, while partially uncovered faces taken into account by
// manipulating value and gradient matrix coefficients in order to ensure
// conservation for both convection and diffusion part across partially
// overlapping faces). VV, 14/Mar/2018.
flux.boundaryField()[patchI] =
matrix.internalCoeffs()[patchI]*this->patchInternalField()
- matrix.boundaryCoeffs()[patchI]*neighbourField;
// Scale the flux on slave patch to ensure global conservation across this
// partially overlapping GGI patch. The slight disbalance can happen since
// we interpolate the matrix coeffs and field separately, and we should do
// it together to ensure conservation. Current code design does not easily
// allow this, so this is a meaningful temporary solution when we have
// partially overlapping faces. VV, 14/Mar/2018.
if (ggiPatch_.bridgeOverlap() && !ggiPatch_.master())
{
// This is slave patch, master already updated the fluxes and we can use
// that info to scale the fluxes on this side
// Get the total flux through master patch
const scalar masterFlux =
gSum(flux.boundaryField()[ggiPatch_.shadowIndex()]);
// Get the total flux through slave patch
const scalar slaveFlux = gSum(flux.boundaryField()[patchI]);
// Calculate the scaling factor. Note: negative sign since master and
// slave fluxes have opposite sign
const scalar scalingFactor = -masterFlux/(slaveFlux + SMALL);
// Scale the slave flux to ensure global patch conservation
flux.boundaryField()[patchI] *= scalingFactor;
if (debug)
{
Info<< "Scaling flux on patch: " << ggiPatch_.name()
<< " with " << scalingFactor
<< " to ensure conservation." << endl;
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,65 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved
Description
Specilalisation of patchFlux member function for scalars needed to
ensure conservation across GGI patches that are partially covered if the
bridge overlap switch is on.
SourceFiles
ggiFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef ggiFvPatchScalarField_H
#define ggiFvPatchScalarField_H
#include "ggiFvPatchField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<>
void Foam::ggiFvPatchField<scalar>::patchFlux
(
GeometricField<scalar, fvsPatchField, surfaceMesh>& flux,
const fvMatrix<scalar>& matrix
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -229,7 +229,7 @@ tmp<Field<Type> > regionCouplingFvPatchField<Type>::patchNeighbourField() const
Field<Type> bridgeField =
transform(I - 2.0*sqr(nHat), this->patchInternalField());
regionCouplePatch_.bridge(bridgeField, pnf);
regionCouplePatch_.setUncoveredFaces(bridgeField, pnf);
}
return tpnf;
@ -304,7 +304,7 @@ void regionCouplingFvPatchField<Type>::initEvaluate
Field<Type> bridgeField =
0.5*(pif + transform(I - 2.0*sqr(nHat), pif));
regionCouplePatch_.bridge(bridgeField, *this);
regionCouplePatch_.setUncoveredFaces(bridgeField, *this);
}
}
@ -357,7 +357,7 @@ void regionCouplingFvPatchField<Type>::updateCoeffs()
Field<Type> bridgeField =
0.5*(pif + transform(I - 2.0*sqr(nHat), pif));
regionCouplePatch_.bridge(bridgeField, *this);
regionCouplePatch_.setUncoveredFaces(bridgeField, *this);
}
}

View file

@ -248,6 +248,26 @@ void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateValueCoeffs
(
fvMatrix<Type>& matrix
) const
{
// do nothing
}
template<class Type>
void Foam::fvPatchField<Type>::manipulateGradientCoeffs
(
fvMatrix<Type>& matrix
) const
{
// do nothing
}
template<class Type>
void Foam::fvPatchField<Type>::patchInterpolate
(

View file

@ -465,6 +465,16 @@ public:
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
//- Manipulate value and boundary coefficients for convection.
// Called from gaussConvection scheme and does nothing by default.
// Currently overriden in ggiFvPatchField only
virtual void manipulateValueCoeffs(fvMatrix<Type>& matrix) const;
//- Manipulate value and boundary coefficients for diffusion.
// Called from gaussConvection scheme and does nothing by default.
// Currently overriden in ggiFvPatchField only
virtual void manipulateGradientCoeffs(fvMatrix<Type>& matrix) const;
// I-O

View file

@ -98,6 +98,15 @@ gaussConvectionScheme<Type>::fvmDiv
fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);
}
// Manipulate internal and boundary coeffs for convection. Needed for very
// special treatment and is currently used only for ensuring implicit
// conservation across GGI interface that has partially covered faces. Does
// nothing for other fvPatchFields. VV, 8/Mar/2018.
forAll(fvm.psi().boundaryField(), patchI)
{
fvm.psi().boundaryField()[patchI].manipulateValueCoeffs(fvm);
}
if (tinterpScheme_().corrected())
{
fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));

View file

@ -76,6 +76,15 @@ gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
fvm.boundaryCoeffs()[patchI] = -patchGamma*psf.gradientBoundaryCoeffs();
}
// Manipulate internal and boundary coeffs for diffusion. Needed for very
// special treatment and is currently used only for ensuring implicit
// conservation across GGI interface that has partially covered faces. Does
// nothing for other fvPatchFields. VV, 8/Mar/2018.
forAll(fvm.psi().boundaryField(), patchI)
{
fvm.psi().boundaryField()[patchI].manipulateGradientCoeffs(fvm);
}
return tfvm;
}

View file

@ -52,42 +52,48 @@ void Foam::cyclicGgiFvPatch::makeWeights(scalarField& w) const
// HJ, 2/Aug/2007
if (cyclicGgiPolyPatch_.master())
{
vectorField n = nf();
// Master side. No need to scale partially uncovered or set fully
// uncovered faces since delta already takes it into account.
// VV, 25/Feb/2018.
const vectorField n = nf();
// Note: mag in the dot-product.
// For all valid meshes, the non-orthogonality will be less than
// 90 deg and the dot-product will be positive. For invalid
// meshes (d & s <= 0), this will stabilise the calculation
// but the result will be poor. HJ, 24/Aug/2011
scalarField nfc =
mag
(
n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf())
);
const scalarField nfc =
mag(n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf()));
w = nfc/(mag(n & (Cf() - Cn())) + nfc);
if (bridgeOverlap())
{
// Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w);
}
w = nfc/(mag(n & (Cf() - Cn())) + nfc + SMALL);
}
else
{
// Pick up weights from the master side
// Slave side. Interpolate the master side weights, scale them for
// partially covered faces and set weights for fully uncovered faces if
// the bridge overlap is switched on. VV, 15/Feb/2018.
scalarField masterWeights(shadow().size());
shadow().makeWeights(masterWeights);
w = interpolate(1 - masterWeights);
// Interpolate master weights to this side
w = interpolate(masterWeights);
if (bridgeOverlap())
{
// Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w);
// Weights for fully uncovered faces
const scalarField uncoveredWeights(w.size(), 0.5);
// Set weights for uncovered faces
setUncoveredFaces(uncoveredWeights, w);
// Scale partially overlapping faces
scalePartialFaces(w);
}
// Finally construct these weights as 1 - master weights
w = 1 - w;
}
}
@ -97,29 +103,38 @@ void Foam::cyclicGgiFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
if (cyclicGgiPolyPatch_.master())
{
// Master side. No need to scale partially uncovered or set fully
// uncovered faces since delta already takes it into account.
// VV, 25/Feb/2018.
// Stabilised form for bad meshes. HJ, 24/Aug/2011
vectorField d = delta();
const vectorField d = delta();
dc = 1.0/max(nf() & d, 0.05*mag(d));
if (bridgeOverlap())
{
scalarField bridgeDeltas = nf() & fvPatch::delta();
bridge(bridgeDeltas, dc);
}
}
else
{
// Slave side. Interpolate the master side, scale it for partially
// covered faces and set deltaCoeffs for fully uncovered faces if the
// bridge overlap is switched on. VV, 15/Feb/2018.
scalarField masterDeltas(shadow().size());
shadow().makeDeltaCoeffs(masterDeltas);
dc = interpolate(masterDeltas);
if (bridgeOverlap())
{
scalarField bridgeDeltas = nf() & fvPatch::delta();
// Delta coeffs for fully uncovered faces obtained from deltas on
// this side
const vectorField d = delta();
const scalarField uncoveredDeltaCoeffs =
1.0/max(nf() & d, 0.05*mag(d));
bridge(bridgeDeltas, dc);
// Set delta coeffs for uncovered faces
setUncoveredFaces(uncoveredDeltaCoeffs, dc);
// Scale partially overlapping faces
scalePartialFaces(dc);
}
}
}
@ -137,33 +152,39 @@ Foam::tmp<Foam::vectorField> Foam::cyclicGgiFvPatch::delta() const
{
if (cyclicGgiPolyPatch_.master())
{
tmp<vectorField> tDelta =
// Master side. Note: scaling partially covered faces and setting deltas
// to fully uncovered faces correctly taken into account in
// reconFaceCellCentres function. VV, 15/Feb/2018.
tmp<vectorField> tdelta =
cyclicGgiPolyPatch_.reconFaceCellCentres() - Cn();
if (bridgeOverlap())
{
vectorField bridgeDeltas = Cf() - Cn();
bridge(bridgeDeltas, tDelta());
}
return tDelta;
return tdelta;
}
else
{
tmp<vectorField> tDelta = interpolate
// Slave side. Interpolate the master side, scale it for partially
// covered faces and set deltas for fully uncovered faces if the bridge
// overlap is switched on. VV, 15/Feb/2018.
tmp<vectorField> tdelta = interpolate
(
shadow().Cn() - cyclicGgiPolyPatch_.shadow().reconFaceCellCentres()
);
if (bridgeOverlap())
{
vectorField bridgeDeltas = Cf() - Cn();
// Deltas for fully uncovered faces
const vectorField uncoveredDeltas(2.0*fvPatch::delta());
bridge(bridgeDeltas, tDelta());
// Set deltas for fully uncovered faces
setUncoveredFaces(uncoveredDeltas, tdelta());
// Scale for partially covered faces
scalePartialFaces(tdelta());
}
return tDelta;
return tdelta;
}
}

View file

@ -95,9 +95,11 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
// Weights for fully uncovered faces
const scalarField uncoveredWeights(w.size(), 0.5);
// Scale partially overlapping faces and set uncovered weights
// for fully uncovered faces
scaleForPartialCoverage(uncoveredWeights, w);
// Set weights for uncovered faces
setUncoveredFaces(uncoveredWeights, w);
// Scale partially overlapping faces
scalePartialFaces(w);
}
// Finally construct these weights as 1 - master weights
@ -138,9 +140,11 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
const scalarField uncoveredDeltaCoeffs =
1.0/max(nf() & d, 0.05*mag(d));
// Scale partially overlapping faces and set uncovered deltaCoeffs
// for fully uncovered faces.
scaleForPartialCoverage(uncoveredDeltaCoeffs, dc);
// Set delta coeffs for uncovered faces
setUncoveredFaces(uncoveredDeltaCoeffs, dc);
// Scale partially overlapping faces
scalePartialFaces(dc);
}
}
}
@ -152,15 +156,24 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const
// Non-orthogonality correction on a ggi interface
// MB, 7/April/2009
// Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs();
// No non-orthognal correction if the bridge overlap is switched on to
// ensure conservative interpolation for partially overlapping faces
if (bridgeOverlap())
{
cv = vector::zero;
}
else
{
// Calculate correction vectors on coupled patches
const scalarField& patchDeltaCoeffs = deltaCoeffs();
const vectorField patchDeltas = delta();
const vectorField n = nf();
const vectorField patchDeltas = delta();
const vectorField n = nf();
// If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011
cv = pos(patchDeltas & n)*(n - patchDeltas*patchDeltaCoeffs);
// If non-orthogonality is over 90 deg, kill correction vector
// HJ, 6/Jan/2011
cv = pos(patchDeltas & n)*(n - patchDeltas*patchDeltaCoeffs);
}
}
@ -173,9 +186,9 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
// to fully uncovered faces correctly taken into account in
// reconFaceCellCentres function. VV, 15/Feb/2018.
tmp<vectorField> tDelta = ggiPolyPatch_.reconFaceCellCentres() - Cn();
tmp<vectorField> tdelta = ggiPolyPatch_.reconFaceCellCentres() - Cn();
return tDelta;
return tdelta;
}
else
{
@ -183,7 +196,7 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
// covered faces and set deltas for fully uncovered faces if the bridge
// overlap is switched on. VV, 15/Feb/2018.
tmp<vectorField> tDelta = interpolate
tmp<vectorField> tdelta = interpolate
(
shadow().Cn() - ggiPolyPatch_.shadow().reconFaceCellCentres()
);
@ -193,12 +206,14 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
// Deltas for fully uncovered faces
const vectorField uncoveredDeltas(2.0*fvPatch::delta());
// Scale partially overlapping faces and set uncovered deltas for
// fully uncovered faces
scaleForPartialCoverage(uncoveredDeltas, tDelta());
// Set deltas for fully uncovered faces
setUncoveredFaces(uncoveredDeltas, tdelta());
// Scale for partially covered faces
scalePartialFaces(tdelta());
}
return tDelta;
return tdelta;
}
}

View file

@ -115,6 +115,12 @@ public:
// Access
//- Return ggiPolyPatch
const ggiPolyPatch& ggiPatch() const
{
return ggiPolyPatch_;
}
//- Return shadow patch
const ggiFvPatch& shadow() const;
@ -128,7 +134,8 @@ public:
virtual tmp<vectorField> delta() const;
// Interpolation
// Interpolation functions and bridging operations for fully uncovered
// and partially covered faces
//- Interpolate face field
template<class Type>
@ -143,33 +150,48 @@ public:
return ggiPolyPatch_.interpolate(tpf);
}
//- Bridge interpolated face field for uncovered faces
//- Set given field for uncovered faces. Usually used for setting
// mirrored field
template<class Type>
void bridge
void setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
return ggiPolyPatch_.bridge(bridgeField, ff);
return ggiPolyPatch_.setUncoveredFaces(fieldToSet, ff);
}
//- Scale field for partially covered faces and set
// uncoveredFaceField to fully uncovered faces. Needed for correct
// calculation of fvPatch data: delta, deltaCoeffs and weights
//- Set given field for partially uncovered faces
template<class Type>
void scaleForPartialCoverage
void setPartialFaces
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
return
ggiPolyPatch_.scaleForPartialCoverage
(
uncoveredFaceField,
ff
);
return ggiPolyPatch_.setPartialFaces(fieldToSet, ff);
}
//- Scale field for partially covered faces. Needed for correct
// construction of mesh interpolation data (weights and delta
// coeffs) and to ensure fully conservative implicit treatment
template<class Type>
void scalePartialFaces(Field<Type>& ff) const
{
return ggiPolyPatch_.scalePartialFaces(ff);
}
//- Add given field to partially covered faces. Needed to ensure
// fully conservative implicit treatment
template<class Type>
void addToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const
{
return ggiPolyPatch_.addToPartialFaces(fieldToAdd, ff);
}

View file

@ -76,7 +76,7 @@ void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const
{
// Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w);
setUncoveredFaces(scalarField(size(), 0.5), w);
}
}
else
@ -93,7 +93,7 @@ void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const
{
// Set overlap weights to 0.5 and use mirrored neighbour field
// for interpolation. HJ, 21/Jan/2009
bridge(scalarField(size(), 0.5), w);
setUncoveredFaces(scalarField(size(), 0.5), w);
}
}
}
@ -120,7 +120,7 @@ void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
scalarField bridgeDeltas = nf() & fvPatch::delta();
bridge(bridgeDeltas, dc);
setUncoveredFaces(bridgeDeltas, dc);
}
}
else
@ -133,7 +133,7 @@ void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
scalarField bridgeDeltas = nf() & fvPatch::delta();
bridge(bridgeDeltas, dc);
setUncoveredFaces(bridgeDeltas, dc);
}
}
}
@ -184,7 +184,7 @@ Foam::tmp<Foam::vectorField> Foam::regionCoupleFvPatch::delta() const
{
vectorField bridgeDeltas = Cf() - Cn();
bridge(bridgeDeltas, tDelta());
setUncoveredFaces(bridgeDeltas, tDelta());
}
return tDelta;
@ -200,7 +200,7 @@ Foam::tmp<Foam::vectorField> Foam::regionCoupleFvPatch::delta() const
{
vectorField bridgeDeltas = Cf() - Cn();
bridge(bridgeDeltas, tDelta());
setUncoveredFaces(bridgeDeltas, tDelta());
}
return tDelta;

View file

@ -152,15 +152,16 @@ public:
return rcPolyPatch_.interpolate(tpf);
}
//- Bridge interpolated face field for uncovered faces
//- Set given field for uncovered faces. Usually used for setting
// mirrored field
template<class Type>
void bridge
void setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
return rcPolyPatch_.bridge(bridgeField, ff);
return rcPolyPatch_.setUncoveredFaces(fieldToSet, ff);
}

View file

@ -93,97 +93,55 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::bridge
void GGIInterpolation<MasterPatch, SlavePatch>::setFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
const labelList& facesToSet
)
{
// Fully uncovered faces
forAll (addr, faceI)
// Loop throuh fully uncovered faces and set the coefficients
forAll (facesToSet, ftsI)
{
ff[addr[faceI]] = bridgeField[addr[faceI]];
}
// Get face index
const label& faceI = facesToSet[ftsI];
// Loop through partially covered faces and correct them. Note the
// operator+= since we assume that the interpolation part is carried out
// before bridging (see e.g. ggiFvPatchField::patchNeighbourField()) using
// weights that do not sum up to 1
forAll (partiallyUncoveredAddr, pcfI)
{
ff[partiallyUncoveredAddr[pcfI]] +=
uncoveredFractions[pcfI]*bridgeField[partiallyUncoveredAddr[pcfI]];
// Set field for this face
ff[faceI] = fieldToSet[faceI];
}
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
void GGIInterpolation<MasterPatch, SlavePatch>::maskedSetFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
const labelList& facesToSet
)
{
// Note: tricky algorithm
// In order for a face to be bridged it needs to be both in the
// mask and in selection of faces that are bridged (addr).
// In order for a face to be set it needs to be both in the
// mask and in selection of faces that are set (addr).
// This implies an n-squared search, but we can avoid it by
// using the fact that both lists are ordered.
label maskAddrI = 0;
forAll (uncoveredFaces, uncoI)
forAll (facesToSet, ftsI)
{
// Pick the uncovered face
const label faceI = uncoveredFaces[uncoI];
// Pick the face index
const label& faceI = facesToSet[ftsI];
// Search through the mask
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked bridged face
// Put the result into condensed list: masked faces only
ff[maskAddrI] = bridgeField[maskAddrI];
break;
}
else if (mask[maskAddrI] > faceI)
{
// Gone beyond my index: my face is not present in the mask
// Go one back and check for next uncovered face
if (maskAddrI > 0)
{
--maskAddrI;
}
break;
}
}
}
// Reset maskAddrI
maskAddrI = 0;
forAll (partiallyUncoveredAddr, pcfI)
{
// Pick partially covered face
const label faceI = partiallyUncoveredAddr[pcfI];
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked partially covered face
ff[maskAddrI] += uncoveredFractions[pcfI]*bridgeField[maskAddrI];
// Found masked face, set the field
ff[maskAddrI] = fieldToSet[maskAddrI];
break;
}
@ -205,22 +163,13 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scale
void GGIInterpolation<MasterPatch, SlavePatch>::scalePartial
(
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
)
{
// Loop through uncovered faces and set them
forAll (addr, ufI)
{
const label& faceI = addr[ufI];
ff[faceI] = uncoveredFaceField[faceI];
}
// Loop through partially covered faces and scale them up
forAll (partiallyUncoveredAddr, pcfI)
{
@ -232,37 +181,33 @@ void GGIInterpolation<MasterPatch, SlavePatch>::scale
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScale
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScalePartial
(
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
)
{
// Note: tricky algorithm
// In order for a face to be bridged it needs to be both in the
// mask and in selection of faces that are bridged (addr).
// In order for a face to be set it needs to be both in the
// mask and in selection of faces that are set (addr).
// This implies an n-squared search, but we can avoid it by
// using the fact that both lists are ordered.
label maskAddrI = 0;
forAll (uncoveredFaces, uncoI)
forAll (partiallyUncoveredAddr, pcfI)
{
// Pick the uncovered face
const label faceI = uncoveredFaces[uncoI];
// Pick partially covered face
const label& faceI = partiallyUncoveredAddr[pcfI];
// Search through the mask
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked bridged face
// Put the result into condensed list: masked faces only
ff[maskAddrI] = uncoveredFaceField[maskAddrI];
// Found masked partially covered face, scale it up
ff[maskAddrI] /= (1.0 - uncoveredFractions[pcfI]);
break;
}
@ -279,21 +224,61 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScale
}
}
}
}
// Reset maskAddrI
maskAddrI = 0;
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::addToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
)
{
// Loop through partially covered faces and add the weighted field
forAll (partiallyUncoveredAddr, pcfI)
{
// Get face index
const label& faceI = partiallyUncoveredAddr[pcfI];
// Add to partially covered face
ff[faceI] += uncoveredFractions[pcfI]*fieldToAdd[faceI];
}
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedAddToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
)
{
// Note: tricky algorithm
// In order for a face to be partially covered it needs to be both in the
// mask and in selection of faces that are partially covered. This implies
// an n-squared search, but we can avoid it by using the fact that both
// lists are ordered.
label maskAddrI = 0;
forAll (partiallyUncoveredAddr, pcfI)
{
// Pick partially covered face
const label faceI = partiallyUncoveredAddr[pcfI];
const label& faceI = partiallyUncoveredAddr[pcfI];
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked partially covered face
ff[maskAddrI] /= (1.0 - uncoveredFractions[pcfI]);
// Found masked partially covered face, add to it
ff[maskAddrI] += uncoveredFractions[pcfI]*fieldToAdd[maskAddrI];
break;
}
@ -658,48 +643,47 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedSlaveToMaster
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
void GGIInterpolation<MasterPatch, SlavePatch>::setUncoveredFacesMaster
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
if
(
(bridgeField.size() != masterPatch_.size())
(fieldToSet.size() != masterPatch_.size())
|| (ff.size() != masterPatch_.size())
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster\n"
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"setUncoveredFacesMaster\n"
"(\n"
" const Field<Type>& bridgeField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< abort(FatalError);
}
bridge
setFaces
(
bridgeField,
fieldToSet,
ff,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
uncoveredMasterFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
void GGIInterpolation<MasterPatch, SlavePatch>::maskedSetUncoveredFacesMaster
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const
@ -709,77 +693,73 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedBridgeMaster\n"
"maskedSetUncoveredFacesMaster\n"
"(\n"
" const Field<Type>& bridgeField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedBridge
maskedSetFaces
(
bridgeField,
fieldToSet,
ff,
mask,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
uncoveredMasterFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
void GGIInterpolation<MasterPatch, SlavePatch>::setUncoveredFacesSlave
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
if
(
(bridgeField.size() != slavePatch_.size())
(fieldToSet.size() != slavePatch_.size())
|| (ff.size() != slavePatch_.size())
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"bridgeSlave\n"
"setUncoveredFacesSlave\n"
"(\n"
" const Field<Type>& bridgeField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff"
") const"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< abort(FatalError);
}
bridge
setFaces
(
bridgeField,
fieldToSet,
ff,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
uncoveredSlaveFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
void GGIInterpolation<MasterPatch, SlavePatch>::maskedSetUncoveredFacesSlave
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const
@ -789,76 +769,73 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedBridgeSlave\n"
"maskedSetUncoveredFacesSlave\n"
"(\n"
" const Field<Type>& bridgeField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " bridge field size: " << bridgeField.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedBridge
maskedSetFaces
(
bridgeField,
fieldToSet,
ff,
mask,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
uncoveredSlaveFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scaleMaster
void GGIInterpolation<MasterPatch, SlavePatch>::setPartialFacesMaster
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
if
(
uncoveredFaceField.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size()
(fieldToSet.size() != masterPatch_.size())
|| (ff.size() != masterPatch_.size())
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::scaleMaster\n"
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"setPartialFacesMaster\n"
"(\n"
" const Field<Type>& uncoveredFaceField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< abort(FatalError);
}
scale
setFaces
(
uncoveredFaceField,
fieldToSet,
ff,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
partiallyUncoveredMasterFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleMaster
void GGIInterpolation<MasterPatch, SlavePatch>::maskedSetPartialFacesMaster
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const
@ -868,26 +845,131 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleMaster
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScaleMaster\n"
"maskedSetPartialFacesMaster\n"
"(\n"
" const Field<Type>& uncoveredFaceField\n,"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScale
maskedSetFaces
(
uncoveredFaceField,
fieldToSet,
ff,
mask,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::setPartialFacesSlave
(
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
if
(
(fieldToSet.size() != slavePatch_.size())
|| (ff.size() != slavePatch_.size())
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"setPartialFacesSlave\n"
"(\n"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff"
") const"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< abort(FatalError);
}
setFaces
(
fieldToSet,
ff,
partiallyUncoveredSlaveFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedSetPartialFacesSlave
(
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const
{
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedSetPartialFacesSlave\n"
"(\n"
" const Field<Type>& fieldToSet\n,"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " field to set size: " << fieldToSet.size()
<< " field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedSetFaces
(
fieldToSet,
ff,
mask,
partiallyUncoveredSlaveFaces()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scalePartialMaster
(
Field<Type>& ff
) const
{
if (ff.size() != masterPatch_.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"scalePartialMaster\n"
"(\n"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
scalePartial
(
ff,
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
);
@ -896,38 +978,216 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleMaster
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scaleSlave
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScalePartialMaster
(
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask
) const
{
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScalePartialMaster\n"
"(\n"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScalePartial
(
ff,
mask,
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scalePartialSlave
(
Field<Type>& ff
) const
{
if (ff.size() != slavePatch_.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"scalePartialSlave\n"
"(\n"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
scalePartial
(
ff,
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScalePartialSlave
(
Field<Type>& ff,
const labelList& mask
) const
{
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScalePartialSlave\n"
"(\n"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScalePartial
(
ff,
mask,
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::addToPartialFacesMaster
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const
{
if
(
uncoveredFaceField.size() != slavePatch_.size()
fieldToAdd.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size()
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"addToPartialFacesMaster\n"
"(\n"
" const Field<Type>& fieldToAdd\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size()
<< " field to add size: " << fieldToAdd.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
addToPartialFaces
(
fieldToAdd,
ff,
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedAddToPartialFacesMaster
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask
) const
{
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedAddToPartialFacesMaster\n"
"(\n"
" const Field<Type>& fieldToAdd\n,"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " field to set size: " << fieldToAdd.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedAddToPartialFaces
(
fieldToAdd,
ff,
mask,
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
);
}
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::addToPartialFacesSlave
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const
{
if
(
fieldToAdd.size() != slavePatch_.size()
|| ff.size() != slavePatch_.size()
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"scaleSlave\n"
"addToPartialFacesSlave\n"
"(\n"
" const Field<Type>& uncoveredFaceField\n,"
" const Field<Type>& fieldToAdd\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " field to add size: " << fieldToAdd.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
scale
addToPartialFaces
(
uncoveredFaceField,
fieldToAdd,
ff,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
);
@ -936,9 +1196,9 @@ void GGIInterpolation<MasterPatch, SlavePatch>::scaleSlave
template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleSlave
void GGIInterpolation<MasterPatch, SlavePatch>::maskedAddToPartialFacesSlave
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask
) const
@ -948,26 +1208,25 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleSlave
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScaleSlave\n"
"maskedAddToPartialFacesSlave\n"
"(\n"
" const Field<Type>& uncoveredFaceField\n,"
" const Field<Type>& fieldToAdd\n,"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " field to add size: " << fieldToAdd.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScale
maskedAddToPartialFaces
(
uncoveredFaceField,
fieldToAdd,
ff,
mask,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
);

View file

@ -478,56 +478,69 @@ class GGIInterpolation
);
// Bridging operations, taking into account partially covered faces
// Operations for setting face coeffs
//- Bridge uncovered faces given addressing
//- Set face coeffs given addressing
template<class Type>
static void bridge
static void setFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
const labelList& facesToSet
);
//- Bridge uncovered faces given addressing
// for masked faces only
//- Set face coeffs given masked addressing
template<class Type>
static void maskedBridge
static void maskedSetFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
const labelList& facesToSet
);
// Scaling operations for partially covered faces
//- Scale partially covered faces and set fully uncovered faces
// given addressing
// Scaling operations for correction of partially covered faces
//- Scale partially covered faces given addressing
template<class Type>
static void scale
static void scalePartial
(
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
);
//- Scale partially covered faces and set fully uncovered faces
// given addressing for masked faces only
//- Scale partially covered faces given masked addressing
template<class Type>
static void maskedScale
static void maskedScalePartial
(
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
);
// Adding operations for partially covered faces
//- Add field weighted by uncoverage factor given addressing
template<class Type>
static void addToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
);
//- Add field weighted by uncoverage factor given masked addressing
template<class Type>
static void maskedAddToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
);
@ -680,81 +693,142 @@ public:
) const;
// Bridging operations taking into account partially uncovered faces
// Operations for setting fully uncovered faces
//- Bridge uncovered master patch field
//- Set uncovered faces for master
template<class Type>
void bridgeMaster
void setUncoveredFacesMaster
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Bridge uncovered master patch field, only for marked master
// faces
//- Set uncovered faces for master, only for marked faces
template<class Type>
void maskedBridgeMaster
void maskedSetUncoveredFacesMaster
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const;
//- Bridge uncovered slave patch field
//- Set uncovered faces for slave
template<class Type>
void bridgeSlave
void setUncoveredFacesSlave
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Bridge uncovered slave patch field, only for marked slave faces
//- Set uncovered faces for slave, only for marked faces
template<class Type>
void maskedBridgeSlave
void maskedSetUncoveredFacesSlave
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const;
// Scaling operations for correction of partially uncovered faces
//- Scale partially uncovered master patch field and set fully
// uncovered face field
// Operations for setting partially covered faces
//- Set partially covered faces for master
template<class Type>
void scaleMaster
void setPartialFacesMaster
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Scale partially uncovered master patch field and set fully
// uncovered face field, only for marked master faces
//- Set partially covered faces for master, only for marked faces
template<class Type>
void maskedScaleMaster
void maskedSetPartialFacesMaster
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const;
//- Scale partially uncovered slave patch field and set fully
// uncovered face field
//- Set partially covered faces for slave
template<class Type>
void scaleSlave
void setPartialFacesSlave
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Scale partially uncovered slave patch field and set fully
// uncovered, only for marked slave faces
//- Set partially covered faces for slave, only for marked faces
template<class Type>
void maskedScaleSlave
void maskedSetPartialFacesSlave
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff,
const labelList& mask
) const;
// Scaling operations for correction of partially covered faces
//- Scale partially covered master faces
template<class Type>
void scalePartialMaster(Field<Type>& ff) const;
//- Scale partially covered master faces, only for marked faces
template<class Type>
void maskedScalePartialMaster
(
Field<Type>& ff,
const labelList& mask
) const;
//- Scale partially covered slave faces
template<class Type>
void scalePartialSlave(Field<Type>& ff) const;
//- Scale partially covered slave face, only for marked slave faces
template<class Type>
void maskedScalePartialSlave
(
Field<Type>& ff,
const labelList& mask
) const;
// Operations for adding a field for partially covered faces
//- Add to partially covered faces for master
template<class Type>
void addToPartialFacesMaster
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const;
//- Add to partially covered faces for master, only for marked faces
template<class Type>
void maskedAddToPartialFacesMaster
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask
) const;
//- Add to partially covered faces for slave
template<class Type>
void addToPartialFacesSlave
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const;
//- Add to partially covered faces for slave, only for marked faces
template<class Type>
void maskedAddToPartialFacesSlave
(
const Field<Type>& fieldToAdd,
Field<Type>& ff,
const labelList& mask
) const;

View file

@ -294,25 +294,26 @@ void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
- boundaryMesh()[shadowID].faceCentres()
);
// Get face centres on master side
const vectorField::subField cf = faceCentres();
if (bridgeOverlap_)
{
// Get necessary mesh data from polyPatch on this (master) side
const vectorField::subField cf = faceCentres();
const vectorField::subField Sf = faceAreas();
// Get face cell centres on master side
const vectorField ccf = faceCellCentres();
const vectorField nf = Sf/mag(Sf);
// Deltas for fully uncovered faces
const vectorField uncoveredDeltas(2.0*(cf - ccf));
const vectorField uncoveredDeltas(cf - ccf);
// Scale partially overlapping faces and set uncovered deltas to
// fully uncovered faces
scaleForPartialCoverage(uncoveredDeltas, tdf());
// Set uncovered deltas to fully uncovered faces
setUncoveredFaces(uncoveredDeltas, tdf());
// Scale partially overlapping faces
scalePartialFaces(tdf());
}
// Calculate the reconstructed cell centres
reconFaceCellCentresPtr_ = new vectorField(tdf() + faceCentres());
reconFaceCellCentresPtr_ = new vectorField(tdf() + cf);
}
else
{

View file

@ -356,7 +356,8 @@ public:
const ggiZoneInterpolation& patchToPatch() const;
// Interpolation functions
// Interpolation functions and bridging operations for fully uncovered
// and partially covered faces
//- Expand face field to full zone size, including reduction
template<class Type>
@ -370,21 +371,31 @@ public:
template<class Type>
tmp<Field<Type> > interpolate(const tmp<Field<Type> >& tpf) const;
//- Bridge interpolated face field for uncovered faces
//- Set given field for fully uncovered faces
template<class Type>
void bridge
void setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Scale field for partially covered faces and set
// uncoveredFaceField to fully uncovered faces. Needed for correct
// calculation of fvPatchData: delta, deltaCoeffs and weights
//- Set given field for partially covered faces
template<class Type>
void scaleForPartialCoverage
void setPartialFaces
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;
//- Scale field for partially covered faces
template<class Type>
void scalePartialFaces(Field<Type>& ff) const;
//- Add given field to partially covered faces
template<class Type>
void addToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const;

View file

@ -201,9 +201,9 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::interpolate
template<class Type>
void Foam::ggiPolyPatch::bridge
void Foam::ggiPolyPatch::setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
@ -212,12 +212,12 @@ void Foam::ggiPolyPatch::bridge
{
FatalErrorIn
(
"template<class Type> void ggiPolyPatch::bridge\n"
"template<class Type> void ggiPolyPatch::setUncoveredFaces\n"
"(\n"
" const Field<Type>& ff,\n"
" const Field<Type>& fieldToSet,\n"
" Field<Type>& ff\n"
") const"
) << "Incorrect patch field size for bridge. Field size: "
) << "Incorrect patch field size for setting. Field size: "
<< ff.size() << " patch size: " << size()
<< abort(FatalError);
}
@ -234,11 +234,11 @@ void Foam::ggiPolyPatch::bridge
{
if (master())
{
patchToPatch().bridgeMaster(bridgeField, ff);
patchToPatch().setUncoveredFacesMaster(fieldToSet, ff);
}
else
{
patchToPatch().bridgeSlave(bridgeField, ff);
patchToPatch().setUncoveredFacesSlave(fieldToSet, ff);
}
}
else
@ -246,18 +246,18 @@ void Foam::ggiPolyPatch::bridge
// Note: since bridging is only a local operation
if (master())
{
patchToPatch().maskedBridgeMaster
patchToPatch().maskedSetUncoveredFacesMaster
(
bridgeField,
fieldToSet,
ff,
zoneAddressing()
);
}
else
{
patchToPatch().maskedBridgeSlave
patchToPatch().maskedSetUncoveredFacesSlave
(
bridgeField,
fieldToSet,
ff,
zoneAddressing()
);
@ -268,9 +268,9 @@ void Foam::ggiPolyPatch::bridge
template<class Type>
void Foam::ggiPolyPatch::scaleForPartialCoverage
void Foam::ggiPolyPatch::setPartialFaces
(
const Field<Type>& uncoveredFaceField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
@ -279,14 +279,13 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
{
FatalErrorIn
(
"template<class Type> ggiPolyPatch::scaleForPartialCoverage\n"
"template<class Type> void ggiPolyPatch::setPartialFaces\n"
"(\n"
" const Field<Type>& uncoveredFaceField,\n"
" Field<Type>& ff,\n"
" const Field<Type>& fieldToSet,\n"
" Field<Type>& ff\n"
") const"
) << "Incorrect patch field size for scaling. Field size: "
<< ff.size() << " uncovered field size: "
<< uncoveredFaceField.size() << " patch size: " << size()
) << "Incorrect patch field size for setting. Field size: "
<< ff.size() << " patch size: " << size()
<< abort(FatalError);
}
@ -302,11 +301,11 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
{
if (master())
{
patchToPatch().scaleMaster(uncoveredFaceField, ff);
patchToPatch().setPartialFacesMaster(fieldToSet, ff);
}
else
{
patchToPatch().scaleSlave(uncoveredFaceField, ff);
patchToPatch().setPartialFacesSlave(fieldToSet, ff);
}
}
else
@ -314,18 +313,146 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
// Note: since bridging is only a local operation
if (master())
{
patchToPatch().maskedScaleMaster
patchToPatch().maskedSetPartialFacesMaster
(
uncoveredFaceField,
fieldToSet,
ff,
zoneAddressing()
);
}
else
{
patchToPatch().maskedScaleSlave
patchToPatch().maskedSetPartialFacesSlave
(
uncoveredFaceField,
fieldToSet,
ff,
zoneAddressing()
);
}
}
}
}
template<class Type>
void Foam::ggiPolyPatch::scalePartialFaces(Field<Type>& ff) const
{
// Check and expand the field from patch size to zone size
if (ff.size() != size())
{
FatalErrorIn
(
"template<class Type> ggiPolyPatch::scalePartialFaces\n"
"(\n"
" Field<Type>& ff,\n"
") const"
) << "Incorrect patch field size for scaling. Field size: "
<< ff.size() << " patch size: " << size()
<< abort(FatalError);
}
if (bridgeOverlap())
{
if (empty())
{
// Patch empty, no bridging
return;
}
if (localParallel())
{
if (master())
{
patchToPatch().scalePartialMaster(ff);
}
else
{
patchToPatch().scalePartialSlave(ff);
}
}
else
{
// Note: since bridging is only a local operation
if (master())
{
patchToPatch().maskedScalePartialMaster
(
ff,
zoneAddressing()
);
}
else
{
patchToPatch().maskedScalePartialSlave
(
ff,
zoneAddressing()
);
}
}
}
}
template<class Type>
void Foam::ggiPolyPatch::addToPartialFaces
(
const Field<Type>& fieldToAdd,
Field<Type>& ff
) const
{
// Check and expand the field from patch size to zone size
if (ff.size() != size())
{
FatalErrorIn
(
"template<class Type> ggiPolyPatch::addToPartialFaces\n"
"(\n"
" const Field<Type>& fieldToAdd,\n"
" Field<Type>& ff,\n"
") const"
) << "Incorrect patch field size for adding. Field size: "
<< ff.size() << " field to add size: "
<< fieldToAdd.size() << " patch size: " << size()
<< abort(FatalError);
}
if (bridgeOverlap())
{
if (empty())
{
// Patch empty, no bridging
return;
}
if (localParallel())
{
if (master())
{
patchToPatch().addToPartialFacesMaster(fieldToAdd, ff);
}
else
{
patchToPatch().addToPartialFacesSlave(fieldToAdd, ff);
}
}
else
{
// Note: since bridging is only a local operation
if (master())
{
patchToPatch().maskedAddToPartialFacesMaster
(
fieldToAdd,
ff,
zoneAddressing()
);
}
else
{
patchToPatch().maskedAddToPartialFacesSlave
(
fieldToAdd,
ff,
zoneAddressing()
);

View file

@ -409,11 +409,11 @@ public:
template<class Type>
tmp<Field<Type> > interpolate(const tmp<Field<Type> >& tpf) const;
//- Bridge interpolated face field for uncovered faces
//- Set given field for fully uncovered faces
template<class Type>
void bridge
void setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const;

View file

@ -200,9 +200,9 @@ Foam::tmp<Foam::Field<Type> > Foam::regionCouplePolyPatch::interpolate
template<class Type>
void Foam::regionCouplePolyPatch::bridge
void Foam::regionCouplePolyPatch::setUncoveredFaces
(
const Field<Type>& bridgeField,
const Field<Type>& fieldToSet,
Field<Type>& ff
) const
{
@ -211,12 +211,12 @@ void Foam::regionCouplePolyPatch::bridge
{
FatalErrorIn
(
"tmp<Field<Type> > regionCouplePolyPatch::bridge\n"
"tmp<Field<Type> > regionCouplePolyPatch::setUncoveredFaces\n"
"(\n"
" const Field<Type>& ff,\n"
" const Field<Type>& fieldToSet,\n"
" Field<Type>& ff\n"
") const"
) << "Incorrect patch field size for bridge. Field size: "
) << "Incorrect patch field size for setting. Field size: "
<< ff.size() << " patch size: " << size()
<< abort(FatalError);
}
@ -233,11 +233,11 @@ void Foam::regionCouplePolyPatch::bridge
{
if (master())
{
patchToPatch().bridgeMaster(bridgeField, ff);
patchToPatch().setUncoveredFacesMaster(fieldToSet, ff);
}
else
{
patchToPatch().bridgeSlave(bridgeField, ff);
patchToPatch().setUncoveredFacesSlave(fieldToSet, ff);
}
}
else
@ -245,18 +245,18 @@ void Foam::regionCouplePolyPatch::bridge
// Note: since bridging is only a local operation
if (master())
{
patchToPatch().maskedBridgeMaster
patchToPatch().maskedSetUncoveredFacesMaster
(
bridgeField,
fieldToSet,
ff,
zoneAddressing()
);
}
else
{
patchToPatch().maskedBridgeSlave
patchToPatch().maskedSetUncoveredFacesSlave
(
bridgeField,
fieldToSet,
ff,
zoneAddressing()
);

View file

@ -10,6 +10,7 @@ field/magGrad/magGrad.C
field/div/div.C
field/randomise/randomise.C
field/interpolate/interpolate.C
field/helicity/helicity.C
basic/addSubtract/addSubtract.C
basic/scalarMult/scalarMult.C

View file

@ -0,0 +1,150 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "helicity.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace calcTypes
{
defineTypeNameAndDebug(helicity, 0);
addToRunTimeSelectionTable(calcType, helicity, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::calcTypes::helicity::helicity()
:
calcType()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::calcTypes::helicity::~helicity()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::calcTypes::helicity::init()
{
argList::validArgs.append("helicity");
argList::validArgs.append("fieldName");
}
void Foam::calcTypes::helicity::preCalc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
)
{}
void Foam::calcTypes::helicity::calc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
)
{
const word& fieldName = args.additionalArgs()[1];
IOobject fieldHeader
(
fieldName,
runTime.timeName(),
mesh,
IOobject::MUST_READ
);
// Check field exists
if (fieldHeader.headerOk())
{
bool processed = false;
writeHelicityField(fieldHeader, mesh, processed);
if (!processed)
{
FatalError
<< "Unable to process " << fieldName << nl
<< "No call to helicity for fields of type "
<< fieldHeader.headerClassName() << nl << nl
<< exit(FatalError);
}
}
else
{
Info<< " No " << fieldName << endl;
}
}
void Foam::calcTypes::helicity::writeHelicityField
(
const IOobject& header,
const fvMesh& mesh,
bool& processed
)
{
if (header.headerClassName() == volVectorField::typeName)
{
Info<< " Reading " << header.name() << endl;
volVectorField field(header, mesh);
Info<< " Calculating helicity" << header.name() << endl;
volScalarField helicityField
(
IOobject
(
"helicity" + header.name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ
),
field & (fvc::curl(field))
);
Info << "helicity(" << header.name() << "): max: "
<< gMax(helicityField.internalField())
<< " min: " << gMin(helicityField.internalField()) << endl;
helicityField.write();
processed = true;
}
}
// ************************************************************************* //

View file

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::calcTypes::helicity
Description
Calculates and writes the helicity of a vector field for each time, defines
as: U & rot(U)
SourceFiles
helicity.C
\*---------------------------------------------------------------------------*/
#ifndef helicity_H
#define helicity_H
#include "calcType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace calcTypes
{
/*---------------------------------------------------------------------------*\
Class helicity Declaration
\*---------------------------------------------------------------------------*/
class helicity
:
public calcType
{
// Private Member Functions
//- Disallow default bitwise copy construct
helicity(const helicity&);
//- Disallow default bitwise assignment
void operator=(const helicity&);
protected:
// Member Functions
// Calculation routines
//- Initialise - typically setting static variables,
// e.g. command line arguments
virtual void init();
//- Pre-time loop calculations
virtual void preCalc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
);
//- Time loop calculations
virtual void calc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
);
// I-O
//- Write component fields
void writeHelicityField
(
const IOobject& header,
const fvMesh& mesh,
bool& processed
);
public:
//- Runtime type information
TypeName("helicity");
// Constructors
//- Construct null
helicity();
// Destructor
virtual ~helicity();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace calcTypes
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -86,6 +86,7 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
if (mesh.foundObject<surfaceScalarField>(phiName_))
{
Info<< "Checking flux " << phiName_ << " GGI balance." << endl;
const surfaceScalarField& phi =
mesh.lookupObject<surfaceScalarField>(phiName_);
@ -110,24 +111,22 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
visited[patchI] = true;
visited[shadowPatchI] = true;
// Calculate local and shadow flux
scalar localArea =
// Calculate local flux and area
const scalar localArea =
gSum(phi.mesh().magSf().boundaryField()[patchI]);
scalar localFlux = gSum(phi.boundaryField()[patchI]);
scalar localFluxMag = gSumMag(phi.boundaryField()[patchI]);
const scalar localFlux = gSum(phi.boundaryField()[patchI]);
const scalar localFluxMag = mag(localFlux);
scalar shadowArea =
// Calculate shadow flux and area
const scalar shadowArea =
gSum(phi.mesh().magSf().boundaryField()[shadowPatchI]);
scalar shadowFlux =
const scalar shadowFlux =
gSum(phi.boundaryField()[shadowPatchI]);
const scalar shadowFluxMag = mag(shadowFlux);
scalar shadowFluxMag =
gSumMag(phi.boundaryField()[shadowPatchI]);
// Report
Info<< "GGI pair (" << ggiPatch.name() << ", "
<< ggiPatch.shadow().name() << ")" << nl
<< "Area: " << localArea << " " << shadowArea
@ -138,7 +137,7 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
<< "Flux: " << localFluxMag << " " << shadowFluxMag
<< " Diff = " << localFlux + shadowFlux << " or "
<< mag(localFlux + shadowFlux)/
(localFluxMag + SMALL)*100
(Foam::max(localFluxMag, shadowFluxMag) + SMALL)*100
<< " %" << endl;
}
else if
@ -157,22 +156,22 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
visited[patchI] = true;
visited[shadowPatchI] = true;
// Calculate local and shadow flux
scalar localArea =
// Calculate local flux and area
const scalar localArea =
gSum(phi.mesh().magSf().boundaryField()[patchI]);
scalar localFlux = gSum(phi.boundaryField()[patchI]);
scalar localFluxMag = gSumMag(phi.boundaryField()[patchI]);
const scalar localFlux = gSum(phi.boundaryField()[patchI]);
const scalar localFluxMag = mag(localFlux);
scalar shadowArea =
// Calculate shadow flux and area
const scalar shadowArea =
gSum(phi.mesh().magSf().boundaryField()[shadowPatchI]);
scalar shadowFlux =
const scalar shadowFlux =
gSum(phi.boundaryField()[shadowPatchI]);
const scalar shadowFluxMag = mag(shadowFlux);
scalar shadowFluxMag =
gSumMag(phi.boundaryField()[shadowPatchI]);
// Report
Info<< "Cyclic GGI pair (" << ggiPatch.name() << ", "
<< ggiPatch.shadow().name() << ")" << nl
<< "Area: " << localArea << " " << shadowArea
@ -183,7 +182,7 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
<< "Flux: " << localFluxMag << " " << shadowFluxMag
<< " Diff = " << localFlux + shadowFlux << " or "
<< mag(localFlux + shadowFlux)/
(localFluxMag + SMALL)*100
(Foam::max(localFluxMag, shadowFluxMag) + SMALL)*100
<< " %" << endl;
}
else if
@ -208,30 +207,25 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
visited[patchI] = true;
visited[shadowPatchI] = true;
// Calculate local and shadow flux
scalar localArea =
// Calculate local flux and area
const scalar localArea =
gSum(phi.mesh().magSf().boundaryField()[patchI])
*ggiPatch.nCopies();
scalar localFlux =
const scalar localFlux =
gSum(phi.boundaryField()[patchI])
*ggiPatch.nCopies();
const scalar localFluxMag = mag(localFlux);
scalar localFluxMag =
gSumMag(phi.boundaryField()[patchI])
*ggiPatch.nCopies();
scalar shadowArea =
// Calculate shadow flux and area
const scalar shadowArea =
gSum(phi.mesh().magSf().boundaryField()[shadowPatchI])
*ggiShadowPatch.nCopies();
scalar shadowFlux =
const scalar shadowFlux =
gSum(phi.boundaryField()[shadowPatchI])
*ggiShadowPatch.nCopies();
scalar shadowFluxMag =
gSumMag(phi.boundaryField()[shadowPatchI])
*ggiShadowPatch.nCopies();
const scalar shadowFluxMag = mag(shadowFlux);
Info<< "Overlap GGI pair (" << ggiPatch.name() << ", "
<< ggiPatch.shadow().name() << ")" << nl
@ -253,7 +247,7 @@ bool Foam::ggiCheckFunctionObject::execute(const bool forceWrite)
else
{
InfoIn("bool ggiCheckFunctionObject::execute(const bool forceWrite)")
<< "Cannot find flux field phi"
<< "Cannot find flux field " << phiName_
<< endl;
return false;