meshToMeshTemplates.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 3 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, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "fvMesh.H"
#include "volFields.H"
#include "directFvPatchFieldMapper.H"
#include "calculatedFvPatchField.H"
#include "fvcGrad.H"
#include "distributedWeightedFvPatchFieldMapper.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    //- Helper class for list
    template<class Type>
    class ListPlusEqOp
    {
        public:
        void operator()(List<Type>& x, const List<Type> y) const
        {
            if (y.size())
            {
                if (x.size())
                {
                    label sz = x.size();
                    x.setSize(sz + y.size());
                    forAll(y, i)
                    {
                        x[sz++] = y[i];
                    }
                }
                else
                {
                    x = y;
                }
            }
        }
    };
}
 
 
// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
template<class Type>
void Foam::meshToMesh::add
(
    UList<Type>& fld,
    const label offset
) const
{
    forAll(fld, i)
    {
        fld[i] += offset;
    }
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapSrcToTgt
(
    const UList<Type>& srcField,
    const CombineOp& cop,
    List<Type>& result
) const
{
    if (result.size() != tgtToSrcCellAddr_.size())
    {
        FatalErrorInFunction
            << "Supplied field size is not equal to target mesh size" << nl
            << "    source mesh    = " << srcToTgtCellAddr_.size() << nl
            << "    target mesh    = " << tgtToSrcCellAddr_.size() << nl
            << "    supplied field = " << result.size()
            << abort(FatalError);
    }
 
    multiplyWeightedOp<Type, CombineOp> cbop(cop);
 
    if (singleMeshProc_ == -1)
    {
        const mapDistribute& map = srcMapPtr_();
 
        List<Type> work(srcField);
        map.distribute(work);
 
        forAll(result, cellI)
        {
            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
 
            if (srcAddress.size())
            {
                result[cellI] *= (1.0 - sum(srcWeight));
                forAll(srcAddress, i)
                {
                    label srcI = srcAddress[i];
                    scalar w = srcWeight[i];
                    cbop(result[cellI], cellI, work[srcI], w);
                }
            }
        }
    }
    else
    {
        forAll(result, cellI)
        {
            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
 
            if (srcAddress.size())
            {
                result[cellI] *= (1.0 - sum(srcWeight));
                forAll(srcAddress, i)
                {
                    label srcI = srcAddress[i];
                    scalar w = srcWeight[i];
                    cbop(result[cellI], cellI, srcField[srcI], w);
                }
            }
        }
    }
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapSrcToTgt
(
    const UList<Type>& srcField,
    const UList<typename outerProduct<vector, Type>::type>& srcGradField,
    const CombineOp& cop,
    List<Type>& result
) const
{
    if (result.size() != tgtToSrcCellAddr_.size())
    {
        FatalErrorInFunction
            << "Supplied field size is not equal to target mesh size" << nl
            << "    source mesh    = " << srcToTgtCellAddr_.size() << nl
            << "    target mesh    = " << tgtToSrcCellAddr_.size() << nl
            << "    supplied field = " << result.size()
            << abort(FatalError);
    }
 
    multiplyWeightedOp<Type, CombineOp> cbop(cop);
 
    if (singleMeshProc_ == -1)
    {
        if (returnReduce(tgtToSrcCellVec_.size(), sumOp<label>()) == 0)
        {
            // No correction vectors calculated. Fall back to first order.
            mapSrcToTgt(srcField, cop, result);
            return;
        }
 
        const mapDistribute& map = srcMapPtr_();
 
        List<Type> work(srcField);
        map.distribute(work);
 
        List<typename outerProduct<vector, Type>::type> workGrad
        (
            srcGradField
        );
        map.distribute(workGrad);
 
        forAll(result, cellI)
        {
            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
            const pointList& srcVec = tgtToSrcCellVec_[cellI];
 
            if (srcAddress.size())
            {
                result[cellI] *= (1.0 - sum(srcWeight));
                forAll(srcAddress, i)
                {
                    label srcI = srcAddress[i];
                    scalar w = srcWeight[i];
                    const vector& v = srcVec[i];
                    const Type srcVal = work[srcI]+(workGrad[srcI]&v);
                    cbop(result[cellI], cellI, srcVal, w);
                }
            }
        }
    }
    else
    {
        if (tgtToSrcCellVec_.empty())
        {
            // No correction vectors calculated. Fall back to first order.
            mapSrcToTgt(srcField, cop, result);
            return;
        }
 
        forAll(result, cellI)
        {
            const labelList& srcAddress = tgtToSrcCellAddr_[cellI];
            const scalarList& srcWeight = tgtToSrcCellWght_[cellI];
            const pointList& srcVec = tgtToSrcCellVec_[cellI];
 
            if (srcAddress.size())
            {
                // Do non-conservative interpolation
                result[cellI] *= (1.0 - sum(srcWeight));
                forAll(srcAddress, i)
                {
                    label srcI = srcAddress[i];
                    scalar w = srcWeight[i];
                    const vector& v = srcVec[i];
                    const Type srcVal = srcField[srcI]+(srcGradField[srcI]&v);
                    cbop(result[cellI], cellI, srcVal, w);
                }
            }
        }
    }
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapSrcToTgt
(
    const Field<Type>& srcField,
    const CombineOp& cop
) const
{
    tmp<Field<Type> > tresult
    (
        new Field<Type>
        (
            tgtToSrcCellAddr_.size(),
            pTraits<Type>::zero
        )
    );
 
    mapSrcToTgt(srcField, cop, tresult());
 
    return tresult;
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapSrcToTgt
(
    const tmp<Field<Type> >& tsrcField,
    const CombineOp& cop
) const
{
    return mapSrcToTgt(tsrcField(), cop);
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapSrcToTgt
(
    const Field<Type>& srcField
) const
{
    return mapSrcToTgt(srcField, plusEqOp<Type>());
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapSrcToTgt
(
    const tmp<Field<Type> >& tsrcField
) const
{
    return mapSrcToTgt(tsrcField());
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapTgtToSrc
(
    const UList<Type>& tgtField,
    const CombineOp& cop,
    List<Type>& result
) const
{
    if (result.size() != srcToTgtCellAddr_.size())
    {
        FatalErrorInFunction
            << "Supplied field size is not equal to source mesh size" << nl
            << "    source mesh    = " << srcToTgtCellAddr_.size() << nl
            << "    target mesh    = " << tgtToSrcCellAddr_.size() << nl
            << "    supplied field = " << result.size()
            << abort(FatalError);
    }
 
    multiplyWeightedOp<Type, CombineOp> cbop(cop);
 
    if (singleMeshProc_ == -1)
    {
        const mapDistribute& map = tgtMapPtr_();
 
        List<Type> work(tgtField);
        map.distribute(work);
 
        forAll(result, cellI)
        {
            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
 
            if (tgtAddress.size())
            {
                result[cellI] *= (1.0 - sum(tgtWeight));
                forAll(tgtAddress, i)
                {
                    label tgtI = tgtAddress[i];
                    scalar w = tgtWeight[i];
                    cbop(result[cellI], cellI, work[tgtI], w);
                }
            }
        }
    }
    else
    {
        forAll(result, cellI)
        {
            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
 
            if (tgtAddress.size())
            {
                result[cellI] *= (1.0 - sum(tgtWeight));
                forAll(tgtAddress, i)
                {
                    label tgtI = tgtAddress[i];
                    scalar w = tgtWeight[i];
                    cbop(result[cellI], cellI, tgtField[tgtI], w);
                }
            }
        }
    }
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapTgtToSrc
(
    const UList<Type>& tgtField,
    const UList<typename outerProduct<vector, Type>::type>& tgtGradField,
    const CombineOp& cop,
    List<Type>& result
) const
{
    if (result.size() != srcToTgtCellAddr_.size())
    {
        FatalErrorInFunction
            << "Supplied field size is not equal to source mesh size" << nl
            << "    source mesh    = " << srcToTgtCellAddr_.size() << nl
            << "    target mesh    = " << tgtToSrcCellAddr_.size() << nl
            << "    supplied field = " << result.size()
            << abort(FatalError);
    }
 
    multiplyWeightedOp<Type, CombineOp> cbop(cop);
 
    if (singleMeshProc_ == -1)
    {
        if (returnReduce(srcToTgtCellVec_.size(), sumOp<label>()) == 0)
        {
            // No correction vectors calculated. Fall back to first order.
            mapTgtToSrc(tgtField, cop, result);
            return;
        }
 
        const mapDistribute& map = tgtMapPtr_();
 
        List<Type> work(tgtField);
        map.distribute(work);
 
        List<typename outerProduct<vector, Type>::type> workGrad
        (
            tgtGradField
        );
        map.distribute(workGrad);
 
        forAll(result, cellI)
        {
            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
            const pointList& tgtVec = srcToTgtCellVec_[cellI];
 
            if (tgtAddress.size())
            {
                result[cellI] *= (1.0 - sum(tgtWeight));
                forAll(tgtAddress, i)
                {
                    label tgtI = tgtAddress[i];
                    scalar w = tgtWeight[i];
                    const vector& v = tgtVec[i];
                    const Type tgtVal = work[tgtI]+(workGrad[tgtI]&v);
                    cbop(result[cellI], cellI, tgtVal, w);
                }
            }
        }
    }
    else
    {
        forAll(result, cellI)
        {
            const labelList& tgtAddress = srcToTgtCellAddr_[cellI];
            const scalarList& tgtWeight = srcToTgtCellWght_[cellI];
            const pointList& tgtVec = srcToTgtCellVec_[cellI];
 
            if (tgtAddress.size())
            {
                result[cellI] *= (1.0 - sum(tgtWeight));
                forAll(tgtAddress, i)
                {
                    label tgtI = tgtAddress[i];
                    scalar w = tgtWeight[i];
                    const vector& v = tgtVec[i];
                    const Type tgtVal = tgtField[tgtI]+(tgtGradField[tgtI]&v);
                    cbop(result[cellI], cellI, tgtVal, w);
                }
            }
        }
    }
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapTgtToSrc
(
    const Field<Type>& tgtField,
    const CombineOp& cop
) const
{
    tmp<Field<Type> > tresult
    (
        new Field<Type>
        (
            srcToTgtCellAddr_.size(),
            pTraits<Type>::zero
        )
    );
 
    mapTgtToSrc(tgtField, cop, tresult());
 
    return tresult;
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapTgtToSrc
(
    const tmp<Field<Type> >& ttgtField,
    const CombineOp& cop
) const
{
    return mapTgtToSrc(ttgtField(), cop);
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapTgtToSrc
(
    const Field<Type>& tgtField
) const
{
    return mapTgtToSrc(tgtField, plusEqOp<Type>());
}
 
 
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::meshToMesh::mapTgtToSrc
(
    const tmp<Field<Type> >& ttgtField
) const
{
    return mapTgtToSrc(ttgtField(), plusEqOp<Type>());
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapInternalSrcToTgt
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    GeometricField<Type, fvPatchField, volMesh>& result,
    const bool secondOrder
) const
{
    if (secondOrder && returnReduce(tgtToSrcCellVec_.size(), sumOp<label>()))
    {
        mapSrcToTgt
        (
            field,
            fvc::grad(field)().internalField(),
            cop,
            result.internalField()
        );
    }
    else
    {
        mapSrcToTgt(field, cop, result.internalField());
    }
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapAndOpSrcToTgt
(
    const AMIPatchToPatchInterpolation& AMI,
    const Field<Type>& srcField,
    Field<Type>& tgtField,
    const CombineOp& cop
) const
{
    tgtField = pTraits<Type>::zero;
 
    AMI.interpolateToTarget
    (
        srcField,
        multiplyWeightedOp<Type, CombineOp>(cop),
        tgtField,
        UList<Type>::null()
    );
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapSrcToTgt
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    GeometricField<Type, fvPatchField, volMesh>& result,
    const bool secondOrder
) const
{
    mapInternalSrcToTgt(field, cop, result, secondOrder);
 
 
    const PtrList<AMIPatchToPatchInterpolation>& AMIList = patchAMIs();
 
    forAll(AMIList, i)
    {
        label srcPatchI = srcPatchID_[i];
        label tgtPatchI = tgtPatchID_[i];
 
        const fvPatchField<Type>& srcField = field.boundaryField()[srcPatchI];
        fvPatchField<Type>& tgtField = result.boundaryField()[tgtPatchI];
 
 
        // Clone and map (since rmap does not do general mapping)
        tmp<fvPatchField<Type> > tnewTgt
        (
            fvPatchField<Type>::New
            (
                srcField,
                tgtField.patch(),
                result.dimensionedInternalField(),
                distributedWeightedFvPatchFieldMapper
                (
                    AMIList[i].singlePatchProc(),
                    (
                        AMIList[i].singlePatchProc() == -1
                      ? &AMIList[i].srcMap()
                      : NULL
                    ),
                    AMIList[i].tgtAddress(),
                    AMIList[i].tgtWeights()
                )
            )
        );
 
        // Transfer all mapped quantities (value and e.g. gradient) onto
        // tgtField
        tgtField.rmap(tnewTgt(), identity(tgtField.size()));
 
 
        // Override value to account for CombineOp (note: is dummy template
        // specialisation for plusEqOp)
        mapAndOpSrcToTgt(AMIList[i], srcField, tgtField, cop);
    }
 
    forAll(cuttingPatches_, i)
    {
        label patchI = cuttingPatches_[i];
        fvPatchField<Type>& pf = result.boundaryField()[patchI];
        pf == pf.patchInternalField();
    }
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapSrcToTgt
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    const bool secondOrder
) const
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
    const fvMesh& tgtMesh = static_cast<const fvMesh&>(tgtRegion_);
 
    const fvBoundaryMesh& tgtBm = tgtMesh.boundary();
    const typename fieldType::GeometricBoundaryField& srcBfld =
        field.boundaryField();
 
    PtrList<fvPatchField<Type> > tgtPatchFields(tgtBm.size());
 
    // constuct tgt boundary patch types as copy of 'field' boundary types
    // note: this will provide place holders for fields with additional
    // entries, but these values will need to be reset
    forAll(tgtPatchID_, i)
    {
        label srcPatchI = srcPatchID_[i];
        label tgtPatchI = tgtPatchID_[i];
 
        if (!tgtPatchFields.set(tgtPatchI))
        {
            tgtPatchFields.set
            (
                tgtPatchI,
                fvPatchField<Type>::New
                (
                    srcBfld[srcPatchI],
                    tgtMesh.boundary()[tgtPatchI],
                    DimensionedField<Type, volMesh>::null(),
                    directFvPatchFieldMapper
                    (
                        labelList(tgtMesh.boundary()[tgtPatchI].size(), -1)
                    )
                )
            );
        }
    }
 
    // Any unset tgtPatchFields become calculated
    forAll(tgtPatchFields, tgtPatchI)
    {
        if (!tgtPatchFields.set(tgtPatchI))
        {
            // Note: use factory New method instead of direct generation of
            //       calculated so we keep constraints
            tgtPatchFields.set
            (
                tgtPatchI,
                fvPatchField<Type>::New
                (
                    calculatedFvPatchField<Type>::typeName,
                    tgtMesh.boundary()[tgtPatchI],
                    DimensionedField<Type, volMesh>::null()
                )
            );
        }
    }
 
    tmp<fieldType> tresult
    (
        new fieldType
        (
            IOobject
            (
                type() + ":interpolate(" + field.name() + ")",
                tgtMesh.time().timeName(),
                tgtMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            tgtMesh,
            field.dimensions(),
            Field<Type>(tgtMesh.nCells(), pTraits<Type>::zero),
            tgtPatchFields
        )
    );
 
    mapSrcToTgt(field, cop, tresult(), secondOrder);
 
    return tresult;
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapSrcToTgt
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
    const CombineOp& cop,
    const bool secondOrder
) const
{
    return mapSrcToTgt(tfield(), cop, secondOrder);
}
 
 
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapSrcToTgt
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const bool secondOrder
) const
{
    return mapSrcToTgt(field, plusEqOp<Type>(), secondOrder);
}
 
 
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapSrcToTgt
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
    const bool secondOrder
) const
{
    return mapSrcToTgt(tfield(), plusEqOp<Type>(), secondOrder);
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapInternalTgtToSrc
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    GeometricField<Type, fvPatchField, volMesh>& result,
    const bool secondOrder
) const
{
    if (secondOrder && returnReduce(srcToTgtCellVec_.size(), sumOp<label>()))
    {
        mapTgtToSrc
        (
            field,
            fvc::grad(field)().internalField(),
            cop,
            result.internalField()
        );
    }
    else
    {
        mapTgtToSrc(field, cop, result.internalField());
    }
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapAndOpTgtToSrc
(
    const AMIPatchToPatchInterpolation& AMI,
    Field<Type>& srcField,
    const Field<Type>& tgtField,
    const CombineOp& cop
) const
{
    srcField = pTraits<Type>::zero;
 
    AMI.interpolateToSource
    (
        tgtField,
        multiplyWeightedOp<Type, CombineOp>(cop),
        srcField,
        UList<Type>::null()
    );
}
 
 
template<class Type, class CombineOp>
void Foam::meshToMesh::mapTgtToSrc
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    GeometricField<Type, fvPatchField, volMesh>& result,
    const bool secondOrder
) const
{
    mapInternalTgtToSrc(field, cop, result, secondOrder);
 
 
    const PtrList<AMIPatchToPatchInterpolation>& AMIList = patchAMIs();
 
    forAll(AMIList, i)
    {
        label srcPatchI = srcPatchID_[i];
        label tgtPatchI = tgtPatchID_[i];
 
        fvPatchField<Type>& srcField = result.boundaryField()[srcPatchI];
        const fvPatchField<Type>& tgtField = field.boundaryField()[tgtPatchI];
 
 
        // Clone and map (since rmap does not do general mapping)
        tmp<fvPatchField<Type> > tnewSrc
        (
            fvPatchField<Type>::New
            (
                tgtField,
                srcField.patch(),
                result.dimensionedInternalField(),
                distributedWeightedFvPatchFieldMapper
                (
                    AMIList[i].singlePatchProc(),
                    (
                        AMIList[i].singlePatchProc() == -1
                      ? &AMIList[i].tgtMap()
                      : NULL
                    ),
                    AMIList[i].srcAddress(),
                    AMIList[i].srcWeights()
                )
            )
        );
 
        // Transfer all mapped quantities (value and e.g. gradient) onto
        // tgtField
        srcField.rmap(tnewSrc(), identity(srcField.size()));
 
 
        // Override value to account for CombineOp (could be dummy for
        // plusEqOp)
        mapAndOpTgtToSrc(AMIList[i], srcField, tgtField, cop);
    }
 
    forAll(cuttingPatches_, i)
    {
        label patchI = cuttingPatches_[i];
        fvPatchField<Type>& pf = result.boundaryField()[patchI];
        pf == pf.patchInternalField();
    }
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapTgtToSrc
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const CombineOp& cop,
    const bool secondOrder
) const
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
    const fvMesh& srcMesh = static_cast<const fvMesh&>(srcRegion_);
 
    const fvBoundaryMesh& srcBm = srcMesh.boundary();
    const typename fieldType::GeometricBoundaryField& tgtBfld =
        field.boundaryField();
 
    PtrList<fvPatchField<Type> > srcPatchFields(srcBm.size());
 
    // constuct src boundary patch types as copy of 'field' boundary types
    // note: this will provide place holders for fields with additional
    // entries, but these values will need to be reset
    forAll(srcPatchID_, i)
    {
        label srcPatchI = srcPatchID_[i];
        label tgtPatchI = tgtPatchID_[i];
 
        if (!srcPatchFields.set(tgtPatchI))
        {
            srcPatchFields.set
            (
                srcPatchI,
                fvPatchField<Type>::New
                (
                    tgtBfld[srcPatchI],
                    srcMesh.boundary()[tgtPatchI],
                    DimensionedField<Type, volMesh>::null(),
                    directFvPatchFieldMapper
                    (
                        labelList(srcMesh.boundary()[srcPatchI].size(), -1)
                    )
                )
            );
        }
    }
 
    // Any unset srcPatchFields become calculated
    forAll(srcPatchFields, srcPatchI)
    {
        if (!srcPatchFields.set(srcPatchI))
        {
            // Note: use factory New method instead of direct generation of
            //       calculated so we keep constraints
            srcPatchFields.set
            (
                srcPatchI,
                fvPatchField<Type>::New
                (
                    calculatedFvPatchField<Type>::typeName,
                    srcMesh.boundary()[srcPatchI],
                    DimensionedField<Type, volMesh>::null()
                )
            );
        }
    }
 
    tmp<fieldType> tresult
    (
        new fieldType
        (
            IOobject
            (
                type() + ":interpolate(" + field.name() + ")",
                srcMesh.time().timeName(),
                srcMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            srcMesh,
            field.dimensions(),
            Field<Type>(srcMesh.nCells(), pTraits<Type>::zero),
            srcPatchFields
        )
    );
 
    mapTgtToSrc(field, cop, tresult(), secondOrder);
 
    return tresult;
}
 
 
template<class Type, class CombineOp>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapTgtToSrc
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
    const CombineOp& cop,
    const bool secondOrder
) const
{
    return mapTgtToSrc(tfield(), cop, secondOrder);
}
 
 
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapTgtToSrc
(
    const GeometricField<Type, fvPatchField, volMesh>& field,
    const bool secondOrder
) const
{
    return mapTgtToSrc(field, plusEqOp<Type>(), secondOrder);
}
 
 
template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::meshToMesh::mapTgtToSrc
(
    const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfield,
    const bool secondOrder
) const
{
    return mapTgtToSrc(tfield(), plusEqOp<Type>(), secondOrder);
}
 
 
// ************************************************************************* //