From 617fcb0d3259db013b67be80e65624429a1ad0a3 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Thu, 29 Sep 2011 02:12:36 +0100 Subject: [PATCH] Improvements in conjugate heat solver: non-matching interface, parallelisation --- .../polyPatches/constraint/ggi/ggiPolyPatch.C | 9 +- .../polyPatches/constraint/ggi/ggiPolyPatch.H | 188 ++-- .../regionCouple/regionCouplePolyPatch.C | 862 +++++++++++++++--- .../regionCouple/regionCouplePolyPatch.H | 307 +++++-- .../regionCouplePolyPatchTemplates.C | 350 +++++++ src/finiteVolume/Make/files | 5 +- .../constraint/ggi/ggiFvPatchField.C | 12 + .../constraint/ggi/ggiFvPatchField.H | 6 + .../regionCouple/regionCoupleFvPatchField.C | 275 ------ .../regionCouplingFvPatchField.C | 468 ++++++++++ .../regionCouplingFvPatchField.H} | 127 ++- .../regionCouplingFvPatchFields.C} | 4 +- .../regionCouplingFvPatchFields.H} | 12 +- .../regionCouplingFvPatchFieldsFwd.H} | 10 +- .../regionCouplingFvsPatchField.C} | 28 +- .../regionCouplingFvsPatchField.H} | 41 +- .../regionCouplingFvsPatchFields.C} | 4 +- .../regionCouplingFvsPatchFields.H} | 12 +- .../regionCouplingFvsPatchFieldsFwd.H} | 10 +- .../fvPatches/constraint/ggi/ggiFvPatch.C | 17 +- .../regionCouple/regionCoupleFvPatch.C | 208 ++++- .../regionCouple/regionCoupleFvPatch.H | 124 ++- .../schemes/harmonic/harmonic.H | 71 +- .../schemes/harmonic/magLongDelta.C | 246 +++++ .../schemes/harmonic/magLongDelta.H | 124 +++ 25 files changed, 2765 insertions(+), 755 deletions(-) create mode 100644 src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatchTemplates.C delete mode 100644 src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C create mode 100644 src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.C rename src/finiteVolume/fields/fvPatchFields/constraint/{regionCouple/regionCoupleFvPatchField.H => regionCoupling/regionCouplingFvPatchField.H} (60%) rename src/finiteVolume/fields/fvPatchFields/constraint/{regionCouple/regionCoupleFvPatchFields.C => regionCoupling/regionCouplingFvPatchFields.C} (95%) rename src/finiteVolume/fields/{fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.H => fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.H} (88%) rename src/finiteVolume/fields/fvPatchFields/constraint/{regionCouple/regionCoupleFvPatchFieldsFwd.H => regionCoupling/regionCouplingFvPatchFieldsFwd.H} (89%) rename src/finiteVolume/fields/fvsPatchFields/constraint/{regionCouple/regionCoupleFvsPatchField.C => regionCoupling/regionCouplingFvsPatchField.C} (84%) rename src/finiteVolume/fields/fvsPatchFields/constraint/{regionCouple/regionCoupleFvsPatchField.H => regionCoupling/regionCouplingFvsPatchField.H} (80%) rename src/finiteVolume/fields/fvsPatchFields/constraint/{regionCouple/regionCoupleFvsPatchFields.C => regionCoupling/regionCouplingFvsPatchFields.C} (95%) rename src/finiteVolume/fields/{fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.H => fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.H} (88%) rename src/finiteVolume/fields/fvsPatchFields/constraint/{regionCouple/regionCoupleFvsPatchFieldsFwd.H => regionCoupling/regionCouplingFvsPatchFieldsFwd.H} (89%) create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.C create mode 100644 src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.H diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index 102d462af..87e20aabb 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -72,12 +72,6 @@ void Foam::ggiPolyPatch::calcZoneAddressing() const << abort(FatalError); } - if (debug) - { - Pout<< "ggiPolyPatch::calcZoneAddressing() const for patch " - << index() << endl; - } - // Calculate patch-to-zone addressing zoneAddressingPtr_ = new labelList(size()); labelList& zAddr = *zoneAddressingPtr_; @@ -195,7 +189,7 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const new ggiZoneInterpolation ( zone()(), // This zone reference - shadow().zone()(), // This shadow zone reference + shadow().zone()(), // Shadow zone reference forwardT(), reverseT(), shadow().separation(), // Slave-to-master separation. Bug fix @@ -640,6 +634,7 @@ Foam::label Foam::ggiPolyPatch::zoneIndex() const { FatalErrorIn("label ggiPolyPatch::zoneIndex() const") << "Face zone name " << zoneName_ + << " for GGI patch " << name() << " not found. Please check your GGI interface definition." << abort(FatalError); } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H index 85edd7abf..fe40633e7 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H @@ -279,116 +279,128 @@ public: // Member functions - //- Return shadow patch name - const word& shadowName() const - { - return shadowName_; - } + // Basic info - //- Return name of interpolation face zone - const word& zoneName() const - { - return zoneName_; - } + //- Return shadow patch name + const word& shadowName() const + { + return shadowName_; + } - //- Return shadow patch index - label shadowIndex() const; + //- Return name of interpolation face zone + const word& zoneName() const + { + return zoneName_; + } - //- Return zone patch index - label zoneIndex() const; + //- Return shadow patch index + label shadowIndex() const; - //- Return shadow patch - const ggiPolyPatch& shadow() const; + //- Return zone patch index + label zoneIndex() const; - //- Return interpolation face zone - const faceZone& zone() const; + //- Return shadow patch + const ggiPolyPatch& shadow() const; - //- Is this the master side? - bool master() const - { - return index() < shadowIndex(); - } + //- Return interpolation face zone + const faceZone& zone() const; - //- Is this the slave side? - bool slave() const - { - return !master(); - } + // Interpolation data - //- Use bridging to fix overlap error in interpolation - bool bridgeOverlap() const - { - return bridgeOverlap_; - } + //- Is this the master side? + bool master() const + { + return index() < shadowIndex(); + } - //- Return patch-to-zone addressing - const labelList& zoneAddressing() const; + //- Is this the slave side? + bool slave() const + { + return !master(); + } - //- Return remote patch-to-zone addressing - const labelList& remoteZoneAddressing() const; + //- Use bridging to fix overlap error in interpolation + bool bridgeOverlap() const + { + return bridgeOverlap_; + } - //- Is the patch localised on a single processor - bool localParallel() const; + //- Return patch-to-zone addressing + const labelList& zoneAddressing() const; - //- Return reference to patch-to-patch interpolation - // Used only for addressing - const ggiZoneInterpolation& patchToPatch() const; + //- Return remote patch-to-zone addressing + const labelList& remoteZoneAddressing() const; - //- Expand face field to full zone size, including reduction - // New. HJ, 12/Jun/2011 - template - tmp > fastExpand(const Field& pf) const; + //- Is the patch localised on a single processor + bool localParallel() const; - //- Expand face field to full zone size, including reduction - // Obsolete. HJ, 12/Jun/2011 - template - tmp > expand(const Field& pf) const; - - //- Filter zone field to patch size - // Obsolete. HJ, 12/Jun/2011 - template - tmp > filter(const Field& zf) const; + //- Return reference to patch-to-patch interpolation + // Used only for addressing + const ggiZoneInterpolation& patchToPatch() const; - //- Interpolate face field: given field on a the shadow side, - // create an interpolated field on this side - template - tmp > interpolate(const Field& pf) const; + // Interpolation functions - template - tmp > interpolate(const tmp >& tpf) const; + //- Expand face field to full zone size, including reduction + // New. HJ, 12/Jun/2011 + template + tmp > fastExpand(const Field& pf) const; - //- Bridge interpolated face field for uncovered faces - template - void bridge - ( - const Field& bridgeField, - Field& ff - ) const; + //- Expand face field to full zone size, including reduction + // Obsolete. HJ, 12/Jun/2011 + template + tmp > expand(const Field& pf) const; - //- Return reconstructed cell centres - const vectorField& reconFaceCellCentres() const; + //- Filter zone field to patch size + // Obsolete. HJ, 12/Jun/2011 + template + tmp > filter(const Field& zf) const; - //- Initialize ordering for primitivePatch. Does not - // refer to *this (except for name() and type() etc.) - virtual void initOrder(const primitivePatch&) const; + //- Interpolate face field: given field on a the shadow side, + // create an interpolated field on this side + template + tmp > interpolate(const Field& pf) const; - //- Return new ordering for primitivePatch. - // Ordering is -faceMap: for every face - // index of the new face -rotation: for every new face the clockwise - // shift of the original face. Return false if nothing changes - // (faceMap is identity, rotation is 0), true otherwise. - virtual bool order - ( - const primitivePatch&, - labelList& faceMap, - labelList& rotation - ) const; + template + tmp > interpolate(const tmp >& tpf) const; - //- Synchronise communications of ordering for primitivePatch - // Used in cases when no topological change happens locally, - // but is happening on other processors - virtual void syncOrder() const; + //- Bridge interpolated face field for uncovered faces + template + void bridge + ( + const Field& bridgeField, + Field& ff + ) const; + + + // Geometric data + + //- Return reconstructed cell centres + const vectorField& reconFaceCellCentres() const; + + + // Patch ordering + + //- Initialize ordering for primitivePatch. Does not + // refer to *this (except for name() and type() etc.) + virtual void initOrder(const primitivePatch&) const; + + //- Return new ordering for primitivePatch. + // Ordering is -faceMap: for every face + // index of the new face -rotation: for every new face the + // clockwise shift of the original face. Return false if nothing + // changes (faceMap is identity, rotation is 0), true otherwise. + virtual bool order + ( + const primitivePatch&, + labelList& faceMap, + labelList& rotation + ) const; + + //- Synchronise communications of ordering for primitivePatch + // Used in cases when no topological change happens locally, + // but is happening on other processors + virtual void syncOrder() const; //- Write diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C index 9dd5b0842..4883d2905 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C @@ -34,6 +34,7 @@ Author #include "Time.H" #include "polyMesh.H" #include "polyPatchID.H" +#include "ZoneIDs.H" #include "objectRegistry.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -42,56 +43,415 @@ namespace Foam { defineTypeNameAndDebug(regionCouplePolyPatch, 0); + addToRunTimeSelectionTable(polyPatch, regionCouplePolyPatch, word); addToRunTimeSelectionTable(polyPatch, regionCouplePolyPatch, dictionary); } // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::regionCouplePolyPatch::calcInterpolation() const +bool Foam::regionCouplePolyPatch::active() const +{ + // Try to find face zone + faceZoneID zone(zoneName_, boundaryMesh().mesh().faceZones()); + + if (!zone.active()) + { + return false; + } + + // Try to find shadow region + if + ( + boundaryMesh().mesh().db().parent().foundObject + ( + shadowRegionName_ + ) + ) + { + // Shadow region present + const polyMesh& sr = boundaryMesh().mesh().db().parent(). + objectRegistry::lookupObject + ( + shadowRegionName_ + ); + + polyPatchID shadowPatch(shadowPatchName_, sr.boundaryMesh()); + + // If shadow patch is active, all components are ready + return shadowPatch.active(); + } + else + { + // No shadow region + return false; + } +} + + +void Foam::regionCouplePolyPatch::calcZoneAddressing() const +{ + // Calculate patch-to-zone addressing + if (zoneAddressingPtr_) + { + FatalErrorIn("void regionCouplePolyPatch::calcZoneAddressing() const") + << "Patch to zone addressing already calculated" + << abort(FatalError); + } + + if (debug) + { + Pout<< "regionCouplePolyPatch::calcZoneAddressing() const for patch " + << index() << endl; + } + + // Calculate patch-to-zone addressing + zoneAddressingPtr_ = new labelList(size()); + labelList& zAddr = *zoneAddressingPtr_; + const faceZone& myZone = zone(); + + for (label i = 0; i < size(); i++) + { + zAddr[i] = myZone.whichFace(start() + i); + } + + // Check zone addressing + if (zAddr.size() > 0 && min(zAddr) < 0) + { + FatalErrorIn("void regionCouplePolyPatch::calcZoneAddressing() const") + << "Problem with patch-to zone addressing: some patch faces " + << "not found in interpolation zone" + << abort(FatalError); + } +} + + +void Foam::regionCouplePolyPatch::calcRemoteZoneAddressing() const +{ + // Calculate patch-to-zone addressing + if (remoteZoneAddressingPtr_) + { + FatalErrorIn + ( + "void regionCouplePolyPatch::calcRemoteZoneAddressing() const" + ) << "Patch to remote zone addressing already calculated" + << abort(FatalError); + } + + // Once zone addressing is established, visit the opposite side and find + // out which face data is needed for interpolation + boolList usedShadows(shadow().zone().size(), false); + + const labelList& zAddr = zoneAddressing(); + + if (master()) + { + const labelListList& addr = patchToPatch().masterAddr(); + + forAll (zAddr, mfI) + { + const labelList& nbrs = addr[zAddr[mfI]]; + + forAll (nbrs, nbrI) + { + usedShadows[nbrs[nbrI]] = true; + } + } + } + else + { + const labelListList& addr = patchToPatch().slaveAddr(); + + forAll (zAddr, mfI) + { + const labelList& nbrs = addr[zAddr[mfI]]; + + forAll (nbrs, nbrI) + { + usedShadows[nbrs[nbrI]] = true; + } + } + } + + // Count and pick up shadow indices + label nShadows = 0; + + forAll (usedShadows, sI) + { + if (usedShadows[sI]) + { + nShadows++; + } + } + + remoteZoneAddressingPtr_ = new labelList(nShadows); + labelList& rza = *remoteZoneAddressingPtr_; + + // Reset counter for re-use + nShadows = 0; + + forAll (usedShadows, sI) + { + if (usedShadows[sI]) + { + rza[nShadows] = sI; + nShadows++; + } + } +} + + +void Foam::regionCouplePolyPatch::calcPatchToPatch() const { // Create patch-to-patch interpolation if (patchToPatchPtr_) { - FatalErrorIn("void regionCouplePolyPatch::calcInterpolation() const") + FatalErrorIn("void regionCouplePolyPatch::calcPatchToPatch() const") << "Patch to patch interpolation already calculated" << abort(FatalError); } - // Get shadow region - const polyMesh& sr = shadowRegion(); - - // Grab shadow patch index - polyPatchID shadow(shadowPatchName_, sr.boundaryMesh()); - - if (!shadow.active()) + if (master()) { - FatalErrorIn("void regionCouplePolyPatch::calcInterpolation() const") - << "Shadow patch name " << shadowPatchName_ - << " not found. Please check your regionCouple interface." - << abort(FatalError); - } + // Create interpolation for zones + patchToPatchPtr_ = + new ggiZoneInterpolation + ( + zone()(), // This zone reference + shadow().zone()(), // Shadow zone reference + forwardT(), + reverseT(), + shadow().separation(), // Slave-to-master separation. Bug fix + 0, // Non-overlapping face tolerances + 0, + true, // Rescale weighting factors +// ggiInterpolation::AABB + ggiInterpolation::BB_OCTREE // Octree search, MB. + ); - shadowIndex_ = shadow.index(); - - // Check the other side is a regionCouple - if (!isType(sr.boundaryMesh()[shadowIndex_])) - { - FatalErrorIn("void regionCouplePolyPatch::calcInterpolation() const") - << "Shadow of regionCouple patch " << name() - << " named " << shadowPatchName() - << " is not a regionCouple." << nl - << "This is not allowed. Please check your mesh definition." - << abort(FatalError); - } - - patchToPatchPtr_ = - new patchToPatchInterpolation + // Abort immediately if uncovered faces are present and the option + // bridgeOverlap is not set. + if ( - sr.boundaryMesh()[shadowIndex_], - *this, - intersection::VISIBLE - ); + ( + patchToPatch().uncoveredMasterFaces().size() > 0 + && !bridgeOverlap() + ) + || ( + patchToPatch().uncoveredSlaveFaces().size() > 0 + && !shadow().bridgeOverlap() + ) + ) + { + FatalErrorIn + ( + "void regionCouplePolyPatch::calcPatchToPatch() const" + ) << "Found uncovered faces for GGI interface " + << name() << "/" << shadowPatchName() + << " while the bridgeOverlap option is not set " + << "in the boundary file." << endl + << "This is an unrecoverable error. Aborting." + << abort(FatalError); + } + } + else + { + FatalErrorIn("void regionCouplePolyPatch::calcPatchToPatch() const") + << "Attempting to create GGIInterpolation on a shadow" + << abort(FatalError); + } +} + + +void Foam::regionCouplePolyPatch::calcReconFaceCellCentres() const +{ + if (reconFaceCellCentresPtr_) + { + FatalErrorIn + ( + "void regionCouplePolyPatch::calcReconFaceCellCentres() const" + ) << "Reconstructed cell centres already calculated" + << abort(FatalError); + } + + // Create neighbouring face centres using interpolation + if (master()) + { + const label shadowID = shadowIndex(); + + const polyMesh& sr = shadowRegion(); + + // Get the transformed and interpolated shadow face cell centers + reconFaceCellCentresPtr_ = + new vectorField + ( + interpolate + ( + sr.boundaryMesh()[shadowID].faceCellCentres() + - sr.boundaryMesh()[shadowID].faceCentres() + ) + + faceCentres() + ); + } + else + { + FatalErrorIn + ( + "void regionCouplePolyPatch::calcReconFaceCellCentres() const" + ) << "Attempting to create reconFaceCellCentres on a shadow" + << abort(FatalError); + } +} + + +void Foam::regionCouplePolyPatch::calcLocalParallel() const +{ + // Calculate patch-to-zone addressing + if (localParallelPtr_) + { + FatalErrorIn("void regionCouplePolyPatch::calcLocalParallel() const") + << "Local parallel switch already calculated" + << abort(FatalError); + } + + localParallelPtr_ = new bool(false); + bool& emptyOrComplete = *localParallelPtr_; + + // If running in serial, all GGIs are expanded to zone size. + // This happens on decomposition and reconstruction where + // size and shadow size may be zero, but zone size may not + // HJ, 1/Jun/2011 + if (!Pstream::parRun()) + { + emptyOrComplete = false; + } + else + { + // Calculate localisation on master and shadow + emptyOrComplete = + ( + zone().size() == size() + && shadow().zone().size() == shadow().size() + ) + || (size() == 0 && shadow().size() == 0); + + reduce(emptyOrComplete, andOp()); + } + + if (debug && Pstream::parRun()) + { + Info<< "regionCouple patch Master: " << name() + << " Slave: " << shadowPatchName() << " is "; + + if (emptyOrComplete) + { + Info<< "local parallel" << endl; + } + else + { + Info<< "split between multiple processors" << endl; + } + } +} + + +void Foam::regionCouplePolyPatch::calcSendReceive() const +{ + // Note: all processors will execute calcSendReceive but only master will + // hold the information. Therefore, pointers on slave processors + // will remain meaningless, but for purposes of consistency + // (of the calc-call) they will be set to zero-sized array + // HJ, 4/Jun/2011 + + if (receiveAddrPtr_ || sendAddrPtr_) + { + FatalErrorIn("void regionCouplePolyPatch::calcSendReceive() const") + << "Send-receive addressing already calculated" + << abort(FatalError); + } + + if (debug) + { + Pout<< "regionCouplePolyPatch::calcSendReceive() const for patch " + << index() << endl; + } + + if (!Pstream::parRun()) + { + FatalErrorIn("void regionCouplePolyPatch::calcSendReceive() const") + << "Requested calculation of send-receive addressing for a " + << "serial run. This is not allowed" + << abort(FatalError); + } + + // Master will receive and store the maps + if (Pstream::master()) + { + receiveAddrPtr_ = new labelListList(Pstream::nProcs()); + labelListList& rAddr = *receiveAddrPtr_; + + sendAddrPtr_ = new labelListList(Pstream::nProcs()); + labelListList& sAddr = *sendAddrPtr_; + + // Insert master + rAddr[0] = zoneAddressing(); + + for (label procI = 1; procI < Pstream::nProcs(); procI++) + { + // Note: must use normal comms because the size of the + // communicated lists is unknown on the receiving side + // HJ, 4/Jun/2011 + + // Opt: reconsider mode of communication + IPstream ip(Pstream::scheduled, procI); + + rAddr[procI] = labelList(ip); + + sAddr[procI] = labelList(ip); + } + } + else + { + // Create dummy pointers: only master processor stores maps + receiveAddrPtr_ = new labelListList(); + sendAddrPtr_ = new labelListList(); + + // Send information to master + const labelList& za = zoneAddressing(); + const labelList& ra = remoteZoneAddressing(); + + // Note: must use normal comms because the size of the + // communicated lists is unknown on the receiving side + // HJ, 4/Jun/2011 + + // Opt: reconsider mode of communication + OPstream op(Pstream::scheduled, Pstream::masterNo()); + + // Send local and remote addressing to master + op << za << ra; + } +} + + +const Foam::labelListList& Foam::regionCouplePolyPatch::receiveAddr() const +{ + if (!receiveAddrPtr_) + { + calcSendReceive(); + } + + return *receiveAddrPtr_; +} + + +const Foam::labelListList& Foam::regionCouplePolyPatch::sendAddr() const +{ + if (!sendAddrPtr_) + { + calcSendReceive(); + } + + return *sendAddrPtr_; } @@ -107,70 +467,62 @@ Foam::regionCouplePolyPatch& Foam::regionCouplePolyPatch::shadow() } -const Foam::patchToPatchInterpolation& -Foam::regionCouplePolyPatch::patchToPatch() const -{ - if (!attached_) - { - FatalErrorIn - ( - "const patchToPatchInterpolation& " - "regionCouplePolyPatch::patchToPatch() const" - ) << "Requesting patchToPatchInterpolation in detached state" - << abort(FatalError); - } - - if (!patchToPatchPtr_) - { - calcInterpolation(); - } - - return *patchToPatchPtr_; -} - - -void Foam::regionCouplePolyPatch::calcReconFaceCellCentres() const -{ - // Create neighbouring face centres using interpolation - - // Get shadow region - const polyMesh& sr = shadowRegion(); - - const label shadowID = shadowIndex(); - - // Reconstruct the shadow cell face centres - vectorField localCtrs = faceCellCentres(); - vectorField reconCtrs = - patchToPatch().faceInterpolate - ( - sr.boundaryMesh()[shadowID].faceCellCentres() - ); - - // Calculate reconstructed centres by eliminating non-orthogonality - const vectorField& n = faceNormals(); - - reconFaceCellCentresPtr_ = - new vectorField(localCtrs + n*(n & (reconCtrs - localCtrs))); -} - - void Foam::regionCouplePolyPatch::clearGeom() const { deleteDemandDrivenData(reconFaceCellCentresPtr_); + + // Remote addressing and send-receive maps depend on the local + // position. Therefore, it needs to be recalculated at mesh motion. + // Local zone addressing does not change with mesh motion + // HJ, 23/Jun/2011 + deleteDemandDrivenData(remoteZoneAddressingPtr_); + + deleteDemandDrivenData(receiveAddrPtr_); + deleteDemandDrivenData(sendAddrPtr_); } -void Foam::regionCouplePolyPatch::clearOut() const +void Foam::regionCouplePolyPatch::clearOut() { clearGeom(); + deleteDemandDrivenData(zoneAddressingPtr_); deleteDemandDrivenData(patchToPatchPtr_); + deleteDemandDrivenData(localParallelPtr_); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -// Construct from components +Foam::regionCouplePolyPatch::regionCouplePolyPatch +( + const word& name, + const label size, + const label start, + const label index, + const polyBoundaryMesh& bm +) +: + coupledPolyPatch(name, size, start, index, bm), + shadowRegionName_(word::null), + shadowPatchName_(word::null), + zoneName_(word::null), + attached_(false), + master_(false), + isWall_(false), + bridgeOverlap_(false), + shadowIndex_(-1), + zoneIndex_(-1), + patchToPatchPtr_(NULL), + zoneAddressingPtr_(NULL), + remoteZoneAddressingPtr_(NULL), + reconFaceCellCentresPtr_(NULL), + localParallelPtr_(NULL), + receiveAddrPtr_(NULL), + sendAddrPtr_(NULL) +{} + + Foam::regionCouplePolyPatch::regionCouplePolyPatch ( const word& name, @@ -180,22 +532,33 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch const polyBoundaryMesh& bm, const word& shadowRegionName, const word& shadowPatchName, + const word& zoneName, const bool attached, - const bool isWall + const bool master, + const bool isWall, + const bool bridgeOverlap ) : coupledPolyPatch(name, size, start, index, bm), shadowRegionName_(shadowRegionName), shadowPatchName_(shadowPatchName), + zoneName_(zoneName), attached_(attached), + master_(master), isWall_(isWall), + bridgeOverlap_(bridgeOverlap), shadowIndex_(-1), + zoneIndex_(-1), patchToPatchPtr_(NULL), - reconFaceCellCentresPtr_(NULL) + zoneAddressingPtr_(NULL), + remoteZoneAddressingPtr_(NULL), + reconFaceCellCentresPtr_(NULL), + localParallelPtr_(NULL), + receiveAddrPtr_(NULL), + sendAddrPtr_(NULL) {} -// Construct from dictionary Foam::regionCouplePolyPatch::regionCouplePolyPatch ( const word& name, @@ -207,15 +570,23 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch coupledPolyPatch(name, dict, index, bm), shadowRegionName_(dict.lookup("shadowRegion")), shadowPatchName_(dict.lookup("shadowPatch")), + zoneName_(dict.lookup("zone")), attached_(dict.lookup("attached")), + master_(dict.lookup("master")), isWall_(dict.lookup("isWall")), + bridgeOverlap_(dict.lookup("bridgeOverlap")), shadowIndex_(-1), + zoneIndex_(-1), patchToPatchPtr_(NULL), - reconFaceCellCentresPtr_(NULL) + zoneAddressingPtr_(NULL), + remoteZoneAddressingPtr_(NULL), + reconFaceCellCentresPtr_(NULL), + localParallelPtr_(NULL), + receiveAddrPtr_(NULL), + sendAddrPtr_(NULL) {} -//- Construct as copy, resetting the boundary mesh Foam::regionCouplePolyPatch::regionCouplePolyPatch ( const regionCouplePolyPatch& pp, @@ -225,15 +596,22 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch coupledPolyPatch(pp, bm), shadowRegionName_(pp.shadowRegionName_), shadowPatchName_(pp.shadowPatchName_), + zoneName_(pp.zoneName_), attached_(pp.attached_), + master_(pp.master_), isWall_(pp.isWall_), shadowIndex_(-1), + zoneIndex_(-1), patchToPatchPtr_(NULL), - reconFaceCellCentresPtr_(NULL) + zoneAddressingPtr_(NULL), + remoteZoneAddressingPtr_(NULL), + reconFaceCellCentresPtr_(NULL), + localParallelPtr_(NULL), + receiveAddrPtr_(NULL), + sendAddrPtr_(NULL) {} -//- Construct as copy, resetting the face list and boundary mesh data Foam::regionCouplePolyPatch::regionCouplePolyPatch ( const regionCouplePolyPatch& pp, @@ -246,11 +624,19 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch coupledPolyPatch(pp, bm, index, newSize, newStart), shadowRegionName_(pp.shadowRegionName_), shadowPatchName_(pp.shadowPatchName_), + zoneName_(pp.zoneName_), attached_(pp.attached_), + master_(pp.master_), isWall_(pp.isWall_), shadowIndex_(-1), + zoneIndex_(-1), patchToPatchPtr_(NULL), - reconFaceCellCentresPtr_(NULL) + zoneAddressingPtr_(NULL), + remoteZoneAddressingPtr_(NULL), + reconFaceCellCentresPtr_(NULL), + localParallelPtr_(NULL), + receiveAddrPtr_(NULL), + sendAddrPtr_(NULL) {} @@ -272,11 +658,134 @@ bool Foam::regionCouplePolyPatch::coupled() const const Foam::polyMesh& Foam::regionCouplePolyPatch::shadowRegion() const { - return boundaryMesh().mesh().db().parent(). - objectRegistry::lookupObject + if (shadowRegionName_ != Foam::word::null) + { + return boundaryMesh().mesh().db().parent(). + objectRegistry::lookupObject + ( + shadowRegionName() + ); + } + else + { + FatalErrorIn ( - shadowRegionName() - ); + "const polyMesh& regionCouplePolyPatch::shadowRegion() const" + ) << "Requested shadowRegion which is not available" + << abort(FatalError); + + // Dummy return + return boundaryMesh().mesh(); + } +} + + +Foam::label Foam::regionCouplePolyPatch::shadowIndex() const +{ + if + ( + shadowIndex_ == -1 + && shadowRegionName_ != Foam::word::null + && shadowPatchName_ != Foam::word::null + ) + { + // Grab shadow patch index from shadow region + const polyMesh& sr = shadowRegion(); + + polyPatchID shadow(shadowPatchName_, sr.boundaryMesh()); + + if (!shadow.active()) + { + FatalErrorIn("label regionCouplePolyPatch::shadowIndex() const") + << "Shadow patch name " << shadowPatchName_ + << " not found. Please check your region couple " + << "interface definition." + << abort(FatalError); + } + + shadowIndex_ = shadow.index(); + + // Check the other side is a region couple + if (!isA(sr.boundaryMesh()[shadowIndex_])) + { + FatalErrorIn("label regionCouplePolyPatch::shadowIndex() const") + << "Shadow of region couple patch " << name() + << " named " << shadowPatchName() + << " on region " << shadowRegionName() + << " is not a region couple. Type: " + << boundaryMesh()[shadowIndex_].type() << nl + << "This is not allowed. Please check your mesh definition." + << abort(FatalError); + } + + // Check for region couple onto self + if (index() == shadowIndex_ && &sr == &boundaryMesh().mesh()) + { + FatalErrorIn("label regionCouplePolyPatch::shadowIndex() const") + << "region couple patch " << name() + << " created as its own shadow" + << abort(FatalError); + } + + // Check definition of master and slave side + const regionCouplePolyPatch& sp = + refCast + ( + sr.boundaryMesh()[shadowIndex_] + ); + + if (master() == sp.master()) + { + FatalErrorIn("label regionCouplePolyPatch::shadowIndex() const") + << "Region couple patch " << name() + << " and its shadow " << shadowPatchName() + << " on region " << shadowRegionName() + << ". Clash on master-slave definition." << nl + << "This is not allowed. Please check your mesh definition." + << abort(FatalError); + } + } + + return shadowIndex_; +} + + +Foam::label Foam::regionCouplePolyPatch::zoneIndex() const +{ + if (zoneIndex_ == -1 && zoneName_ != Foam::word::null) + { + // Grab zone patch index + faceZoneID zone(zoneName_, boundaryMesh().mesh().faceZones()); + + if (!zone.active()) + { + FatalErrorIn("label regionCouplePolyPatch::zoneIndex() const") + << "Face zone name " << zoneName_ + << " for region couple patch " << name() + << " not found. Please check your region couple " + << "interface definition." + << abort(FatalError); + } + + zoneIndex_ = zone.index(); + } + + return zoneIndex_; +} + + +const Foam::regionCouplePolyPatch& Foam::regionCouplePolyPatch::shadow() const +{ + return refCast + ( + shadowRegion().boundaryMesh()[shadowIndex()] + ); +} + + +const Foam::faceZone& Foam::regionCouplePolyPatch::zone() const +{ + return boundaryMesh().mesh().faceZones()[zoneIndex()]; } @@ -310,23 +819,73 @@ void Foam::regionCouplePolyPatch::detach() const } -Foam::label Foam::regionCouplePolyPatch::shadowIndex() const +const Foam::labelList& Foam::regionCouplePolyPatch::zoneAddressing() const { - if (shadowIndex_ == -1) + if (!zoneAddressingPtr_) { - calcInterpolation(); + calcZoneAddressing(); } - return shadowIndex_; + return *zoneAddressingPtr_; } -const Foam::regionCouplePolyPatch& Foam::regionCouplePolyPatch::shadow() const +const Foam::labelList& +Foam::regionCouplePolyPatch::remoteZoneAddressing() const { - return refCast - ( - shadowRegion().boundaryMesh()[shadowIndex()] - ); + if (!remoteZoneAddressingPtr_) + { + calcRemoteZoneAddressing(); + } + + return *remoteZoneAddressingPtr_; +} + + +bool Foam::regionCouplePolyPatch::localParallel() const +{ + // Calculate patch-to-zone addressing + if (!localParallelPtr_) + { + calcLocalParallel(); + } + + return *localParallelPtr_; +} + + +const Foam::ggiZoneInterpolation& +Foam::regionCouplePolyPatch::patchToPatch() const +{ + if (!attached_) + { + FatalErrorIn + ( + "const patchToPatchInterpolation& " + "regionCouplePolyPatch::patchToPatch() const" + ) << "Requesting patchToPatchInterpolation in detached state" + << abort(FatalError); + } + + if (master()) + { + if (!patchToPatchPtr_) + { + Info<< "Initializing the region couple interpolator between " + << "master/shadow patches: " + << name() << " and " << shadowPatchName() + << " on region " << shadowRegionName() + << endl; + + calcPatchToPatch(); + } + + return *patchToPatchPtr_; + } + else + { + return shadow().patchToPatch(); + } } @@ -354,6 +913,24 @@ Foam::regionCouplePolyPatch::reconFaceCellCentres() const void Foam::regionCouplePolyPatch::initAddressing() { + if (active()) + { + // Calculate transforms for correct GGI cut + calcTransforms(); + + // Force zone addressing and remote zone addressing + // (uses GGI interpolator) + zoneAddressing(); + remoteZoneAddressing(); + + // Force local parallel + if (Pstream::parRun() && !localParallel()) + { + // Calculate send addressing + sendAddr(); + } + } + polyPatch::initAddressing(); } @@ -366,6 +943,17 @@ void Foam::regionCouplePolyPatch::calcAddressing() void Foam::regionCouplePolyPatch::initGeometry() { + // Communication is allowed either before or after processor + // patch comms. HJ, 11/Jul/2011 + if (active()) + { + // Note: Only master calculates recon; slave uses master interpolation + if (master()) + { + reconFaceCellCentres(); + } + } + polyPatch::initGeometry(); } @@ -373,13 +961,45 @@ void Foam::regionCouplePolyPatch::initGeometry() void Foam::regionCouplePolyPatch::calcGeometry() { polyPatch::calcGeometry(); - // Reconstruct the cell face centres -// reconFaceCellCentres(); + + // Note: Calculation of transforms must be forced before the + // reconFaceCellCentres in order to correctly set the transformation + // in the interpolation routines. + // HJ, 3/Jul/2009 } void Foam::regionCouplePolyPatch::initMovePoints(const pointField& p) { + clearGeom(); + + // Calculate transforms on mesh motion? + calcTransforms(); + + // Update interpolation for new relative position of GGI interfaces + if (patchToPatchPtr_) + { + patchToPatchPtr_->movePoints(); + } + + // Recalculate send and receive maps + if (active()) + { + // Force zone addressing first + zoneAddressing(); + remoteZoneAddressing(); + + if (Pstream::parRun() && !localParallel()) + { + sendAddr(); + } + } + + if (active() && master()) + { + reconFaceCellCentres(); + } + polyPatch::initMovePoints(p); } @@ -387,14 +1007,6 @@ void Foam::regionCouplePolyPatch::initMovePoints(const pointField& p) void Foam::regionCouplePolyPatch::movePoints(const pointField& p) { polyPatch::movePoints(p); - - // Clear reconstructed face centres - deleteDemandDrivenData(reconFaceCellCentresPtr_); - - if (patchToPatchPtr_) - { - patchToPatchPtr_->movePoints(); - } } @@ -411,21 +1023,16 @@ void Foam::regionCouplePolyPatch::updateMesh() } -void Foam::regionCouplePolyPatch::calcTransformTensors -( - const vectorField& Cf, - const vectorField& Cr, - const vectorField& nf, - const vectorField& nry -) const +void Foam::regionCouplePolyPatch::calcTransforms() { - FatalErrorIn("void regionCouplePolyPatch::calcTransformTensors") - << "Not ready" - << abort(FatalError); + // No transform or separation + forwardT_.setSize(0); + reverseT_.setSize(0); + separation_.setSize(0); } -void Foam::regionCouplePolyPatch::initOrder(const primitivePatch& pp) const +void Foam::regionCouplePolyPatch::initOrder(const primitivePatch&) const {} @@ -459,10 +1066,17 @@ void Foam::regionCouplePolyPatch::write(Ostream& os) const << shadowRegionName_ << token::END_STATEMENT << nl; os.writeKeyword("shadowPatch") << shadowPatchName_ << token::END_STATEMENT << nl; + os.writeKeyword("zone") << zoneName_ + << token::END_STATEMENT << nl; + os.writeKeyword("attached") << attached_ << token::END_STATEMENT << nl; + os.writeKeyword("master") + << master_ << token::END_STATEMENT << nl; os.writeKeyword("isWall") << isWall_ << token::END_STATEMENT << nl; + os.writeKeyword("bridgeOverlap") << bridgeOverlap_ + << token::END_STATEMENT << nl; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H index 1933db121..970295cd9 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H @@ -40,7 +40,8 @@ SourceFiles #define regionCouplePolyPatch_H #include "coupledPolyPatch.H" -#include "patchToPatchInterpolation.H" +#include "ggiInterpolation.H" +#include "faceZone.H" #include "Switch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -66,33 +67,108 @@ class regionCouplePolyPatch //- Shadow patch name const word shadowPatchName_; + //- Interpolation zone name + const word zoneName_; + //- Are the regions attached mutable Switch attached_; + //- Master side. Defined by user, checked for consistency + Switch master_; + //- Are the region attached walls Switch isWall_; + //- Use bridging to fix overlap error in interpolation + Switch bridgeOverlap_; + //- Shadow patch index. Delayed evaluation for construction mutable label shadowIndex_; + //- Interpolation zone index. Delayed evaluation for construction + mutable label zoneIndex_; + //- Patch-to-patch interpolation - mutable patchToPatchInterpolation* patchToPatchPtr_; + mutable ggiZoneInterpolation* patchToPatchPtr_; + + //- Patch-to-zone addressing + mutable labelList* zoneAddressingPtr_; + + //- Remote zone addressing: data needed for interpolation + mutable labelList* remoteZoneAddressingPtr_; //- Reconstructed patch neighbour cell centres mutable vectorField* reconFaceCellCentresPtr_; + // Parallel communication optimisation, stored on master processor + + //- Is the patch localised on a single processor + // (single processor in a parallel run)? + // Used for parallel optimisation + mutable bool* localParallelPtr_; + + //- List of zone faces indices received from each processor + mutable labelListList* receiveAddrPtr_; + + //- List of zone faces indices to send to each processor + mutable labelListList* sendAddrPtr_; + + // Private member functions + //- Calculate patch-to-zone addressing + virtual void calcZoneAddressing() const; + + //- Calculate remote patch-to-zone addressing + virtual void calcRemoteZoneAddressing() const; + //- Calculate interpolation - void calcInterpolation() const; + virtual void calcPatchToPatch() const; + + + // Geometry + + //- Calculate reconstructed cell centres + void calcReconFaceCellCentres() const; + + //- Force calculation of transformation tensors + virtual void calcTransforms(); + + + // Parallel communication optimisation, stored on master processor + + //- Calculate local parallel switch + void calcLocalParallel() const; + + //- Calculate send and receive addressing + void calcSendReceive() const; + + + //- Return receive addressing + const labelListList& receiveAddr() const; + + //- Return send addressing + const labelListList& sendAddr() const; + + + // Memory management + + //- Clear geometry + void clearGeom() const; + + //- Clear out + void clearOut(); protected: // Protected Member functions + //- Is the region couple active? (zone and shadow present) + bool active() const; + //- Initialise the calculation of the patch addressing virtual void initAddressing(); @@ -121,18 +197,6 @@ protected: //- Return non-constant access to shadow patch regionCouplePolyPatch& shadow(); - //- Return reference to patch-to-patch interpolation - const patchToPatchInterpolation& patchToPatch() const; - - //- Calculate reconstructed cell centres - void calcReconFaceCellCentres() const; - - //- Clear geometry - void clearGeom() const; - - //- Clear out - void clearOut() const; - public: @@ -142,6 +206,16 @@ public: // Constructors + //- Construct from components + regionCouplePolyPatch + ( + const word& name, + const label size, + const label start, + const label index, + const polyBoundaryMesh& bm + ); + //- Construct from components regionCouplePolyPatch ( @@ -152,8 +226,11 @@ public: const polyBoundaryMesh& bm, const word& shadowRegionName, const word& shadowPatchName, + const word& zoneName, const bool attached, - const bool attachedWall + const bool master, + const bool isWall, + const bool bridgeOverlap ); //- Construct from dictionary @@ -220,93 +297,143 @@ public: // Member functions - //- Return true if coupled - virtual bool coupled() const; + // Basic info - //- Return shadow region name - const word& shadowRegionName() const - { - return shadowRegionName_; - } + //- Return true if coupled + virtual bool coupled() const; - //- Return shadow patch name - const word& shadowPatchName() const - { - return shadowPatchName_; - } + //- Return shadow region name + const word& shadowRegionName() const + { + return shadowRegionName_; + } - //- Return shadow region - const polyMesh& shadowRegion() const; + //- Return shadow patch name + const word& shadowPatchName() const + { + return shadowPatchName_; + } - //- Return shadow patch index - int shadowIndex() const; + //- Return name of interpolation face zone + const word& zoneName() const + { + return zoneName_; + } - //- Return true if regions are attached - bool attached() const - { - return attached_; - } + //- Return shadow region + const polyMesh& shadowRegion() const; - // Return true if patch is considered a wall - bool isWall() const - { - return isWall_; - } + //- Return shadow patch index + label shadowIndex() const; - //- Attach regions - void attach() const; + //- Return zone patch index + label zoneIndex() const; - //- Attach regions - void detach() const; + // Return true if patch is considered a wall + bool isWall() const + { + return isWall_; + } - //- Return shadow patch - const regionCouplePolyPatch& shadow() const; + //- Return shadow patch + const regionCouplePolyPatch& shadow() const; - //- Interpolate face field - template - tmp > interpolate(const Field& pf) const - { - return patchToPatch().faceInterpolate(pf); - } + //- Return interpolation face zone + const faceZone& zone() const; - template - tmp > interpolate(const tmp >& tpf) const - { - return patchToPatch().faceInterpolate(tpf); - } - //- Return reconstructed cell centres - const vectorField& reconFaceCellCentres() const; + // Edit state - //- Force calculation of transformation tensors - void calcTransformTensors - ( - const vectorField& Cf, - const vectorField& Cr, - const vectorField& nf, - const vectorField& nr - ) const; + //- Attach regions + void attach() const; - //- Initialize ordering for primitivePatch. Does not - // refer to *this (except for name() and type() etc.) - virtual void initOrder(const primitivePatch&) const; + //- Attach regions + void detach() const; - //- Return new ordering for primitivePatch. - // Ordering is -faceMap: for every face - // index of the new face -rotation: for every new face the clockwise - // shift of the original face. Return false if nothing changes - // (faceMap is identity, rotation is 0), true otherwise. - virtual bool order - ( - const primitivePatch&, - labelList& faceMap, - labelList& rotation - ) const; - //- Synchronise communications of ordering for primitivePatch - // Used in cases when no topological change happens locally, - // but is happening on other processors - virtual void syncOrder() const; + // Interpolation data + + //- Is this the master side? + bool master() const + { + return master_; + } + + //- Is this the slave side? + bool slave() const + { + return !master(); + } + + //- Use bridging to fix overlap error in interpolation + bool bridgeOverlap() const + { + return bridgeOverlap_; + } + + //- Return patch-to-zone addressing + const labelList& zoneAddressing() const; + + //- Return remote patch-to-zone addressing + const labelList& remoteZoneAddressing() const; + + //- Is the patch localised on a single processor + bool localParallel() const; + + //- Return reference to patch-to-patch interpolation + const ggiZoneInterpolation& patchToPatch() const; + + + // Interpolation functions + + //- Expand face field to full zone size, including reduction + template + tmp > fastExpand(const Field& pf) const; + + //- Interpolate face field + template + tmp > interpolate(const Field& pf) const; + + template + tmp > interpolate(const tmp >& tpf) const; + + //- Bridge interpolated face field for uncovered faces + template + void bridge + ( + const Field& bridgeField, + Field& ff + ) const; + + + // Geometric data + + //- Return reconstructed cell centres + const vectorField& reconFaceCellCentres() const; + + + // Patch ordering + + //- Initialize ordering for primitivePatch. Does not + // refer to *this (except for name() and type() etc.) + virtual void initOrder(const primitivePatch&) const; + + //- Return new ordering for primitivePatch. + // Ordering is -faceMap: for every face + // index of the new face -rotation: for every new face the + // clockwise shift of the original face. Return false if nothing + // changes (faceMap is identity, rotation is 0), true otherwise. + virtual bool order + ( + const primitivePatch&, + labelList& faceMap, + labelList& rotation + ) const; + + //- Synchronise communications of ordering for primitivePatch + // Used in cases when no topological change happens locally, + // but is happening on other processors + virtual void syncOrder() const; //- Write @@ -320,6 +447,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "regionCouplePolyPatchTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatchTemplates.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatchTemplates.C new file mode 100644 index 000000000..617af8cf0 --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatchTemplates.C @@ -0,0 +1,350 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Hrvoje Jasak, Wikki Ltd. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "regionCouplePolyPatch.H" +#include "OSspecific.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::tmp > Foam::regionCouplePolyPatch::fastExpand +( + const Field& ff +) const +{ + // Check and expand the field from patch size to zone size + // with communication + + // Algorithm: + // 1) Master processor holds maps of all zone addressing (data provided) + // and all remote zone addressing (data required) + // 2) Each processor will send the locally active data to the master + // 3) Master assembles all the data + // 4) Master sends to all processors the data they need to receive + // + // Notes: + // A) If the size of zone addressing is zero, data is not sent + // B) Communicated data on each processor has the size of live faces + // C) Expanded data will be equal to actual data fronm other processors + // only for the faces marked in remote; for other faces, it will be + // equal to zero + // D) On processor zero, complete data is available + // HJ, 4/Jun/2011 + if (ff.size() != size()) + { + FatalErrorIn + ( + "tmp > regionCouplePolyPatch::fastExpand\n" + "(\n" + " const Field& ff\n" + ") const" + ) << "Incorrect patch field size. Field size: " + << ff.size() << " patch size: " << size() + << abort(FatalError); + } + + if (localParallel()) + { + FatalErrorIn + ( + "tmp > regionCouplePolyPatch::fastExpand" + "(" + " const Field& ff" + ") const" + ) << "Requested expand on local parallel. This is not allowed" + << abort(FatalError); + } + + // Expand the field to zone size + tmp > texpandField + ( + new Field(zone().size(), pTraits::zero) + ); + + Field& expandField = texpandField(); + + if (Pstream::master()) + { + // Insert master processor + const labelList& za = zoneAddressing(); + + forAll (za, i) + { + expandField[za[i]] = ff[i]; + } + + // Master receives and inserts data from all processors for which + // receiveAddr contains entries + for (label procI = 1; procI < Pstream::nProcs(); procI++) + { + const labelList& curRAddr = receiveAddr()[procI]; + + if (!curRAddr.empty()) + { + Field receiveBuf(curRAddr.size()); + + // Opt: reconsider mode of communication + IPstream::read + ( + Pstream::blocking, + procI, + reinterpret_cast(receiveBuf.begin()), + receiveBuf.byteSize() + ); + + // Insert received information + forAll (curRAddr, i) + { + expandField[curRAddr[i]] = receiveBuf[i]; + } + } + } + + // Expanded field complete, send required data to other processors + for (label procI = 1; procI < Pstream::nProcs(); procI++) + { + const labelList& curSAddr = shadow().sendAddr()[procI]; + + if (!curSAddr.empty()) + { + Field sendBuf(curSAddr.size()); + + forAll (curSAddr, i) + { + sendBuf[i] = expandField[curSAddr[i]]; + } + + // Opt: reconsider mode of communication + OPstream::write + ( + Pstream::blocking, + procI, + reinterpret_cast(sendBuf.begin()), + sendBuf.byteSize() + ); + } + } + } + else + { + // Send local data to master and receive remote data + // If patch is empty, communication is avoided + // HJ, 4/Jun/2011 + if (size()) + { + // Opt: reconsider mode of communication + OPstream::write + ( + Pstream::blocking, + Pstream::masterNo(), + reinterpret_cast(ff.begin()), + ff.byteSize() + ); + } + + // Prepare to receive remote data + const labelList& rza = shadow().remoteZoneAddressing(); + + if (!rza.empty()) + { + Field receiveBuf(rza.size()); + + // Opt: reconsider mode of communication + IPstream::read + ( + Pstream::blocking, + Pstream::masterNo(), + reinterpret_cast(receiveBuf.begin()), + receiveBuf.byteSize() + ); + + // Insert the data into expanded field + forAll (rza, i) + { + expandField[rza[i]] = receiveBuf[i]; + } + } + } + + return texpandField; +} + + +template +Foam::tmp > Foam::regionCouplePolyPatch::interpolate +( + const Field& ff +) const +{ + // Check and expand the field from patch size to zone size + if (ff.size() != shadow().size()) + { + FatalErrorIn + ( + "tmp > regionCouplePolyPatch::interpolate\n" + "(\n" + " const Field& ff\n" + ") const" + ) << "Incorrect slave patch field size. Field size: " + << ff.size() << " patch size: " << shadow().size() + << abort(FatalError); + } + + if (localParallel()) + { + // No expansion or filtering needed. HJ, 4/Jun/2011 + + if (empty()) + { + // Patch empty, no interpolation + return tmp >(new Field()); + } + + // Interpolate field + if (master()) + { + return patchToPatch().slaveToMaster(ff); + } + else + { + return patchToPatch().masterToSlave(ff); + } + } + else + { + // Expand shadow + Field expandField = shadow().fastExpand(ff); + + tmp > tresult(new Field(size())); + Field& result = tresult(); + + if (master()) + { + patchToPatch().maskedSlaveToMaster + ( + expandField, + result, + zoneAddressing() + ); + } + else + { + patchToPatch().maskedMasterToSlave + ( + expandField, + result, + zoneAddressing() + ); + } + + return tresult; + } +} + + +template +Foam::tmp > Foam::regionCouplePolyPatch::interpolate +( + const tmp >& tff +) const +{ + tmp > tint = interpolate(tff()); + tff.clear(); + return tint; +} + + +template +void Foam::regionCouplePolyPatch::bridge +( + const Field& bridgeField, + Field& ff +) const +{ + // Check and expand the field from patch size to zone size + if (ff.size() != size()) + { + FatalErrorIn + ( + "tmp > regionCouplePolyPatch::bridge\n" + "(\n" + " const Field& ff,\n" + " Field& ff\n" + ") const" + ) << "Incorrect patch field size for bridge. Field size: " + << ff.size() << " patch size: " << size() + << abort(FatalError); + } + + if (bridgeOverlap()) + { + if (empty()) + { + // Patch empty, no bridging + return; + } + + if (localParallel()) + { + if (master()) + { + patchToPatch().bridgeMaster(bridgeField, ff); + } + else + { + patchToPatch().bridgeSlave(bridgeField, ff); + } + } + else + { + // Note: since bridging is only a local operation + if (master()) + { + patchToPatch().maskedBridgeMaster + ( + bridgeField, + ff, + zoneAddressing() + ); + } + else + { + patchToPatch().maskedBridgeSlave + ( + bridgeField, + ff, + zoneAddressing() + ); + } + } + } +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 5cdcfe904..0f85e7702 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -110,7 +110,7 @@ $(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C $(constraintFvPatchFields)/ggi/ggiFvPatchFields.C $(constraintFvPatchFields)/cyclicGgi/cyclicGgiFvPatchFields.C $(constraintFvPatchFields)/overlapGgi/overlapGgiFvPatchFields.C -$(constraintFvPatchFields)/regionCouple/regionCoupleFvPatchFields.C +$(constraintFvPatchFields)/regionCoupling/regionCouplingFvPatchFields.C derivedFvPatchFields = $(fvPatchFields)/derived $(derivedFvPatchFields)/activeBaffleVelocity/activeBaffleVelocityFvPatchVectorField.C @@ -185,7 +185,7 @@ $(constraintFvsPatchFields)/wedge/wedgeFvsPatchFields.C $(constraintFvsPatchFields)/ggi/ggiFvsPatchFields.C $(constraintFvsPatchFields)/cyclicGgi/cyclicGgiFvsPatchFields.C $(constraintFvsPatchFields)/overlapGgi/overlapGgiFvsPatchFields.C -$(constraintFvsPatchFields)/regionCouple/regionCoupleFvsPatchFields.C +$(constraintFvsPatchFields)/regionCoupling/regionCouplingFvsPatchFields.C fields/volFields/volFields.C @@ -227,6 +227,7 @@ $(schemes)/skewCorrected/skewCorrected.C $(schemes)/outletStabilised/outletStabilised.C $(schemes)/reverseLinear/reverseLinear.C $(schemes)/clippedLinear/clippedLinear.C +$(schemes)/harmonic/magLongDelta.C $(schemes)/harmonic/harmonic.C $(schemes)/fixedBlended/fixedBlended.C $(schemes)/localBlended/localBlended.C diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C index f697e46d5..2903116c0 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C @@ -140,6 +140,15 @@ ggiFvPatchField::ggiFvPatchField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// Return shadow field +template +const Foam::ggiFvPatchField& +ggiFvPatchField::shadowPatchField() const +{ + return this->boundaryField()[ggiPatch_.shadowIndex()]; +} + + // Return neighbour field template tmp > ggiFvPatchField::patchNeighbourField() const @@ -147,6 +156,9 @@ tmp > ggiFvPatchField::patchNeighbourField() const const Field& iField = this->internalField(); // Get shadow face-cells and assemble shadow field + // This is a patchInternalField of neighbour but access is inconvenient. + // Assemble by hand. + // HJ, 27/Sep/2011 const unallocLabelList& sfc = ggiPatch_.shadow().faceCells(); Field sField(sfc.size()); diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H index 32be042a8..d2a88a80d 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.H @@ -127,6 +127,12 @@ public: // Member functions + // Access + + //- Return shadow patch field + const ggiFvPatchField& shadowPatchField() const; + + // Evaluation functions //- Return neighbour field given internal cell data diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C deleted file mode 100644 index 33a7d3d7c..000000000 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C +++ /dev/null @@ -1,275 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / 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 - Hrvoje Jasak, Wikki Ltd. All rights reserved - -\*---------------------------------------------------------------------------*/ - -#include "regionCoupleFvPatchField.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -template -regionCoupleFvPatchField::regionCoupleFvPatchField -( - const fvPatch& p, - const DimensionedField& iF -) -: - coupledFvPatchField(p, iF), - regionCouplePatch_(refCast(p)), - remoteFieldName_(iF.name()), - matrixUpdateBuffer_() -{} - - -template -regionCoupleFvPatchField::regionCoupleFvPatchField -( - const fvPatch& p, - const DimensionedField& iF, - const dictionary& dict -) -: - coupledFvPatchField(p, iF, dict), - regionCouplePatch_(refCast(p)), - remoteFieldName_(dict.lookup("remoteField")), - matrixUpdateBuffer_() -{ - if (!isType(p)) - { - FatalIOErrorIn - ( - "regionCoupleFvPatchField::regionCoupleFvPatchField\n" - "(\n" - " const fvPatch& p,\n" - " const DimensionedField& iF,\n" - " const dictionary& dict\n" - ")\n", - dict - ) << "patch " << this->patch().index() << " not regionCouple type. " - << "Patch type = " << p.type() - << exit(FatalIOError); - } -} - - -template -regionCoupleFvPatchField::regionCoupleFvPatchField -( - const regionCoupleFvPatchField& ptf, - const fvPatch& p, - const DimensionedField& iF, - const fvPatchFieldMapper& mapper -) -: - coupledFvPatchField(ptf, p, iF, mapper), - regionCouplePatch_(refCast(p)), - remoteFieldName_(ptf.remoteFieldName_), - matrixUpdateBuffer_() -{ - if (!isType(this->patch())) - { - FatalErrorIn - ( - "regionCoupleFvPatchField::regionCoupleFvPatchField\n" - "(\n" - " const regionCoupleFvPatchField& ptf,\n" - " const fvPatch& p,\n" - " const DimensionedField& iF,\n" - " const fvPatchFieldMapper& mapper\n" - ")\n" - ) << "Field type does not correspond to patch type for patch " - << this->patch().index() << "." << endl - << "Field type: " << typeName << endl - << "Patch type: " << this->patch().type() - << exit(FatalError); - } -} - - -template -regionCoupleFvPatchField::regionCoupleFvPatchField -( - const regionCoupleFvPatchField& ptf, - const DimensionedField& iF -) -: - regionCoupleLduInterfaceField(), - coupledFvPatchField(ptf, iF), - regionCouplePatch_(refCast(ptf.patch())), - remoteFieldName_(ptf.remoteFieldName_), - matrixUpdateBuffer_() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -// Return neighbour field -template -const regionCoupleFvPatchField& -regionCoupleFvPatchField::shadowPatchField() const -{ - // Lookup neighbour field - typedef GeometricField GeoField; - - const GeoField& coupleField = - regionCouplePatch_.shadowRegion(). - objectRegistry::lookupObject(remoteFieldName_); - - return refCast > - ( - coupleField.boundaryField()[regionCouplePatch_.shadowIndex()] - ); -} - - -// Return neighbour field -template -tmp > regionCoupleFvPatchField::patchNeighbourField() const -{ - return regionCouplePatch_.interpolate - ( - shadowPatchField().patchInternalField() - ); -} - - -template -void regionCoupleFvPatchField::evaluate -( - const Pstream::commsTypes -) -{ - // Implement weights-based stabilised harmonic interpolation using - // magnitude of type - // Algorithm: - // 1) calculate magnitude of internal field and neighbour field - // 2) calculate harmonic mean magnitude - // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei - // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei) - // 5) Use weights to interpolate values - - Field fOwn = this->patchInternalField(); - Field fNei = this->patchNeighbourField(); - - scalarField magFOwn = mag(fOwn); - scalarField magFNei = mag(fNei); - - // Calculate internal weights using field magnitude - scalarField weights(fOwn.size()); - - forAll (weights, faceI) - { - scalar mOwn = magFOwn[faceI]; - scalar mNei = magFNei[faceI]; - - scalar den1 = mOwn + mNei; - scalar den2 = mOwn - mNei; - - // Strange optimisation check: floating point failure if only - // checking for den2. den1 check should not be needed - // HJ, 24/Aug/2011 - if (mag(den1) > SMALL && mag(den2) > SMALL) - { - scalar mean = 2.0*mOwn*mNei/den1; - weights[faceI] = (mean - mNei)/den2; - } - else - { - weights[faceI] = 0.5; - } - } - - // Do interpolation - Field::operator=(weights*fOwn + (1.0 - weights)*fNei); -} - - -// Initialise neighbour processor internal cell data -template -void regionCoupleFvPatchField::initInterfaceMatrixUpdate -( - const scalarField& psiInternal, - scalarField&, - const lduMatrix&, - const scalarField&, - const direction, - const Pstream::commsTypes -) const -{ - matrixUpdateBuffer_ = this->patch().patchInternalField(psiInternal); -} - - -// Return matrix product for coupled boundary -template -void regionCoupleFvPatchField::updateInterfaceMatrix -( - const scalarField& psiInternal, - scalarField& result, - const lduMatrix&, - const scalarField& coeffs, - const direction cmpt, - const Pstream::commsTypes -) const -{ - scalarField pnf = - regionCouplePatch_.interpolate - ( - this->shadowPatchField().matrixUpdateBuffer() - ); - - // Multiply the field by coefficients and add into the result - const unallocLabelList& fc = regionCouplePatch_.faceCells(); - - forAll(fc, elemI) - { - result[fc[elemI]] -= coeffs[elemI]*pnf[elemI]; - } -} - - -// Write -template -void regionCoupleFvPatchField::write(Ostream& os) const -{ - fvPatchField::write(os); - os.writeKeyword("remoteField") - << remoteFieldName_ << token::END_STATEMENT << nl; - this->writeEntry("value", os); -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.C new file mode 100644 index 000000000..b47ccecd7 --- /dev/null +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.C @@ -0,0 +1,468 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + Hrvoje Jasak, Wikki Ltd. All rights reserved + +\*---------------------------------------------------------------------------*/ + +#include "regionCouplingFvPatchField.H" +#include "symmTransformField.H" +#include "magLongDelta.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +regionCouplingFvPatchField::regionCouplingFvPatchField +( + const fvPatch& p, + const DimensionedField& iF +) +: + coupledFvPatchField(p, iF), + regionCouplePatch_(refCast(p)), + remoteFieldName_(iF.name()), + matrixUpdateBuffer_(), + originalPatchField_(), + curTimeIndex_(-1) +{} + + +template +regionCouplingFvPatchField::regionCouplingFvPatchField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + coupledFvPatchField(p, iF, dict), + regionCouplePatch_(refCast(p)), + remoteFieldName_(dict.lookup("remoteField")), + matrixUpdateBuffer_(), + originalPatchField_(), + curTimeIndex_(-1) +{ + if (!isType(p)) + { + FatalIOErrorIn + ( + "regionCouplingFvPatchField::regionCouplingFvPatchField\n" + "(\n" + " const fvPatch& p,\n" + " const DimensionedField& iF,\n" + " const dictionary& dict\n" + ")\n", + dict + ) << "patch " << this->patch().index() << " not regionCouple type. " + << "Patch type = " << p.type() + << exit(FatalIOError); + } + + if (!dict.found("value")) + { + // Grab the internal value for initialisation. (?) HJ, 27/Feb/2009 + fvPatchField::operator=(this->patchInternalField()()); + } +} + + +template +regionCouplingFvPatchField::regionCouplingFvPatchField +( + const regionCouplingFvPatchField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + coupledFvPatchField(ptf, p, iF, mapper), + regionCouplePatch_(refCast(p)), + remoteFieldName_(ptf.remoteFieldName_), + matrixUpdateBuffer_(), + originalPatchField_(), + curTimeIndex_(-1) +{ + if (!isType(this->patch())) + { + FatalErrorIn + ( + "regionCouplingFvPatchField::regionCouplingFvPatchField\n" + "(\n" + " const regionCouplingFvPatchField& ptf,\n" + " const fvPatch& p,\n" + " const DimensionedField& iF,\n" + " const fvPatchFieldMapper& mapper\n" + ")\n" + ) << "Field type does not correspond to patch type for patch " + << this->patch().index() << "." << endl + << "Field type: " << typeName << endl + << "Patch type: " << this->patch().type() + << exit(FatalError); + } +} + + +template +regionCouplingFvPatchField::regionCouplingFvPatchField +( + const regionCouplingFvPatchField& ptf, + const DimensionedField& iF +) +: + ggiLduInterfaceField(), + coupledFvPatchField(ptf, iF), + regionCouplePatch_(refCast(ptf.patch())), + remoteFieldName_(ptf.remoteFieldName_), + matrixUpdateBuffer_(), + originalPatchField_(), + curTimeIndex_(-1) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +// Return a named shadow patch field +template +template +const typename GeometricField::PatchFieldType& +regionCouplingFvPatchField::lookupShadowPatchField +( + const word& name, + const GeometricField*, + const Type2* +) const +{ + // Lookup neighbour field + const GeometricField& shadowField = + regionCouplePatch_.shadowRegion(). + objectRegistry::lookupObject(name); + + return shadowField.boundaryField()[regionCouplePatch_.shadowIndex()]; +} + + +// Return shadow patch field +template +const regionCouplingFvPatchField& +regionCouplingFvPatchField::shadowPatchField() const +{ + // Lookup neighbour field + typedef GeometricField GeoField; + + return refCast > + ( + lookupShadowPatchField(remoteFieldName_) + ); +} + + +// Return neighbour field +template +tmp > regionCouplingFvPatchField::patchNeighbourField() const +{ + Field sField = shadowPatchField().patchInternalField(); + + tmp > tpnf + ( + regionCouplePatch_.interpolate + ( + shadowPatchField().patchInternalField() + ) + ); + + Field& pnf = tpnf(); + + if (regionCouplePatch_.bridgeOverlap()) + { + // Symmetry treatment used for overlap + vectorField nHat = this->patch().nf(); + + // Use mirrored neighbour field for interpolation + // HJ, 21/Jan/2009 + Field bridgeField = + transform(I - 2.0*sqr(nHat), this->patchInternalField()); + + regionCouplePatch_.bridge(bridgeField, pnf); + } + + return tpnf; +} + + +// Return neighbour field given internal cell data +template +tmp > regionCouplingFvPatchField::patchNeighbourField +( + const word& name +) const +{ + // Lookup neighbour field + typedef GeometricField GeoField; + + return regionCouplePatch_.interpolate + ( + lookupShadowPatchField(name).patchInternalField() + ); + + // Note: this field is not bridged because local data does not exist + // for named field. HJ, 27/Sep/2011 +} + + +template +void regionCouplingFvPatchField::initEvaluate +( + const Pstream::commsTypes commsType +) +{ + if (!this->updated()) + { + this->updateCoeffs(); + } + + // Interpolation must happen at init + // Implement weights-based stabilised harmonic interpolation using + // magnitude of type + // Algorithm: + // 1) calculate magnitude of internal field and neighbour field + // 2) calculate harmonic mean magnitude + // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei + // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei) + // 5) Use weights to interpolate values + + const Field& fOwn = this->originalPatchField(); + Field fNei = regionCouplePatch_.interpolate + ( + this->shadowPatchField().originalPatchField() + ); + + // Larger small for complex arithmetic accuracy + const scalar kSmall = 1000*SMALL; + +# if 0 + // Hrv's treatment + scalarField mOwn = mag(fOwn); + scalarField mNei = mag(fNei); + scalarField mean = 2*(mOwn*mNei)/(mOwn + mNei); + + scalarField weights(fOwn.size(), 0.5); + + scalar den; + + forAll (weights, faceI) + { + den = (mNei[faceI] - mOwn[faceI]); + + // Note: complex arithmetic requires extra accuracy + // This is a division of two close subtractions + // HJ, 28/Sep/2011 + if (mag(den) > kSmall) + { + weights[faceI] = (mNei[faceI] - mean[faceI])/den; + } + else + { + // Use 0.5 weights + } + } + +# else + + // Henrik's treatment + const fvPatch& p = this->patch(); + + // Note: for interpolation, work with face fields, to allow wall-corrected + // diffusivity (eg wall functions) to operate correctly. + // HJ, 28/Sep/2011 + Field f = *this; + const regionCouplingFvPatchField& spf = shadowPatchField(); + const fvPatch& sp = spf.patch(); + + // Mag long deltas are identical on both sides. HJ, 28/Sep/2011 + const magLongDelta& mld = magLongDelta::New(p.boundaryMesh().mesh()); + + scalarField magPhiOwn = mag(fOwn); + scalarField magPhiNei = mag(fNei); + + const scalarField& pWeights = p.weights(); + const scalarField& pDeltaCoeffs = p.deltaCoeffs(); + const scalarField& pLongDelta = mld.magDelta(p.index()); + + // Calculate internal weights using field magnitude + scalarField weights(fOwn.size()); + + forAll (weights, faceI) + { + scalar mOwn = magPhiOwn[faceI]/(1 - pWeights[faceI]); + scalar mNei = magPhiNei[faceI]/pWeights[faceI]; + + scalar den = magPhiNei[faceI] - magPhiOwn[faceI]; + + // Note: complex arithmetic requires extra accuracy + // This is a division of two close subtractions + // HJ, 28/Sep/2011 + if (mag(den) > kSmall) + { + scalar mean = mOwn*mNei/ + ( + (mOwn + mNei)* + pLongDelta[faceI]* + pDeltaCoeffs[faceI] + ); + + weights[faceI] = (magPhiNei[faceI] - mean)/den; + } + else + { + weights[faceI] = 0.5; + } + } + +#endif + + // Do interpolation + Field::operator=(weights*fOwn + (1.0 - weights)*fNei); +} + + +template +void regionCouplingFvPatchField::evaluate +( + const Pstream::commsTypes +) +{ + // No interpolation allowed + + fvPatchField::evaluate(); +} + + +template +void regionCouplingFvPatchField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + // Store original field for symmetric evaluation + // Henrik Rusche, Aug/2011 + if (curTimeIndex_ != this->db().time().timeIndex()) + { + originalPatchField_ = *this; + curTimeIndex_ = this->db().time().timeIndex(); + } +} + + +template +tmp > regionCouplingFvPatchField::snGrad() const +{ + if (regionCouplePatch_.coupled()) + { + return coupledFvPatchField::snGrad(); + } + else + { + return fvPatchField::snGrad(); + } +} + + +// Initialise neighbour processor internal cell data +template +void regionCouplingFvPatchField::initInterfaceMatrixUpdate +( + const scalarField& psiInternal, + scalarField& result, + const lduMatrix&, + const scalarField& coeffs, + const direction cmpt, + const Pstream::commsTypes +) const +{ + // Prepare local matrix update buffer for the remote side. + // Note that only remote side will have access to its psiInternal + // as they are on different regions + + // Since interpolation needs to happen on the shadow, and within the + // init, prepare interpolation for the other side. + matrixUpdateBuffer_ = + this->shadowPatchField().regionCouplePatch().interpolate + ( + this->patch().patchInternalField(psiInternal) + ); +} + + +// Return matrix product for coupled boundary +template +void regionCouplingFvPatchField::updateInterfaceMatrix +( + const scalarField& psiInternal, + scalarField& result, + const lduMatrix&, + const scalarField& coeffs, + const direction , + const Pstream::commsTypes +) const +{ + // Note: interpolation involves parallel communications and needs to + // happen during init. This changes the use of matrix update buffer + // compared to earlier versions + // HJ, 28/Sep/2011 + scalarField pnf = this->shadowPatchField().matrixUpdateBuffer(); + + // Multiply the field by coefficients and add into the result + const unallocLabelList& fc = regionCouplePatch_.faceCells(); + + forAll(fc, elemI) + { + result[fc[elemI]] -= coeffs[elemI]*pnf[elemI]; + } +} + + +// Write +template +void regionCouplingFvPatchField::write(Ostream& os) const +{ + fvPatchField::write(os); + os.writeKeyword("remoteField") + << remoteFieldName_ << token::END_STATEMENT << nl; + this->writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.H similarity index 60% rename from src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.H rename to src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.H index 8d608d1e8..bce85b7c3 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCoupleFvPatchField + regionCouplingFvPatchField Description Region couple patch field @@ -32,15 +32,15 @@ Author Hrvoje Jasak, Wikki Ltd. All rights reserved SourceFiles - regionCoupleFvPatchField.C + regionCouplingFvPatchField.C \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvPatchField_H -#define regionCoupleFvPatchField_H +#ifndef regionCouplingFvPatchField_H +#define regionCouplingFvPatchField_H #include "coupledFvPatchField.H" -#include "regionCoupleLduInterfaceField.H" +#include "ggiLduInterfaceField.H" #include "regionCoupleFvPatch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,13 +49,13 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class regionCoupleFvPatchField Declaration + Class regionCouplingFvPatchField Declaration \*---------------------------------------------------------------------------*/ template -class regionCoupleFvPatchField +class regionCouplingFvPatchField : - public regionCoupleLduInterfaceField, + public ggiLduInterfaceField, public coupledFvPatchField { // Private data @@ -69,9 +69,29 @@ class regionCoupleFvPatchField //- Matrix update buffer mutable scalarField matrixUpdateBuffer_; + //- Original patch field. Required for correct evaluation + // in harmonic averaging + Field originalPatchField_; + + //- Current time index used to store originalPatchField_ + label curTimeIndex_; + protected: + //- Return original patch field + const Field& originalPatchField() const + { + if (curTimeIndex_ == this->db().time().timeIndex()) + { + return originalPatchField_; + } + else + { + return *this; + } + } + //- Return contents of a matrix update buffer const scalarField& matrixUpdateBuffer() const { @@ -82,40 +102,40 @@ protected: public: //- Runtime type information -// TypeName(regionCoupleFvPatch::typeName_()); TypeName("regionCoupling"); // Constructors //- Construct from patch and internal field - regionCoupleFvPatchField + regionCouplingFvPatchField ( const fvPatch&, const DimensionedField& ); //- Construct from patch, internal field and dictionary - regionCoupleFvPatchField + regionCouplingFvPatchField ( const fvPatch&, const DimensionedField&, const dictionary& ); - //- Construct by mapping given regionCoupleFvPatchField onto a new patch - regionCoupleFvPatchField + //- Construct by mapping given regionCouplingFvPatchField + // onto a new patch + regionCouplingFvPatchField ( - const regionCoupleFvPatchField&, + const regionCouplingFvPatchField&, const fvPatch&, const DimensionedField&, const fvPatchFieldMapper& ); //- Construct as copy setting internal field reference - regionCoupleFvPatchField + regionCouplingFvPatchField ( - const regionCoupleFvPatchField&, + const regionCouplingFvPatchField&, const DimensionedField& ); @@ -124,7 +144,7 @@ public: { return tmp > ( - new regionCoupleFvPatchField(*this) + new regionCouplingFvPatchField(*this) ); } @@ -136,14 +156,14 @@ public: { return tmp > ( - new regionCoupleFvPatchField(*this, iF) + new regionCouplingFvPatchField(*this, iF) ); } // Member functions - // Evaluation functions + // Access functions //- Return remote field name const word& remoteFieldName() const @@ -151,18 +171,52 @@ public: return remoteFieldName_; } + //- Return a named shadow patch field + template + const typename GeometricField::PatchFieldType& + lookupShadowPatchField + ( + const word& name, + const GeometricField* = NULL, + const Type2* = NULL + ) const; + + + const regionCoupleFvPatch& regionCouplePatch() const + { + return regionCouplePatch_; + } + //- Return shadow patch field - const regionCoupleFvPatchField& shadowPatchField() const; + const regionCouplingFvPatchField& shadowPatchField() const; + + + // Evaluation functions //- Return neighbour field given internal cell data virtual tmp > patchNeighbourField() const; + //- Return neighbour field given internal cell data + virtual tmp > patchNeighbourField + ( + const word& name + ) const; + + //- Initialise the evaluation of the patch field + virtual void initEvaluate(const Pstream::commsTypes commsType); + //- Evaluate the patch field virtual void evaluate ( const Pstream::commsTypes commsType ); + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Return patch-normal gradient + virtual tmp > snGrad() const; + //- Initialise neighbour matrix update virtual void initInterfaceMatrixUpdate ( @@ -186,6 +240,37 @@ public: ) const; + // GGI coupled interface functions + + //- Does the patch field perform the transfromation + virtual bool doTransform() const + { + return + !( + regionCouplePatch_.parallel() + || pTraits::rank == 0 + ); + } + + //- Return face transformation tensor + virtual const tensorField& forwardT() const + { + return regionCouplePatch_.forwardT(); + } + + //- Return neighbour-cell transformation tensor + virtual const tensorField& reverseT() const + { + return regionCouplePatch_.reverseT(); + } + + //- Return rank of component for transform + virtual int rank() const + { + return pTraits::rank; + } + + //- Write virtual void write(Ostream&) const; @@ -199,7 +284,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository -# include "regionCoupleFvPatchField.C" +# include "regionCouplingFvPatchField.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.C similarity index 95% rename from src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.C rename to src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.C index 927ae6764..7c5b77c0c 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.C @@ -30,7 +30,7 @@ Author \*---------------------------------------------------------------------------*/ -#include "regionCoupleFvPatchFields.H" +#include "regionCouplingFvPatchFields.H" #include "addToRunTimeSelectionTable.H" #include "volFields.H" @@ -41,7 +41,7 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -makePatchFields(regionCouple); +makePatchFields(regionCoupling); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.H b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.H similarity index 88% rename from src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.H rename to src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.H index b7e1359e9..a9768bce8 100644 --- a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFields.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCouplewFvsPatchFields + regionCouplingFvPatchFields Description Region couple patch field @@ -32,14 +32,14 @@ Author Hrvoje Jasak, Wikki Ltd. All rights reserved SourceFiles - regionCoupleFvsPatchFields.C + regionCouplingFvPatchFields.C \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvsPatchFields_H -#define regionCoupleFvsPatchFields_H +#ifndef regionCouplingFvPatchFields_H +#define regionCouplingFvPatchFields_H -#include "regionCoupleFvsPatchField.H" +#include "regionCouplingFvPatchField.H" #include "fieldTypes.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +49,7 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -makeFvsPatchTypeFieldTypedefs(regionCouple) +makePatchTypeFieldTypedefs(regionCoupling) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFieldsFwd.H b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFieldsFwd.H similarity index 89% rename from src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFieldsFwd.H rename to src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFieldsFwd.H index 378447df1..43b764597 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFieldsFwd.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchFieldsFwd.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCoupleFvPatchField + regionCouplingFvPatchField Description Region couple patch field @@ -33,8 +33,8 @@ Author \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvPatchFieldsFwd_H -#define regionCoupleFvPatchFieldsFwd_H +#ifndef regionCouplingFvPatchFieldsFwd_H +#define regionCouplingFvPatchFieldsFwd_H #include "fvPatchField.H" #include "fieldTypes.H" @@ -46,9 +46,9 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -template class regionCoupleFvPatchField; +template class regionCouplingFvPatchField; -makePatchTypeFieldTypedefs(regionCouple) +makePatchTypeFieldTypedefs(regionCoupling) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.C b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.C similarity index 84% rename from src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.C rename to src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.C index 1401b7d9a..ebb73d11b 100644 --- a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.C +++ b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.C @@ -27,7 +27,7 @@ Author \*---------------------------------------------------------------------------*/ -#include "regionCoupleFvsPatchField.H" +#include "regionCouplingFvsPatchField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -37,7 +37,7 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -regionCoupleFvsPatchField::regionCoupleFvsPatchField +regionCouplingFvsPatchField::regionCouplingFvsPatchField ( const fvPatch& p, const DimensionedField& iF @@ -51,7 +51,7 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField template -regionCoupleFvsPatchField::regionCoupleFvsPatchField +regionCouplingFvsPatchField::regionCouplingFvsPatchField ( const fvPatch& p, const DimensionedField& iF, @@ -67,7 +67,7 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField { FatalIOErrorIn ( - "regionCoupleFvsPatchField::regionCoupleFvsPatchField\n" + "regionCouplingFvsPatchField::regionCouplingFvsPatchField\n" "(\n" " const fvPatch& p,\n" " const DimensionedField& iF,\n" @@ -82,9 +82,9 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField template -regionCoupleFvsPatchField::regionCoupleFvsPatchField +regionCouplingFvsPatchField::regionCouplingFvsPatchField ( - const regionCoupleFvsPatchField& ptf, + const regionCouplingFvsPatchField& ptf, const fvPatch& p, const DimensionedField& iF, const fvPatchFieldMapper& mapper @@ -99,9 +99,9 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField { FatalErrorIn ( - "regionCoupleFvsPatchField::regionCoupleFvsPatchField\n" + "regionCouplingFvsPatchField::regionCouplingFvsPatchField\n" "(\n" - " const regionCoupleFvsPatchField& ptf,\n" + " const regionCouplingFvsPatchField& ptf,\n" " const fvPatch& p,\n" " const DimensionedField& iF,\n" " const fvPatchFieldMapper& mapper\n" @@ -116,9 +116,9 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField template -regionCoupleFvsPatchField::regionCoupleFvsPatchField +regionCouplingFvsPatchField::regionCouplingFvsPatchField ( - const regionCoupleFvsPatchField& ptf, + const regionCouplingFvsPatchField& ptf, const DimensionedField& iF ) : @@ -133,8 +133,8 @@ regionCoupleFvsPatchField::regionCoupleFvsPatchField // Return neighbour field template -const regionCoupleFvsPatchField& -regionCoupleFvsPatchField::shadowPatchField() const +const regionCouplingFvsPatchField& +regionCouplingFvsPatchField::shadowPatchField() const { // Lookup neighbour field typedef GeometricField GeoField; @@ -143,7 +143,7 @@ regionCoupleFvsPatchField::shadowPatchField() const regionCouplePatch_.shadowRegion(). objectRegistry::lookupObject(remoteFieldName_); - return refCast > + return refCast > ( coupleField.boundaryField()[regionCouplePatch_.shadowIndex()] ); @@ -152,7 +152,7 @@ regionCoupleFvsPatchField::shadowPatchField() const // Write template -void regionCoupleFvsPatchField::write(Ostream& os) const +void regionCouplingFvsPatchField::write(Ostream& os) const { fvsPatchField::write(os); os.writeKeyword("remoteField") diff --git a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.H b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.H similarity index 80% rename from src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.H rename to src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.H index ea28878cd..ea329ae72 100644 --- a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchField.H +++ b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchField.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCoupleFvsPatchField + regionCouplingFvsPatchField Description Region couple patch field @@ -31,13 +31,17 @@ Description Author Hrvoje Jasak, Wikki Ltd. All rights reserved +Note + The name of the class needs to be different from the name of the patch + (regionCouple) to avoid constraint patch behaviour + SourceFiles - regionCoupleFvsPatchField.C + regionCouplingFvsPatchField.C \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvsPatchField_H -#define regionCoupleFvsPatchField_H +#ifndef regionCouplingFvsPatchField_H +#define regionCouplingFvsPatchField_H #include "coupledFvsPatchField.H" #include "regionCoupleLduInterfaceField.H" @@ -49,17 +53,17 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class regionCoupleFvsPatchField Declaration + Class regionCouplingFvsPatchField Declaration \*---------------------------------------------------------------------------*/ template -class regionCoupleFvsPatchField +class regionCouplingFvsPatchField : public coupledFvsPatchField { // Private data - //- Local reference cast into the regionCouple patch + //- Local reference cast into the regionCoupling patch const regionCoupleFvPatch& regionCouplePatch_; //- Name of remote field to couple to @@ -81,41 +85,40 @@ protected: public: //- Runtime type information -// TypeName(regionCoupleFvPatch::typeName_()); TypeName("regionCoupling"); // Constructors //- Construct from patch and internal field - regionCoupleFvsPatchField + regionCouplingFvsPatchField ( const fvPatch&, const DimensionedField& ); //- Construct from patch, internal field and dictionary - regionCoupleFvsPatchField + regionCouplingFvsPatchField ( const fvPatch&, const DimensionedField&, const dictionary& ); - //- Construct by mapping given regionCoupleFvsPatchField + //- Construct by mapping given regionCouplingFvsPatchField // onto a new patch - regionCoupleFvsPatchField + regionCouplingFvsPatchField ( - const regionCoupleFvsPatchField&, + const regionCouplingFvsPatchField&, const fvPatch&, const DimensionedField&, const fvPatchFieldMapper& ); //- Construct as copy setting internal field reference - regionCoupleFvsPatchField + regionCouplingFvsPatchField ( - const regionCoupleFvsPatchField&, + const regionCouplingFvsPatchField&, const DimensionedField& ); @@ -124,7 +127,7 @@ public: { return tmp > ( - new regionCoupleFvsPatchField(*this) + new regionCouplingFvsPatchField(*this) ); } @@ -136,7 +139,7 @@ public: { return tmp > ( - new regionCoupleFvsPatchField(*this, iF) + new regionCouplingFvsPatchField(*this, iF) ); } @@ -152,7 +155,7 @@ public: } //- Return shadow patch field - const regionCoupleFvsPatchField& shadowPatchField() const; + const regionCouplingFvsPatchField& shadowPatchField() const; //- Write @@ -168,7 +171,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #ifdef NoRepository -# include "regionCoupleFvsPatchField.C" +# include "regionCouplingFvsPatchField.C" #endif // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.C b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.C similarity index 95% rename from src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.C rename to src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.C index a6962c698..ecd0725bb 100644 --- a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFields.C +++ b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.C @@ -30,7 +30,7 @@ Author \*---------------------------------------------------------------------------*/ -#include "regionCoupleFvsPatchFields.H" +#include "regionCouplingFvsPatchFields.H" #include "fvsPatchFields.H" #include "addToRunTimeSelectionTable.H" @@ -41,7 +41,7 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -makeFvsPatchFields(regionCouple); +makeFvsPatchFields(regionCoupling); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.H b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.H similarity index 88% rename from src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.H rename to src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.H index 6a57d58d0..d32d1a879 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchFields.H +++ b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFields.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCouplewFvPatchFields + regionCouplingFvsPatchFields Description Region couple patch field @@ -32,14 +32,14 @@ Author Hrvoje Jasak, Wikki Ltd. All rights reserved SourceFiles - regionCoupleFvPatchFields.C + regionCouplingFvsPatchFields.C \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvPatchFields_H -#define regionCoupleFvPatchFields_H +#ifndef regionCouplingFvsPatchFields_H +#define regionCouplingFvsPatchFields_H -#include "regionCoupleFvPatchField.H" +#include "regionCouplingFvsPatchField.H" #include "fieldTypes.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +49,7 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -makePatchTypeFieldTypedefs(regionCouple) +makeFvsPatchTypeFieldTypedefs(regionCoupling) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFieldsFwd.H b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFieldsFwd.H similarity index 89% rename from src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFieldsFwd.H rename to src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFieldsFwd.H index 9489d7064..385f9da7c 100644 --- a/src/finiteVolume/fields/fvsPatchFields/constraint/regionCouple/regionCoupleFvsPatchFieldsFwd.H +++ b/src/finiteVolume/fields/fvsPatchFields/constraint/regionCoupling/regionCouplingFvsPatchFieldsFwd.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - regionCoupleFvsPatchField + regionCouplingFvsPatchField Description Region couple patch field @@ -33,8 +33,8 @@ Author \*---------------------------------------------------------------------------*/ -#ifndef regionCoupleFvsPatchFieldsFwd_H -#define regionCoupleFvsPatchFieldsFwd_H +#ifndef regionCouplingFvsPatchFieldsFwd_H +#define regionCouplingFvsPatchFieldsFwd_H #include "fvPatchField.H" #include "fieldTypes.H" @@ -46,9 +46,9 @@ namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -template class regionCoupleFvsPatchField; +template class regionCouplingFvsPatchField; -makeFvsPatchTypeFieldTypedefs(regionCouple) +makeFvsPatchTypeFieldTypedefs(regionCoupling) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C index 5881daa64..6fba62205 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C @@ -135,6 +135,7 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const } +// Make patch face non-orthogonality correction vectors void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const { // Non-orthogonality correction on a ggi interface @@ -149,14 +150,6 @@ void Foam::ggiFvPatch::makeCorrVecs(vectorField& cv) const } -const Foam::ggiFvPatch& Foam::ggiFvPatch::shadow() const -{ - const fvPatch& p = this->boundaryMesh()[ggiPolyPatch_.shadowIndex()]; - - return refCast(p); -} - - // Return delta (P to N) vectors across coupled patch Foam::tmp Foam::ggiFvPatch::delta() const { @@ -192,6 +185,14 @@ Foam::tmp Foam::ggiFvPatch::delta() const } +const Foam::ggiFvPatch& Foam::ggiFvPatch::shadow() const +{ + const fvPatch& p = this->boundaryMesh()[ggiPolyPatch_.shadowIndex()]; + + return refCast(p); +} + + bool Foam::ggiFvPatch::master() const { return ggiPolyPatch_.master(); diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C index 57e7df5b5..c67e37eca 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C @@ -45,24 +45,57 @@ namespace Foam } +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::regionCoupleFvPatch::~regionCoupleFvPatch() +{} + + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Make patch weighting factors void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const { - if (rcPolyPatch_.attached()) + if (rcPolyPatch_.coupled()) { - vectorField n = nf(); + if (rcPolyPatch_.master()) + { + vectorField n = nf(); - // Note: mag in the dot-product. - // For all valid meshes, the non-orthogonality will be less that - // 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 & (rcPolyPatch_.reconFaceCellCentres() - Cf())); + // Note: mag in the dot-product. + // For all valid meshes, the non-orthogonality will be less that + // 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 & (rcPolyPatch_.reconFaceCellCentres() - Cf())); - w = nfc/(mag(n & (Cf() - Cn())) + nfc); + 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); + } + } + else + { + // Pick up weights from the master side + scalarField masterWeights(shadow().size()); + shadow().makeWeights(masterWeights); + + scalarField oneMinusW = 1 - masterWeights; + + w = interpolate(oneMinusW); + + 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); + } + } } else { @@ -74,12 +107,35 @@ void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const // Make patch face - neighbour cell distances void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const { - if (rcPolyPatch_.attached()) + if (rcPolyPatch_.coupled()) { - // Stabilised form for bad meshes. HJ, 24/Aug/2011 - vectorField d = delta(); + if (rcPolyPatch_.master()) + { + // Stabilised form for bad meshes. HJ, 24/Aug/2011 + vectorField d = delta(); - dc = 1.0/max(nf() & d, 0.05*mag(d)); + dc = 1.0/max(nf() & d, 0.05*mag(d)); + + if (bridgeOverlap()) + { + scalarField bridgeDeltas = nf() & fvPatch::delta(); + + bridge(bridgeDeltas, dc); + } + } + else + { + scalarField masterDeltas(shadow().size()); + shadow().makeDeltaCoeffs(masterDeltas); + dc = interpolate(masterDeltas); + + if (bridgeOverlap()) + { + scalarField bridgeDeltas = nf() & fvPatch::delta(); + + bridge(bridgeDeltas, dc); + } + } } else { @@ -91,7 +147,7 @@ void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const // Make patch face non-orthogonality correction vectors void Foam::regionCoupleFvPatch::makeCorrVecs(vectorField& cv) const { - if (rcPolyPatch_.attached()) + if (rcPolyPatch_.coupled()) { // Non-orthogonality correction in attached state identical to ggi // interface @@ -111,7 +167,48 @@ void Foam::regionCoupleFvPatch::makeCorrVecs(vectorField& cv) const } -// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +// Return delta (P to N) vectors across coupled patch +Foam::tmp Foam::regionCoupleFvPatch::delta() const +{ + if (rcPolyPatch_.coupled()) + { + if (rcPolyPatch_.master()) + { + tmp tDelta = + rcPolyPatch_.reconFaceCellCentres() - Cn(); + + if (bridgeOverlap()) + { + vectorField bridgeDeltas = Cf() - Cn(); + + bridge(bridgeDeltas, tDelta()); + } + + return tDelta; + } + else + { + tmp tDelta = interpolate + ( + shadow().Cn() - rcPolyPatch_.shadow().reconFaceCellCentres() + ); + + if (bridgeOverlap()) + { + vectorField bridgeDeltas = Cf() - Cn(); + + bridge(bridgeDeltas, tDelta()); + } + + return tDelta; + } + } + else + { + return fvPatch::delta(); + } +} + bool Foam::regionCoupleFvPatch::coupled() const { @@ -138,16 +235,74 @@ const Foam::regionCoupleFvPatch& Foam::regionCoupleFvPatch::shadow() const } -// Return delta (P to N) vectors across coupled patch -Foam::tmp Foam::regionCoupleFvPatch::delta() const +bool Foam::regionCoupleFvPatch::master() const { - if (rcPolyPatch_.attached()) + return rcPolyPatch_.master(); +} + + +bool Foam::regionCoupleFvPatch::fineLevel() const +{ + return true; +} + + +Foam::label Foam::regionCoupleFvPatch::shadowIndex() const +{ + return rcPolyPatch_.shadowIndex(); +} + + +const Foam::ggiLduInterface& +Foam::regionCoupleFvPatch::shadowInterface() const +{ + const fvPatch& p = + shadowRegion().boundary()[rcPolyPatch_.shadowIndex()]; + + return refCast(p); +} + + +Foam::label Foam::regionCoupleFvPatch::zoneSize() const +{ + return rcPolyPatch_.zone().size(); +} + + +const Foam::labelList& Foam::regionCoupleFvPatch::zoneAddressing() const +{ + return rcPolyPatch_.zoneAddressing(); +} + + +const Foam::labelListList& Foam::regionCoupleFvPatch::addressing() const +{ + if (rcPolyPatch_.master()) { - return rcPolyPatch_.reconFaceCellCentres() - Cn(); + return rcPolyPatch_.patchToPatch().masterAddr(); } else { - return fvPatch::delta(); + return rcPolyPatch_.patchToPatch().slaveAddr(); + } +} + + +bool Foam::regionCoupleFvPatch::localParallel() const +{ + return rcPolyPatch_.localParallel(); +} + + +const Foam::scalarListList& Foam::regionCoupleFvPatch::weights() const +{ + if (rcPolyPatch_.master()) + { + return rcPolyPatch_.patchToPatch().masterWeights(); + } + else + { + return rcPolyPatch_.patchToPatch().slaveWeights(); } } @@ -167,7 +322,7 @@ void Foam::regionCoupleFvPatch::initTransfer const unallocLabelList& interfaceData ) const { - transferBuffer_ = interfaceData; + labelTransferBuffer_ = interfaceData; } @@ -177,8 +332,8 @@ Foam::tmp Foam::regionCoupleFvPatch::transfer const unallocLabelList& interfaceData ) const { - //HJ Should this be mapped? 22/Jun/2007 - return shadow().transferBuffer(); + + return shadow().labelTransferBuffer(); } @@ -188,7 +343,7 @@ void Foam::regionCoupleFvPatch::initInternalFieldTransfer const unallocLabelList& iF ) const { - transferBuffer_ = patchInternalField(iF); + labelTransferBuffer_ = patchInternalField(iF); } @@ -198,8 +353,7 @@ Foam::tmp Foam::regionCoupleFvPatch::internalFieldTransfer const unallocLabelList& iF ) const { - //HJ Should this be mapped? 22/Jun/2007 - return shadow().transferBuffer(); + return shadow().labelTransferBuffer(); } diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H index 98d8bae53..c79c0e307 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H @@ -41,7 +41,7 @@ SourceFiles #define regionCoupleFvPatch_H #include "coupledFvPatch.H" -#include "regionCoupleLduInterface.H" +#include "ggiLduInterface.H" #include "regionCouplePolyPatch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -58,16 +58,13 @@ class fvMesh; class regionCoupleFvPatch : public coupledFvPatch, - public regionCoupleLduInterface + public ggiLduInterface { // Private Data //- Reference to polyPatch const regionCouplePolyPatch& rcPolyPatch_; - //- Transfer buffer - mutable labelField transferBuffer_; - protected: @@ -82,12 +79,6 @@ protected: //- Make patch face non-orthogonality correction vectors virtual void makeCorrVecs(vectorField&) const; - //- Return contents of transfer buffer - const labelField& transferBuffer() const - { - return transferBuffer_; - } - public: @@ -101,53 +92,108 @@ public: regionCoupleFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm) : coupledFvPatch(patch, bm), - rcPolyPatch_(refCast(patch)), - transferBuffer_() + rcPolyPatch_(refCast(patch)) {} // Destructor - virtual ~regionCoupleFvPatch() - {} + virtual ~regionCoupleFvPatch(); // Member functions - //- Return true if coupled - virtual bool coupled() const; + // Access - //- Return shadow patch index - int shadowIndex() const - { - return rcPolyPatch_.shadowIndex(); - } + //- Return true if coupled + virtual bool coupled() const; - //- Return shadow region - const fvMesh& shadowRegion() const; + //- Return shadow patch + const fvMesh& shadowRegion() const; - //- Return shadow patch - const regionCoupleFvPatch& shadow() const; + //- Return shadow patch + const regionCoupleFvPatch& shadow() const; - //- Interpolate face field - template - tmp > interpolate(const Field& pf) const - { - return rcPolyPatch_.interpolate(pf); - } + //- Use bridging to fix overlap error in interpolation + bool bridgeOverlap() const + { + return rcPolyPatch_.bridgeOverlap(); + } - template - tmp > interpolate(const tmp >& tpf) const - { - return rcPolyPatch_.interpolate(tpf); - } + //- Return delta (P to N) vectors across coupled patch + virtual tmp delta() const; - //- Return delta (P to N) vectors across coupled patch - virtual tmp delta() const; + + // Interpolation + + //- Interpolate face field + template + tmp > interpolate(const Field& pf) const + { + return rcPolyPatch_.interpolate(pf); + } + + template + tmp > interpolate(const tmp >& tpf) const + { + return rcPolyPatch_.interpolate(tpf); + } + + //- Bridge interpolated face field for uncovered faces + template + void bridge + ( + const Field& bridgeField, + Field& ff + ) const + { + return rcPolyPatch_.bridge(bridgeField, ff); + } // Interface transfer functions + //- Is this the master side? + virtual bool master() const; + + //- Is this the fine level? + virtual bool fineLevel() const; + + //- Return shadow patch index + virtual label shadowIndex() const; + + //- Return shadow interface + virtual const ggiLduInterface& shadowInterface() const; + + //- Return zone size + virtual label zoneSize() const; + + //- Return zone addressing + virtual const labelList& zoneAddressing() const; + + //- Return addressing. Master side returns own addressing and + // slave side returns addressing from master + virtual const labelListList& addressing() const; + + //- Is the patch localised on a single processor + virtual bool localParallel() const; + + //- Return weights. Master side returns own weights and + // slave side returns weights from master + virtual const scalarListList& weights() const; + + //- Return face transformation tensor + virtual const tensorField& forwardT() const + { + return coupledFvPatch::forwardT(); + } + + //- Return neighbour-cell transformation tensor + virtual const tensorField& reverseT() const + { + return coupledFvPatch::reverseT(); + } + //- Return the values of the given internal data adjacent to // the interface as a field virtual tmp interfaceInternalField diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H index e4c201a3a..2fcb18aa4 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H @@ -42,6 +42,7 @@ SourceFiles #include "surfaceInterpolationScheme.H" #include "volFields.H" #include "surfaceFields.H" +#include "magLongDelta.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -124,6 +125,12 @@ public: ) ); + const surfaceScalarField& deltaCoeffs = this->mesh().deltaCoeffs(); + const surfaceScalarField& weights = this->mesh().weights(); + + const magLongDelta& mld = magLongDelta::New(this->mesh()); + const surfaceScalarField& longDelta = mld.magDelta(); + surfaceScalarField& w = tw(); const unallocLabelList& owner = this->mesh().owner(); @@ -133,18 +140,30 @@ public: scalarField& wIn = w.internalField(); + // Larger small for complex arithmetic accuracy + const scalar kSmall = 1000*SMALL; + // Calculate internal weights using field magnitude forAll (owner, faceI) { - scalar mOwn = magPhi[owner[faceI]]; - scalar mNei = magPhi[neighbour[faceI]]; + label own = owner[faceI]; + label nei = neighbour[faceI]; - scalar den = mOwn - mNei; + scalar mOwn = magPhi[own]/(1 - weights[faceI]); + scalar mNei = magPhi[nei]/weights[faceI]; - if (mag(den) > SMALL) + scalar den = magPhi[nei] - magPhi[own]; + + scalar mean = mOwn*mNei/ + ((mOwn + mNei)*longDelta[faceI]*deltaCoeffs[faceI]); + + // Note: complex arithmetic requires extra accuracy + // This is a division of two close subtractions + // HJ, 28/Sep/2011 + if (mag(den) > kSmall) { - scalar mean = 2.0*mOwn*mNei/(mOwn + mNei); - wIn[faceI] = (mean - mNei)/den; + // mOwn > mNei + wIn[faceI] = (magPhi[nei] - mean)/den; } else { @@ -156,29 +175,45 @@ public: { fvsPatchScalarField& wp = w.boundaryField()[pi]; - const fvPatchField& patchPhi = phi.boundaryField()[pi]; + const fvPatchField& pPhi = phi.boundaryField()[pi]; - if (patchPhi.coupled()) + if (pPhi.coupled()) { - scalarField magPhiOwn = mag(patchPhi.patchInternalField()); - scalarField magPhiNei = mag(patchPhi.patchNeighbourField()); + const scalarField& pWeights = weights.boundaryField()[pi]; + + scalarField magPhiOwn = mag(pPhi.patchInternalField()); + scalarField magPhiNei = mag(pPhi.patchNeighbourField()); + + const scalarField& pDeltaCoeffs = + deltaCoeffs.boundaryField()[pi]; + + const scalarField& pLongDelta = mld.magDelta(pi); // Calculate internal weights using field magnitude - forAll (patchPhi, faceI) + forAll (pPhi, faceI) { - scalar mOwn = magPhiOwn[faceI]; - scalar mNei = magPhiNei[faceI]; + scalar mOwn = magPhiOwn[faceI]/(1 - pWeights[faceI]); + scalar mNei = magPhiNei[faceI]/pWeights[faceI]; - scalar den = mOwn - mNei; + scalar den = magPhiNei[faceI] - magPhiOwn[faceI]; - if (mag(den) > SMALL) + // Note: complex arithmetic requires extra accuracy + // This is a division of two close subtractions + // HJ, 28/Sep/2011 + if (mag(den) > kSmall) { - scalar mean = 2.0*mOwn*mNei/(mOwn + mNei); - wp[faceI] = (mean - mNei)/den; + scalar mean = mOwn*mNei/ + ( + (mOwn + mNei)* + pLongDelta[faceI]* + pDeltaCoeffs[faceI] + ); + + wp[faceI] = (magPhiNei[faceI] - mean)/den; } else { - wp[faceI] = 0.5; + wp[faceI] = 0.5; } } } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.C new file mode 100644 index 000000000..7cf447bdf --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.C @@ -0,0 +1,246 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + +\*---------------------------------------------------------------------------*/ + +#include "magLongDelta.H" +#include "surfaceFields.H" +#include "volFields.H" +#include "mapPolyMesh.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(magLongDelta, 0); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::magLongDelta::magLongDelta(const fvMesh& mesh) +: + MeshObject(mesh), + magLongDeltaPtr_(NULL), + magLongDeltaBnd_(mesh.boundary().size(), NULL) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::magLongDelta::~magLongDelta() +{ + clearData(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::magLongDelta::clearData() const +{ + if (magLongDeltaPtr_) + { + delete magLongDeltaPtr_; + } + else + { + forAll (magLongDeltaBnd_, i) + { + if (magLongDeltaBnd_[i]) + { + delete magLongDeltaBnd_[i]; + } + } + } +} + + +void Foam::magLongDelta::makeMagLongDistance() const +{ + if (debug) + { + Info<< "magLongDelta::makeMagLongDistance() :" + << "Constructing magnitude of long cell distance" + << endl; + } + + magLongDeltaPtr_ = new surfaceScalarField + ( + IOobject + ( + "magLongDelta", + mesh().pointsInstance(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh(), + dimensionedScalar("zero", dimLength, 0.0) + ); + surfaceScalarField& mldp = *magLongDeltaPtr_; + + // Set local references to mesh data + const unallocLabelList& owner = mesh().owner(); + const unallocLabelList& neighbour = mesh().neighbour(); + + const vectorField& Cf = mesh().faceCentres(); + const vectorField& C = mesh().cellCentres(); + const vectorField& Sf = mesh().faceAreas(); + const scalarField& magSf = mesh().magSf(); + + forAll (owner, facei) + { + // This must be the same as in surfaceInterpolation.C + scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]])); + scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei])); + mldp[facei] = (SfdOwn + SfdNei)/magSf[facei]; + } + + forAll (mldp.boundaryField(), patchi) + { + const fvPatch& p = mesh().boundary()[patchi]; + + if (p.coupled()) + { + if (!magLongDeltaBnd_[patchi]) + { + makeMagLongDistance(patchi); + } + + mldp.boundaryField()[patchi] = *magLongDeltaBnd_[patchi]; + delete magLongDeltaBnd_[patchi]; + magLongDeltaBnd_[patchi] = &mldp.boundaryField()[patchi]; + } + else + { + delete magLongDeltaBnd_[patchi]; + } + } + + if (debug) + { + Info<< "magLongDelta::makeMagLongDistance() :" + << "Finished magnitude of long cell distance" + << endl; + } +} + + +void Foam::magLongDelta::makeMagLongDistance(label patchi) const +{ + if (debug) + { + Info<< "magLongDelta::makeMagLongDistanceBnd(label patchi) :" + << "Constructing magnitude of long cell distance" + << endl; + } + + const fvPatch& p = mesh().boundary()[patchi]; + + vectorField d = p.fvPatch::delta(); + + magLongDeltaBnd_[patchi] = + new scalarField + ( + (mag(p.Sf() & d) + mag(p.Sf() & (p.delta() - d)))/p.magSf() + ); + + if (debug) + { + Info<< "magLongDelta::makeMagLongDistanceBnd(label patchi) :" + << "Finished magnitude of long cell distance" + << endl; + } +} + + +const Foam::surfaceScalarField& Foam::magLongDelta::magDelta() const +{ + if (!magLongDeltaPtr_) + { + makeMagLongDistance(); + } + + return *magLongDeltaPtr_; +} + + +const Foam::scalarField& Foam::magLongDelta::magDelta +( + const label patchi +) const +{ + if (!mesh().boundary()[patchi].coupled()) + { + FatalErrorIn + ( + "const Foam::scalarField& Foam::magLongDelta::magDelta(" + "const label patchi) const" + ) << "patch is not a coupled. Cannot calculate long distance" + << abort(FatalError); + } + + if (!magLongDeltaBnd_[patchi]) + { + makeMagLongDistance(patchi); + } + + return *magLongDeltaBnd_[patchi]; +} + + +bool Foam::magLongDelta::movePoints() const +{ + if (debug) + { + InfoIn("bool magLongDelta::movePoints() const") + << "Clearing long cell distance data" << endl; + } + + clearData(); + + magLongDeltaBnd_.setSize(mesh().boundary().size(), NULL); + + return true; +} + + +bool Foam::magLongDelta::updateMesh(const mapPolyMesh& mpm) const +{ + if (debug) + { + InfoIn("bool magLongDelta::updateMesh(const mapPolyMesh&) const") + << "Clearing long cell distance data" << endl; + } + + clearData(); + + magLongDeltaBnd_.setSize(mesh().boundary().size(), NULL); + + return true; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.H new file mode 100644 index 000000000..bddfb5146 --- /dev/null +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/magLongDelta.H @@ -0,0 +1,124 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::magLongDelta + +Description + Long distance to the neighbour via face centre + +SourceFiles + magLongDelta.C + +\*---------------------------------------------------------------------------*/ + +#ifndef magLongDelta_H +#define magLongDelta_H + +#include "MeshObject.H" +#include "fvMesh.H" +#include "primitiveFieldsFwd.H" +#include "surfaceFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class mapPolyMesh; + +/*---------------------------------------------------------------------------*\ + Class magLongDelta Declaration +\*---------------------------------------------------------------------------*/ + +class magLongDelta +: + public MeshObject +{ + // Private data + + //- Magnitude of the normal distances to the face + mutable surfaceScalarField* magLongDeltaPtr_; + + //- Magnitude of the normal distances on coupled patches + //- Data in magLongDeltaPtr_ is reused, if it exists + mutable List magLongDeltaBnd_; + + + // Private member functions + + //- Construct Long distance + void makeMagLongDistance() const; + + //- Construct Long distance + void makeMagLongDistance(label patchi) const; + + //- Clear Data + void clearData() const; + +public: + + // Declare name of the class and its debug switch + TypeName("magLongDelta"); + + + // Constructors + + //- Construct given an fvMesh + explicit magLongDelta(const fvMesh&); + + + // Destructor + + virtual ~magLongDelta(); + + + // Member functions + + //- Return reference to the magnitude of the long distance + const surfaceScalarField& magDelta() const; + + //- Return reference to the magnitude of the long distance + //- for patch i + const scalarField& magDelta(const label patchi) const; + + //- Update after mesh motion: + // Delete data when the mesh moves + virtual bool movePoints() const; + + //- Update after topo change: + // Delete data when mesh changes + virtual bool updateMesh(const mapPolyMesh&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //