diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C index e306ccf23..5ff86f070 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/regionCouple/regionCoupleFvPatchField.C @@ -180,25 +180,27 @@ void regionCoupleFvPatchField::evaluate Field fOwn = this->patchInternalField(); Field fNei = this->patchNeighbourField(); - scalarField mOwn = mag(fOwn); - scalarField mNei = mag(fNei); - scalarField mean = 2*(mOwn*mNei)/(mOwn + mNei + SMALL); + scalarField magFOwn = mag(fOwn); + scalarField magFNei = mag(fNei); - scalarField weights(fOwn.size(), 0.5); - - scalar den; + // Calculate internal weights using field magnitude + scalarField weights(fOwn.size()); forAll (weights, faceI) { - den = mOwn[faceI] - mNei[faceI]; + scalar mOwn = magFOwn[faceI]; + scalar mNei = magFNei[faceI]; + + scalar den = mOwn - mNei; if (mag(den) > SMALL) { - weights[faceI] = (mean[faceI] - mNei[faceI])/den; + scalar mean = 2.0*mOwn*mNei/(mOwn + mNei); + weights[faceI] = (mean - mNei)/den; } else { - // Use 0.5 weights + weights[faceI] = 0.5; } } diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H index b7d2e5df5..e4c201a3a 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H @@ -131,33 +131,63 @@ public: scalarField magPhi = mag(phi); - // Initialise weights to 0.5 for uniform field (den = 0) scalarField& wIn = w.internalField(); - wIn = 0.5; // Calculate internal weights using field magnitude - scalar mOwn, mNei, den, mean; - forAll (owner, faceI) { - mOwn = magPhi[owner[faceI]]; - mNei = magPhi[neighbour[faceI]]; + scalar mOwn = magPhi[owner[faceI]]; + scalar mNei = magPhi[neighbour[faceI]]; - mean = 2*(mOwn*mNei)/(mOwn + mNei + SMALL); - den = mOwn - mNei; + scalar den = mOwn - mNei; if (mag(den) > SMALL) { + scalar mean = 2.0*mOwn*mNei/(mOwn + mNei); wIn[faceI] = (mean - mNei)/den; } else { - // Use 0.5 weights + wIn[faceI] = 0.5; } } - // Boundary weights are 1 - w.boundaryField() = 1; + forAll (phi.boundaryField(), pi) + { + fvsPatchScalarField& wp = w.boundaryField()[pi]; + + const fvPatchField& patchPhi = phi.boundaryField()[pi]; + + if (patchPhi.coupled()) + { + scalarField magPhiOwn = mag(patchPhi.patchInternalField()); + scalarField magPhiNei = mag(patchPhi.patchNeighbourField()); + + // Calculate internal weights using field magnitude + forAll (patchPhi, faceI) + { + scalar mOwn = magPhiOwn[faceI]; + scalar mNei = magPhiNei[faceI]; + + scalar den = mOwn - mNei; + + if (mag(den) > SMALL) + { + scalar mean = 2.0*mOwn*mNei/(mOwn + mNei); + wp[faceI] = (mean - mNei)/den; + } + else + { + wp[faceI] = 0.5; + } + } + } + else + { + // Boundary weights for uncoupled patches are 1 + wp = 1; + } + } return tw; }