From 7fb371eb031ada1e8c2d80a79a5b6aef63d18c70 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 4 Jan 2018 15:17:56 +0000 Subject: [PATCH] reactingEulerFoam::phasePair: Added const_iterator which provides access to the current phase and the corresponding other phase for each of the phases in the pair. This allows some simplification of the phase pair loops in several sub-models and avoids the need for pointer swaps. --- .../NonRandomTwoLiquid/NonRandomTwoLiquid.C | 6 +- .../HeatAndMassTransferPhaseSystem.C | 16 ++- .../HeatTransferPhaseSystem.C | 16 ++- ...terfaceCompositionPhaseChangePhaseSystem.C | 4 +- .../MomentumTransferPhaseSystem.C | 29 ++--- .../PopulationBalancePhaseSystem.C | 2 +- .../ThermalPhaseChangePhaseSystem.C | 30 +++-- .../phasePair/phasePair/phasePair.H | 55 ++++++++- .../phasePair/phasePair/phasePairI.H | 115 +++++++++++++++++- .../phasePair/phasePairKey/phasePairKey.C | 2 +- .../populationBalanceModel.C | 4 +- .../reactingMultiphaseEulerFoam/pU/pEqn.H | 2 +- .../reactingMultiphaseEulerFoam/pUf/UEqns.H | 2 +- .../reactingMultiphaseEulerFoam/pUf/pEqn.H | 4 +- 14 files changed, 218 insertions(+), 69 deletions(-) diff --git a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/NonRandomTwoLiquid/NonRandomTwoLiquid.C b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/NonRandomTwoLiquid/NonRandomTwoLiquid.C index f3704e1e39..7f1222f797 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/NonRandomTwoLiquid/NonRandomTwoLiquid.C +++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/interfaceCompositionModels/NonRandomTwoLiquid/NonRandomTwoLiquid.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -216,7 +216,7 @@ Foam::interfaceCompositionModels::NonRandomTwoLiquid::Yf *speciesModel1_->Yf(speciesName, Tf) *gamma1_; } - else if(speciesName == species2Name_) + else if (speciesName == species2Name_) { return this->otherThermo_.composition().Y(speciesName) @@ -248,7 +248,7 @@ YfPrime *speciesModel1_->YfPrime(speciesName, Tf) *gamma1_; } - else if(speciesName == species2Name_) + else if (speciesName == species2Name_) { return this->otherThermo_.composition().Y(speciesName) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C index 8bee082f3d..bbd0a394ab 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -242,9 +242,6 @@ Foam::HeatAndMassTransferPhaseSystem::heatTransfer() const this->phasePairs_[heatTransferModelIter.key()] ); - const phaseModel* phase = &pair.phase1(); - const phaseModel* otherPhase = &pair.phase2(); - const volScalarField& Tf(*Tf_[pair]); const volScalarField K1 @@ -270,14 +267,15 @@ Foam::HeatAndMassTransferPhaseSystem::heatTransfer() const forAllConstIter(phasePair, pair, iter) { - const volScalarField& he(phase->thermo().he()); - volScalarField Cpv(phase->thermo().Cpv()); + const phaseModel& phase = iter(); - *eqns[phase->name()] += - (*K)*(Tf - phase->thermo().T()) + const volScalarField& he(phase.thermo().he()); + volScalarField Cpv(phase.thermo().Cpv()); + + *eqns[phase.name()] += + (*K)*(Tf - phase.thermo().T()) + KEff/Cpv*he - fvm::Sp(KEff/Cpv, he); - Swap(phase, otherPhase); Swap(K, otherK); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C index 3d594234fe..fed441b0fb 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C @@ -146,19 +146,17 @@ Foam::HeatTransferPhaseSystem::heatTransfer() const const phasePair& pair(this->phasePairs_[heatTransferModelIter.key()]); - const phaseModel* phase = &pair.phase1(); - const phaseModel* otherPhase = &pair.phase2(); - forAllConstIter(phasePair, pair, iter) { - const volScalarField& he(phase->thermo().he()); - volScalarField Cpv(phase->thermo().Cpv()); + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); - *eqns[phase->name()] += - K*(otherPhase->thermo().T() - phase->thermo().T() + he/Cpv) + const volScalarField& he(phase.thermo().he()); + volScalarField Cpv(phase.thermo().Cpv()); + + *eqns[phase.name()] += + K*(otherPhase.thermo().T() - phase.thermo().T() + he/Cpv) - fvm::Sp(K/Cpv, he); - - Swap(phase, otherPhase); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C index 64864ab096..a78786fe49 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C @@ -159,7 +159,7 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem::iDmdt { tiDmdt.ref() += this->iDmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) ); } } @@ -208,7 +208,7 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem::dmdt { tDmdt.ref() += this->dmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) ); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 22b46ea69c..7c78c7c5d2 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -554,19 +554,11 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const ) { const volScalarField& K(*KdIter()); - const phasePair& pair(this->phasePairs_[KdIter.key()]); - const phaseModel* phase = &pair.phase1(); - const phaseModel* otherPhase = &pair.phase2(); - forAllConstIter(phasePair, pair, iter) { - const volVectorField& U = phase->U(); - - *eqns[phase->name()] -= fvm::Sp(K, U); - - Swap(phase, otherPhase); + *eqns[iter().name()] -= fvm::Sp(K, iter().U()()); } } @@ -591,28 +583,25 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const ) { const volScalarField& Vm(*VmIter()); - const phasePair& pair(this->phasePairs_[VmIter.key()]); - const phaseModel* phase = &pair.phase1(); - const phaseModel* otherPhase = &pair.phase2(); - forAllConstIter(phasePair, pair, iter) { - const volVectorField& U = phase->U(); - const surfaceScalarField& phi = phase->phi(); + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); - *eqns[phase->name()] -= + const volVectorField& U = phase.U(); + const surfaceScalarField& phi = phase.phi(); + + *eqns[phase.name()] -= Vm *( fvm::ddt(U) + fvm::div(phi, U) - fvm::Sp(fvc::div(phi), U) - - otherPhase->DUDt() + - otherPhase.DUDt() ) - + this->MRF_.DDt(Vm, U - otherPhase->U()); - - Swap(phase, otherPhase); + + this->MRF_.DDt(Vm, U - otherPhase.U()); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C index cf06e62036..9b3030634c 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C @@ -191,7 +191,7 @@ Foam::PopulationBalancePhaseSystem::dmdt { tDmdtFromKey.ref() += this->dmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name(), false) ); } } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C index cc2cdaf528..06de0ed78c 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C @@ -190,7 +190,7 @@ Foam::ThermalPhaseChangePhaseSystem::iDmdt { tiDmdt.ref() += this->iDmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name()) ); } } @@ -252,7 +252,7 @@ Foam::ThermalPhaseChangePhaseSystem::wDmdt { twDmdt.ref() += this->wDmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name()) ); } } @@ -301,7 +301,7 @@ Foam::ThermalPhaseChangePhaseSystem::dmdt { tDmdt.ref() += this->dmdt ( - phasePairKey(phase.name(), pair.other(phase).name(), false) + phasePairKey(phase.name(), pair.otherPhase(phase).name()) ); } } @@ -583,23 +583,24 @@ void Foam::ThermalPhaseChangePhaseSystem::correctThermo() wMDotL *= 0; bool wallBoilingActive = false; - const phaseModel* phase1Ptr = &phase1; - const phaseModel* phase2Ptr = &phase2; forAllConstIter(phasePair, pair, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + if ( - phase1Ptr->mesh().foundObject + phase.mesh().foundObject ( - "alphat." + phase1Ptr->name() + "alphat." + phase.name() ) ) { const volScalarField& alphat = - phase1Ptr->mesh().lookupObject + phase.mesh().lookupObject ( - "alphat." + phase1Ptr->name() + "alphat." + phase.name() ); const fvPatchList& patches = this->mesh().boundary(); @@ -620,11 +621,10 @@ void Foam::ThermalPhaseChangePhaseSystem::correctThermo() ( alphat.boundaryField()[patchi] ); - phasePairKey key - ( - phase1Ptr->name(), phase2Ptr->name(), false - ); - if(PCpatch.activePhasePair(key)) + + phasePairKey key(phase.name(), otherPhase.name()); + + if (PCpatch.activePhasePair(key)) { wallBoilingActive = true; @@ -649,8 +649,6 @@ void Foam::ThermalPhaseChangePhaseSystem::correctThermo() } } } - - Swap(phase1Ptr, phase2Ptr); } if (wallBoilingActive) diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H index 26ef837e8e..423bba0f26 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePair.H @@ -162,10 +162,63 @@ public: //- Return the other phase relative to the given phase // Generates a FatalError if this phasePair does not contain // the given phase - inline const phaseModel& other(const phaseModel& phase) const; + inline const phaseModel& otherPhase(const phaseModel& phase) const; //- Return gravitation acceleration inline const uniformDimensionedVectorField& g() const; + + + //- STL const_iterator + class const_iterator + { + // Private data + + //- Reference to the pair for which this is an iterator + const phasePair& pair_; + + //- Current index + label index_; + + //- Construct an iterator with the given index + inline const_iterator(const phasePair&, const label index); + + public: + + friend class phasePair; + + // Constructors + + //- Construct from pair, moving to its 'begin' position + inline explicit const_iterator(const phasePair&); + + + // Member operators + + inline bool operator==(const const_iterator&) const; + + inline bool operator!=(const const_iterator&) const; + + inline const phaseModel& operator*() const; + inline const phaseModel& operator()() const; + + inline const phaseModel& otherPhase() const; + + inline const_iterator& operator++(); + inline const_iterator operator++(int); + }; + + + //- const_iterator set to the beginning of the pair + inline const_iterator cbegin() const; + + //- const_iterator set to beyond the end of the pair + inline const_iterator cend() const; + + //- const_iterator set to the beginning of the pair + inline const_iterator begin() const; + + //- const_iterator set to beyond the end of the pair + inline const_iterator end() const; }; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H index 6eaf1a612d..5a1373f45e 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePair/phasePairI.H @@ -43,7 +43,7 @@ inline bool Foam::phasePair::contains(const phaseModel& phase) const } -inline const Foam::phaseModel& Foam::phasePair::other +inline const Foam::phaseModel& Foam::phasePair::otherPhase ( const phaseModel& phase ) const @@ -73,4 +73,117 @@ inline const Foam::uniformDimensionedVectorField& Foam::phasePair::g() const } +// * * * * * * * * * * * * * * * * Iterators * * * * * * * * * * * * * * * * // + +inline Foam::phasePair::const_iterator::const_iterator +( + const phasePair& pair, + const label index +) +: + pair_(pair), + index_(index) +{} + + +inline Foam::phasePair::const_iterator::const_iterator(const phasePair& pair) +: + const_iterator(pair, 0) +{} + + +inline bool Foam::phasePair::const_iterator::operator== +( + const const_iterator& iter +) const +{ + return (this->index_ == iter.index_); +} + + +inline bool Foam::phasePair::const_iterator::operator!= +( + const const_iterator& iter +) const +{ + return !(this->operator==(iter)); +} + + +inline const Foam::phaseModel& +Foam::phasePair::const_iterator::operator*() const +{ + if (index_ == 0) + { + return pair_.phase1_; + } + else + { + return pair_.phase2_; + } +} + + +inline const Foam::phaseModel& +Foam::phasePair::const_iterator::operator()() const +{ + return operator*(); +} + + +inline const Foam::phaseModel& +Foam::phasePair::const_iterator::otherPhase() const +{ + if (index_ == 0) + { + return pair_.phase2_; + } + else + { + return pair_.phase1_; + } +} + + +inline Foam::phasePair::const_iterator& +Foam::phasePair::const_iterator::operator++() +{ + index_++; + return *this; +} + + +inline Foam::phasePair::const_iterator +Foam::phasePair::const_iterator::operator++(int) +{ + const_iterator old = *this; + this->operator++(); + return old; +} + + +inline Foam::phasePair::const_iterator Foam::phasePair::cbegin() const +{ + return const_iterator(*this); +} + + +inline Foam::phasePair::const_iterator Foam::phasePair::cend() const +{ + return const_iterator(*this, 2); +} + + +inline Foam::phasePair::const_iterator Foam::phasePair::begin() const +{ + return const_iterator(*this); +} + + +inline Foam::phasePair::const_iterator Foam::phasePair::end() const +{ + return const_iterator(*this, 2); +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePairKey/phasePairKey.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePairKey/phasePairKey.C index b138653e49..0d52f5fffd 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePairKey/phasePairKey.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phasePair/phasePairKey/phasePairKey.C @@ -127,7 +127,7 @@ Foam::Istream& Foam::operator>>(Istream& is, phasePairKey& key) { key.ordered_ = false; } - else if(temp[1] == "in") + else if (temp[1] == "in") { key.ordered_ = true; } diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 6ae600c92f..581beb9751 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -40,12 +40,12 @@ void Foam::diameterModels::populationBalanceModel::registerVelocityGroups() { forAll(fluid_.phases(), phasei) { - if(isA(fluid_.phases()[phasei].dPtr()())) + if (isA(fluid_.phases()[phasei].dPtr()())) { const velocityGroup& velGroup = refCast(fluid_.phases()[phasei].dPtr()()); - if(velGroup.popBalName() == this->name()) + if (velGroup.popBalName() == this->name()) { velocityGroups_.append(&const_cast(velGroup)); diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H index 0fdad237bd..29439be542 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H @@ -229,7 +229,7 @@ while (pimple.correct()) if (pair.contains(phase)) { - const phaseModel& otherPhase = pair.other(phase); + const phaseModel& otherPhase = pair.otherPhase(phase); phiHbyAs[phasei] += fvc::interpolate(rAUs[phasei]*K) diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H index d7a9c06bc1..4b3f036fe9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H @@ -62,7 +62,7 @@ PtrList UEqns(phases.size()); if (pair.contains(phase)) { - const phaseModel& otherPhase = pair.other(phase); + const phaseModel& otherPhase = pair.otherPhase(phase); UEqns[phasei] += Vm diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H index b211c71f02..cfdbfd5da9 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H @@ -197,7 +197,7 @@ while (pimple.correct()) if (pair.contains(phase)) { phiHbyAs[phasei] += - rAUfs[phasei]*Kf*MRF.absolute(pair.other(phase).phi()); + rAUfs[phasei]*Kf*MRF.absolute(pair.otherPhase(phase).phi()); } } @@ -219,7 +219,7 @@ while (pimple.correct()) rAUfs[phasei] *( Vmf*byDt(MRF.absolute(phase.phi().oldTime())) - + Vmf*ddtPhis[pair.other(phase).index()] + + Vmf*ddtPhis[pair.otherPhase(phase).index()] ); } }