diff --git a/.gitignore b/.gitignore index 9ac47f062..af81b4b2f 100644 --- a/.gitignore +++ b/.gitignore @@ -101,4 +101,12 @@ src/lduSolvers/amg/amgPolicy/samgPolicy.H src/lduSolvers/amg/amgPolicy/aamgPolicy.C src/lduSolvers/amg/amgPolicy/aamgPolicy.H +# The following files are blacklisted because of a DMCA complaint by ANSYS. +src/lduSolvers/tools/PriorityArray.C +src/lduSolvers/tools/PriorityArray.H +src/lduSolvers/amg/amgPolicy/samgPolicy.C +src/lduSolvers/amg/amgPolicy/samgPolicy.H +src/lduSolvers/amg/amgPolicy/aamgPolicy.C +src/lduSolvers/amg/amgPolicy/aamgPolicy.H + # end-of-file diff --git a/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H b/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H index 3412e3ebb..ce028cfd0 100644 --- a/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H +++ b/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H @@ -95,7 +95,7 @@ public: ClassName("GGIInterpolation"); //- Quick reject names - static const NamedEnum quickRejectNames_; + static const NamedEnum quickRejectNames_; // Constructors diff --git a/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationName.C b/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationName.C index 0a4b75c92..3b06ad428 100644 --- a/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationName.C +++ b/src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationName.C @@ -41,12 +41,15 @@ defineTypeNameAndDebug(Foam::GGIInterpolationName, 0); template<> const char* -Foam::NamedEnum::names[] = +Foam::NamedEnum::names[] = { - "3DDistance", + "distance3D", "AABB", + "bbOctree", "nSquared" }; +const Foam::NamedEnum + Foam::GGIInterpolationName::quickRejectNames_; // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C index 2d8c79231..1e4612ace 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/BlockLduMatrixATmul.C @@ -42,7 +42,7 @@ void Foam::BlockLduMatrix::Amul // Initialise the update of coupled interfaces initInterfaces(coupleUpper_, Ax, x); - + AmulCore(Ax, x); // Update coupled interfaces @@ -64,14 +64,43 @@ void Foam::BlockLduMatrix::AmulCore const unallocLabelList& u = lduAddr().upperAddr(); const unallocLabelList& l = lduAddr().lowerAddr(); + const TypeCoeffField& Diag = this->diag(); const TypeCoeffField& Upper = this->upper(); - // Diagonal multiplication, no indirection - multiply(Ax, Diag, x); - // Create multiplication function object typename BlockCoeff::multiply mult; + + // AmulCore must be additive to account for initialisation step + // in ldu interfaces. HJ, 6/Nov/2007 + // Fixed by IC 19/Oct/2011 + if (Diag.activeType() == blockCoeffBase::SCALAR) + { + const scalarTypeField& activeDiag = Diag.asScalar(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Ax[cellI] += mult(activeDiag[cellI], x[cellI]); + } + } + else if (Diag.activeType() == blockCoeffBase::LINEAR) + { + const linearTypeField& activeDiag = Diag.asLinear(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Ax[cellI] += mult(activeDiag[cellI], x[cellI]); + } + } + else if (Diag.activeType() == blockCoeffBase::SQUARE) + { + const squareTypeField& activeDiag = Diag.asSquare(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Ax[cellI] += mult(activeDiag[cellI], x[cellI]); + } + } // Lower multiplication @@ -209,12 +238,40 @@ void Foam::BlockLduMatrix::TmulCore const TypeCoeffField& Diag = this->diag(); const TypeCoeffField& Upper = this->upper(); - // Diagonal multiplication, no indirection - multiply(Tx, Diag, x); - // Create multiplication function object typename BlockCoeff::multiply mult; + // AmulCore must be additive to account for initialisation step + // in ldu interfaces. HJ, 6/Nov/2007 + // Fixed by IC 19/Oct/2011 + if (Diag.activeType() == blockCoeffBase::SCALAR) + { + const scalarTypeField& activeDiag = Diag.asScalar(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Tx[cellI] += mult(activeDiag[cellI], x[cellI]); + } + } + else if (Diag.activeType() == blockCoeffBase::LINEAR) + { + const linearTypeField& activeDiag = Diag.asLinear(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Tx[cellI] += mult(activeDiag[cellI], x[cellI]); + } + } + else if (Diag.activeType() == blockCoeffBase::SQUARE) + { + const squareTypeField& activeDiag = Diag.asSquare(); + + for (register label cellI = 0; cellI < u.size(); cellI++) + { + Tx[cellI] += mult(activeDiag[cellI].T(), x[cellI]); + } + } + // Upper multiplication if (Upper.activeType() == blockCoeffBase::SCALAR) diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C index 0f662410e..34703cc24 100644 --- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C +++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduMatrix/scalarBlockLduMatrix.C @@ -248,37 +248,53 @@ void BlockLduMatrix::AmulCore const scalarField& x ) const { - const scalarField& Diag = diag(); - - for (register label rowI = 0; rowI < x.size(); rowI++) - { - Ax[rowI] = Diag[rowI]*x[rowI]; - } - // Note: pointer looping const label* const __restrict__ U = lduAddr().upperAddr().begin(); const label* const __restrict__ L = lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ X = x.begin(); + scalar* __restrict__ AX = Ax.begin(); + + if (thereIsDiag()) + { + const scalar* const __restrict__ diagPtr = diag().begin(); + register const label nCells = diag().size(); + for (register label cell=0; cell::TmulCore const scalarField& x ) const { - const scalarField& Diag = diag(); - - for (register label rowI = 0; rowI < x.size(); rowI++) - { - Tx[rowI] = Diag[rowI]*x[rowI]; - } - // Note: pointer looping const label* const __restrict__ U = lduAddr().upperAddr().begin(); const label* const __restrict__ L = lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ X = x.begin(); + scalar* __restrict__ TX = Tx.begin(); + + if (thereIsDiag()) + { + const scalar* const __restrict__ diagPtr = diag().begin(); + + register const label nCells = diag().size(); + for (register label cell=0; cell BlockLduMatrix::H(const scalarField& x) const const label* const __restrict__ U = lduAddr().upperAddr().begin(); const label* const __restrict__ L = lduAddr().lowerAddr().begin(); - const scalar* const __restrict__ Lower = lower().begin(); - const scalar* const __restrict__ Upper = upper().begin(); const scalar* const __restrict__ X = x.begin(); - - scalar* __restrict__ R = result.begin(); - - for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + + if (symmetric()) { - R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]]; + if (thereIsUpper()) + { + const scalar* const __restrict__ Upper = upper().begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]]; + } + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[L[coeffI]] -= Upper[coeffI]*X[L[coeffI]]; + } + } + else if (thereIsLower()) + { + const scalar* const __restrict__ Lower = lower().begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[U[coeffI]] -= Lower[coeffI]*X[U[coeffI]]; + } + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]]; + } + } } - - for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + else { - R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]]; + const scalar* const __restrict__ Lower = lower().begin(); + const scalar* const __restrict__ Upper = upper().begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[U[coeffI]] -= Upper[coeffI]*X[U[coeffI]]; + } + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[L[coeffI]] -= Lower[coeffI]*X[L[coeffI]]; + } } } - + return tresult; } @@ -397,19 +468,57 @@ tmp BlockLduMatrix::faceH(const scalarField& x) const if (thereIsUpper() || thereIsLower()) { - // Note: pointer looping - const label* const __restrict__ U = lduAddr().upperAddr().begin(); - const label* const __restrict__ L = lduAddr().lowerAddr().begin(); - - const scalar* const __restrict__ Lower = lower().begin(); - const scalar* const __restrict__ Upper = upper().begin(); - const scalar* const __restrict__ X = x.begin(); - - scalar* __restrict__ R = result.begin(); - - for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + if (symmetric()) { - R[coeffI] = Upper[coeffI]*X[U[coeffI]] - Lower[coeffI]*X[L[coeffI]]; + if (thereIsUpper()) + { + // Note: pointer looping + const label* const __restrict__ U = lduAddr().upperAddr().begin(); + const label* const __restrict__ L = lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ Upper = upper().begin(); + const scalar* const __restrict__ X = x.begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[coeffI] = Upper[coeffI]*(X[U[coeffI]] - X[L[coeffI]]); + } + } + else if (thereIsLower()) + { + // Note: pointer looping + const label* const __restrict__ U = lduAddr().upperAddr().begin(); + const label* const __restrict__ L = lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ Lower = lower().begin(); + const scalar* const __restrict__ X = x.begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[coeffI] = Lower[coeffI]*(X[U[coeffI]] - X[L[coeffI]]); + } + } + } + else + { + // Note: pointer looping + const label* const __restrict__ U = lduAddr().upperAddr().begin(); + const label* const __restrict__ L = lduAddr().lowerAddr().begin(); + + const scalar* const __restrict__ Lower = lower().begin(); + const scalar* const __restrict__ Upper = upper().begin(); + const scalar* const __restrict__ X = x.begin(); + + scalar* __restrict__ R = result.begin(); + + for (register label coeffI = 0; coeffI < upper().size(); coeffI++) + { + R[coeffI] = Upper[coeffI]*X[U[coeffI]] - Lower[coeffI]*X[L[coeffI]]; + } } } diff --git a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/ggiLduInterfaceField/ggiLduInterfaceField.C b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/ggiLduInterfaceField/ggiLduInterfaceField.C index a641a4541..198a62da4 100644 --- a/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/ggiLduInterfaceField/ggiLduInterfaceField.C +++ b/src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/ggiLduInterfaceField/ggiLduInterfaceField.C @@ -69,4 +69,5 @@ void Foam::ggiLduInterfaceField::transformCoupleField } + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C index a0330b254..24e6ed35b 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.C @@ -257,7 +257,7 @@ const Foam::cyclicGgiPolyPatch& Foam::cyclicGgiPolyPatch::cyclicShadow() const } -void Foam::cyclicGgiPolyPatch::calcTransforms() +void Foam::cyclicGgiPolyPatch::calcTransforms() const { if (active() && debug) { diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.H index 48034ea98..f686557b2 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclicGgi/cyclicGgiPolyPatch.H @@ -75,7 +75,7 @@ class cyclicGgiPolyPatch //- Calculate cyclic transforms (rotation and translation) // Virtual over-ride for base GGI patch. HJ, 14/Jan/2009 - virtual void calcTransforms(); + virtual void calcTransforms() const; //- Check definition: angles and offsets void checkDefinition() const; diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C index 0ba7f4ec9..2fd735adc 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C @@ -202,8 +202,7 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const 0, // Non-overlapping face tolerances 0, // HJ, 24/Oct/2008 true, // Rescale weighting factors. Bug fix, MB. -// ggiInterpolation::AABB - ggiInterpolation::BB_OCTREE // Octree search, MB. + reject_ // Quick rejection algorithm, default BB_OCTREE ); // Abort immediately if uncovered faces are present and the option @@ -469,6 +468,7 @@ Foam::ggiPolyPatch::ggiPolyPatch shadowName_("initializeMe"), zoneName_("initializeMe"), bridgeOverlap_(false), + reject_(ggiZoneInterpolation::BB_OCTREE), shadowIndex_(-1), zoneIndex_(-1), patchToPatchPtr_(NULL), @@ -490,13 +490,15 @@ Foam::ggiPolyPatch::ggiPolyPatch const polyBoundaryMesh& bm, const word& shadowName, const word& zoneName, - const bool bridgeOverlap + const bool bridgeOverlap, + const ggiZoneInterpolation::quickReject reject ) : coupledPolyPatch(name, size, start, index, bm), shadowName_(shadowName), zoneName_(zoneName), bridgeOverlap_(bridgeOverlap), + reject_(reject), shadowIndex_(-1), zoneIndex_(-1), patchToPatchPtr_(NULL), @@ -521,6 +523,7 @@ Foam::ggiPolyPatch::ggiPolyPatch shadowName_(dict.lookup("shadowPatch")), zoneName_(dict.lookup("zone")), bridgeOverlap_(dict.lookup("bridgeOverlap")), + reject_(ggiZoneInterpolation::BB_OCTREE), shadowIndex_(-1), zoneIndex_(-1), patchToPatchPtr_(NULL), @@ -530,7 +533,15 @@ Foam::ggiPolyPatch::ggiPolyPatch localParallelPtr_(NULL), receiveAddrPtr_(NULL), sendAddrPtr_(NULL) -{} +{ + if (dict.found("quickReject")) + { + reject_ = ggiZoneInterpolation::quickRejectNames_.read + ( + dict.lookup("quickReject") + ); + } +} Foam::ggiPolyPatch::ggiPolyPatch @@ -543,6 +554,7 @@ Foam::ggiPolyPatch::ggiPolyPatch shadowName_(pp.shadowName_), zoneName_(pp.zoneName_), bridgeOverlap_(pp.bridgeOverlap_), + reject_(pp.reject_), shadowIndex_(-1), zoneIndex_(-1), patchToPatchPtr_(NULL), @@ -569,6 +581,7 @@ Foam::ggiPolyPatch::ggiPolyPatch shadowName_(pp.shadowName_), zoneName_(pp.zoneName_), bridgeOverlap_(pp.bridgeOverlap_), + reject_(pp.reject_), shadowIndex_(-1), zoneIndex_(-1), patchToPatchPtr_(NULL), @@ -741,6 +754,11 @@ void Foam::ggiPolyPatch::initAddressing() { // Calculate transforms for correct GGI cut calcTransforms(); + + if (master()) + { + shadow().calcTransforms(); + } // Force zone addressing and remote zone addressing // (uses GGI interpolator) @@ -847,7 +865,7 @@ void Foam::ggiPolyPatch::updateMesh() } -void Foam::ggiPolyPatch::calcTransforms() +void Foam::ggiPolyPatch::calcTransforms() const { // Simplest mixing plane: no transform or separation. HJ, 24/Oct/2008 forwardT_.setSize(0); diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H index d4a9613cf..0576e4f46 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H @@ -73,6 +73,10 @@ class ggiPolyPatch //- Use bridging to fix overlap error in interpolation Switch bridgeOverlap_; + //- Quick reject algorithm + ggiZoneInterpolation::quickReject reject_; + + // Demand-driven data //- Shadow patch index. Delayed evaluation for construction @@ -214,7 +218,9 @@ public: const polyBoundaryMesh& bm, const word& shadowName, const word& zoneName, - const bool bridgeOverlap + const bool bridgeOverlap, + const ggiZoneInterpolation::quickReject + reject = ggiZoneInterpolation::BB_OCTREE ); //- Construct from dictionary @@ -307,6 +313,11 @@ public: //- Return interpolation face zone const faceZone& zone() const; + //- Is this the slave side? + bool slave() const + { + return !master(); + } // Interpolation data diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C index aac0804e3..98625a86c 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.C @@ -289,6 +289,35 @@ void cyclicGgiFvPatchField::updateInterfaceMatrix {} +template +void cyclicGgiFvPatchField::initInterfaceMatrixUpdate +( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType +) const +{ + notImplemented("void cyclicGgiFvPatchField::initInterfaceMatrixUpdate" + " for block-coupled matrices of tensor types" + ) +} + + +template +void cyclicGgiFvPatchField::updateInterfaceMatrix +( + const Field&, + Field&, + const BlockLduMatrix&, + const CoeffField&, + const Pstream::commsTypes commsType +) const +{ +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.H index c24b40d2a..87f9e762b 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchField.H @@ -168,7 +168,26 @@ public: const bool switchToLhs ) const; + //- Initialise neighbour matrix update + virtual void initInterfaceMatrixUpdate + ( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType + ) const; + //- Update result field based on interface functionality + virtual void updateInterfaceMatrix + ( + const Field&, + Field&, + const BlockLduMatrix&, + const CoeffField&, + const Pstream::commsTypes commsType + ) const; + //- GGI coupled interface functions //- Does the patch field perform the transfromation diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.C index 37cf372c4..8bb29e950 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.C @@ -37,12 +37,102 @@ Contributor: #include "cyclicGgiFvPatchFields.H" #include "addToRunTimeSelectionTable.H" #include "volFields.H" +#include "scalarCoeffField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +void cyclicGgiFvPatchField::initInterfaceMatrixUpdate +( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType +) const +{ + // Communication is allowed either before or after processor + // patch comms. HJ, 11/Jul/2011 + + // Get shadow face-cells and assemble shadow field + const unallocLabelList& sfc = cyclicGgiPatch_.shadow().faceCells(); + + scalarField sField(sfc.size()); + + forAll (sField, i) + { + sField[i] = psiInternal[sfc[i]]; + } + + // Note: scalar interpolate does not get a transform, so this is safe + // HJ, 12/Jan/2009 + scalarField pnf = cyclicGgiPatch_.interpolate(sField); + + // Multiply the field by coefficients and add into the result + const unallocLabelList& fc = cyclicGgiPatch_.faceCells(); + + forAll(fc, elemI) + { + result[fc[elemI]] -= coeffs[elemI]*pnf[elemI]; + }; +} + + +template<> +void cyclicGgiFvPatchField::initInterfaceMatrixUpdate +( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType +) const +{ + // Communication is allowed either before or after processor + // patch comms. HJ, 11/Jul/2011 + + // Get shadow face-cells and assemble shadow field + const unallocLabelList& sfc = cyclicGgiPatch_.shadow().faceCells(); + + Field sField(sfc.size()); + + forAll (sField, i) + { + sField[i] = psiInternal[sfc[i]]; + } + + // Transformation is handled in interpolation. HJ, 7/Jan/2009 + Field pnf = cyclicGgiPatch_.interpolate(sField); + + if (coeffs.activeType() == blockCoeffBase::SCALAR) + { + pnf = coeffs.asScalar() * pnf; + } + else if (coeffs.activeType() == blockCoeffBase::LINEAR) + { + pnf = cmptMultiply(coeffs.asLinear(), pnf); + } + else if (coeffs.activeType() == blockCoeffBase::SQUARE) + { + pnf = coeffs.asSquare() & pnf; + } + + // Multiply the field by coefficients and add into the result + const unallocLabelList& fc = cyclicGgiPatch_.faceCells(); + + // Multiply the field by coefficients and add into the result + forAll(fc, elemI) + { + result[fc[elemI]] -= pnf[elemI]; + } +} + + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // makePatchFields(cyclicGgi); diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.H index bd4c18e7f..a4759a6d6 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicGgi/cyclicGgiFvPatchFields.H @@ -51,6 +51,30 @@ SourceFiles namespace Foam { +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +void cyclicGgiFvPatchField::initInterfaceMatrixUpdate +( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType +) const; + + +template<> +void cyclicGgiFvPatchField::initInterfaceMatrixUpdate +( + const Field& psiInternal, + Field& result, + const BlockLduMatrix& m, + const CoeffField& coeffs, + const Pstream::commsTypes commsType +) const; + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // makePatchTypeFieldTypedefs(cyclicGgi)